Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <pup.h>
8 :
9 : #include "DataStructures/CachedTempBuffer.hpp"
10 : #include "DataStructures/ComplexDataVector.hpp"
11 : #include "DataStructures/DataBox/Prefixes.hpp"
12 : #include "DataStructures/DataVector.hpp"
13 : #include "DataStructures/TaggedTuple.hpp"
14 : #include "DataStructures/Tensor/Tensor.hpp"
15 : #include "Elliptic/Systems/Poisson/Tags.hpp"
16 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
17 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
18 : #include "Options/String.hpp"
19 : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
20 : #include "Utilities/ContainerHelpers.hpp"
21 : #include "Utilities/Gsl.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 : namespace Poisson::Solutions {
25 :
26 : namespace detail {
27 : template <typename DataType, size_t Dim>
28 : struct LorentzianVariables {
29 : using Cache = CachedTempBuffer<
30 : Tags::Field<DataType>,
31 : ::Tags::deriv<Tags::Field<DataType>, tmpl::size_t<Dim>, Frame::Inertial>,
32 : ::Tags::Flux<Tags::Field<DataType>, tmpl::size_t<Dim>, Frame::Inertial>,
33 : ::Tags::FixedSource<Tags::Field<DataType>>>;
34 :
35 : // NOLINTNEXTLINE(cppcoreguidelines-avoid-const-or-ref-data-members)
36 : const tnsr::I<DataVector, Dim>& x;
37 : double constant;
38 : double complex_phase;
39 :
40 : void operator()(gsl::not_null<Scalar<DataType>*> field,
41 : gsl::not_null<Cache*> cache,
42 : Tags::Field<DataType> /*meta*/) const;
43 : void operator()(gsl::not_null<tnsr::i<DataType, Dim>*> field_gradient,
44 : gsl::not_null<Cache*> cache,
45 : ::Tags::deriv<Tags::Field<DataType>, tmpl::size_t<Dim>,
46 : Frame::Inertial> /*meta*/) const;
47 : void operator()(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
48 : gsl::not_null<Cache*> cache,
49 : ::Tags::Flux<Tags::Field<DataType>, tmpl::size_t<Dim>,
50 : Frame::Inertial> /*meta*/) const;
51 : void operator()(gsl::not_null<Scalar<DataType>*> fixed_source_for_field,
52 : gsl::not_null<Cache*> cache,
53 : ::Tags::FixedSource<Tags::Field<DataType>> /*meta*/) const;
54 : };
55 : } // namespace detail
56 :
57 : /*!
58 : * \brief A Lorentzian solution to the Poisson equation
59 : *
60 : * \details This implements the Lorentzian solution
61 : * \f$u(\boldsymbol{x})=\left(1+r^2\right)^{-\frac{1}{2}}\f$ to the
62 : * three-dimensional Poisson equation
63 : * \f$-\Delta u(\boldsymbol{x})=f(\boldsymbol{x})\f$, where
64 : * \f$r^2=x^2+y^2+z^2\f$. The corresponding source is
65 : * \f$f(\boldsymbol{x})=3\left(1+r^2\right)^{-\frac{5}{2}}\f$.
66 : *
67 : * If `DataType` is `ComplexDataVector`, the solution is multiplied by
68 : * `exp(i * complex_phase)` to rotate it in the complex plane. This allows to
69 : * use this solution for the complex Poisson equation.
70 : *
71 : * \note Corresponding 1D and 2D solutions are not implemented yet.
72 : */
73 : template <size_t Dim, typename DataType = DataVector>
74 1 : class Lorentzian : public elliptic::analytic_data::AnalyticSolution {
75 : static_assert(
76 : Dim == 3,
77 : "This solution is currently implemented in 3 spatial dimensions only");
78 :
79 : public:
80 0 : struct PlusConstant {
81 0 : using type = double;
82 0 : static constexpr Options::String help{"Constant added to the solution."};
83 : };
84 :
85 0 : struct ComplexPhase {
86 0 : using type = double;
87 0 : static constexpr Options::String help{
88 : "Phase 'phi' of a complex exponential 'exp(i phi)' that rotates the "
89 : "solution in the complex plane."};
90 : };
91 :
92 0 : using options = tmpl::flatten<tmpl::list<
93 : PlusConstant,
94 : tmpl::conditional_t<std::is_same_v<DataType, ComplexDataVector>,
95 : ComplexPhase, tmpl::list<>>>>;
96 0 : static constexpr Options::String help{
97 : "A Lorentzian solution to the Poisson equation."};
98 :
99 0 : Lorentzian() = default;
100 0 : Lorentzian(const Lorentzian&) = default;
101 0 : Lorentzian& operator=(const Lorentzian&) = default;
102 0 : Lorentzian(Lorentzian&&) = default;
103 0 : Lorentzian& operator=(Lorentzian&&) = default;
104 0 : ~Lorentzian() override = default;
105 :
106 0 : explicit Lorentzian(const double constant, const double complex_phase = 0.)
107 : : constant_(constant), complex_phase_(complex_phase) {
108 : ASSERT((std::is_same_v<DataType, ComplexDataVector> or complex_phase == 0.),
109 : "The complex phase is only supported for ComplexDataVector.");
110 : }
111 :
112 0 : double constant() const { return constant_; }
113 0 : double complex_phase() const { return complex_phase_; }
114 :
115 0 : std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
116 : const override {
117 : return std::make_unique<Lorentzian>(*this);
118 : }
119 :
120 : /// \cond
121 : explicit Lorentzian(CkMigrateMessage* m)
122 : : elliptic::analytic_data::AnalyticSolution(m) {}
123 : using PUP::able::register_constructor;
124 : WRAPPED_PUPable_decl_template(Lorentzian); // NOLINT
125 : /// \endcond
126 :
127 : template <typename... RequestedTags>
128 0 : tuples::TaggedTuple<RequestedTags...> variables(
129 : const tnsr::I<DataVector, Dim>& x,
130 : tmpl::list<RequestedTags...> /*meta*/) const {
131 : using VarsComputer = detail::LorentzianVariables<DataType, Dim>;
132 : typename VarsComputer::Cache cache{get_size(*x.begin())};
133 : const VarsComputer computer{x, constant_, complex_phase_};
134 : return {cache.get_var(computer, RequestedTags{})...};
135 : }
136 :
137 : template <typename... RequestedTags>
138 0 : tuples::TaggedTuple<RequestedTags...> variables(
139 : const tnsr::I<DataVector, Dim>& x, const Mesh<Dim>& /*mesh*/,
140 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
141 : Frame::Inertial>& /*inv_jacobian*/,
142 : tmpl::list<RequestedTags...> meta) const {
143 : return variables(x, meta);
144 : }
145 :
146 0 : void pup(PUP::er& p) override {
147 : elliptic::analytic_data::AnalyticSolution::pup(p);
148 : p | constant_;
149 : p | complex_phase_;
150 : }
151 :
152 : private:
153 0 : double constant_ = std::numeric_limits<double>::signaling_NaN();
154 0 : double complex_phase_ = std::numeric_limits<double>::signaling_NaN();
155 : };
156 :
157 : /// \cond
158 : template <size_t Dim, typename DataType>
159 : PUP::able::PUP_ID Lorentzian<Dim, DataType>::my_PUP_ID = 0; // NOLINT
160 : /// \endcond
161 :
162 : template <size_t Dim, typename DataType>
163 0 : bool operator==(const Lorentzian<Dim, DataType>& lhs,
164 : const Lorentzian<Dim, DataType>& rhs) {
165 : return lhs.constant() == rhs.constant() and
166 : lhs.complex_phase() == rhs.complex_phase();
167 : }
168 :
169 : template <size_t Dim, typename DataType>
170 0 : bool operator!=(const Lorentzian<Dim, DataType>& lhs,
171 : const Lorentzian<Dim, DataType>& rhs) {
172 : return not(lhs == rhs);
173 : }
174 :
175 : } // namespace Poisson::Solutions
|
