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/DataBox/Prefixes.hpp"
11 : #include "DataStructures/DataVector.hpp"
12 : #include "DataStructures/TaggedTuple.hpp"
13 : #include "DataStructures/Tensor/Tensor.hpp"
14 : #include "Elliptic/Systems/Poisson/Tags.hpp"
15 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
16 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
17 : #include "Options/String.hpp"
18 : #include "PointwiseFunctions/InitialDataUtilities/AnalyticSolution.hpp"
19 : #include "Utilities/ContainerHelpers.hpp"
20 : #include "Utilities/Gsl.hpp"
21 : #include "Utilities/TMPL.hpp"
22 :
23 : namespace Poisson::Solutions {
24 :
25 : namespace detail {
26 : template <typename DataType, size_t Dim>
27 : struct MoustacheVariables {
28 : using Cache = CachedTempBuffer<
29 : Tags::Field<DataType>,
30 : ::Tags::deriv<Tags::Field<DataType>, tmpl::size_t<Dim>, Frame::Inertial>,
31 : ::Tags::Flux<Tags::Field<DataType>, tmpl::size_t<Dim>, Frame::Inertial>,
32 : ::Tags::FixedSource<Tags::Field<DataType>>>;
33 :
34 : const tnsr::I<DataType, Dim>& x;
35 :
36 : void operator()(gsl::not_null<Scalar<DataType>*> field,
37 : gsl::not_null<Cache*> cache,
38 : Tags::Field<DataType> /*meta*/) const;
39 : void operator()(gsl::not_null<tnsr::i<DataType, Dim>*> field_gradient,
40 : gsl::not_null<Cache*> cache,
41 : ::Tags::deriv<Tags::Field<DataType>, tmpl::size_t<Dim>,
42 : Frame::Inertial> /*meta*/) const;
43 : void operator()(gsl::not_null<tnsr::I<DataType, Dim>*> flux_for_field,
44 : gsl::not_null<Cache*> cache,
45 : ::Tags::Flux<Tags::Field<DataType>, tmpl::size_t<Dim>,
46 : Frame::Inertial> /*meta*/) const;
47 : void operator()(gsl::not_null<Scalar<DataType>*> fixed_source_for_field,
48 : gsl::not_null<Cache*> cache,
49 : ::Tags::FixedSource<Tags::Field<DataType>> /*meta*/) const;
50 : };
51 : } // namespace detail
52 :
53 : /*!
54 : * \brief A solution to the Poisson equation with a discontinuous first
55 : * derivative.
56 : *
57 : * \details This implements the solution \f$u(x,y)=x\left(1-x\right)
58 : * y\left(1-y\right)\left(\left(x-\frac{1}{2}\right)^2+\left(y-
59 : * \frac{1}{2}\right)^2\right)^\frac{3}{2}\f$ to the Poisson equation
60 : * in two dimensions, and
61 : * \f$u(x)=x\left(1-x\right)\left|x-\frac{1}{2}\right|^3\f$ in one dimension.
62 : * Their boundary conditions vanish on the square \f$[0,1]^2\f$ or interval
63 : * \f$[0,1]\f$, respectively.
64 : *
65 : * The corresponding source \f$f=-\Delta u\f$ has a discontinuous first
66 : * derivative at \f$\frac{1}{2}\f$. This accomplishes two things:
67 : *
68 : * - It makes it useful to test the convergence behaviour of our elliptic DG
69 : * solver.
70 : * - It makes it look like a moustache (at least in 1D).
71 : *
72 : * This solution is taken from \cite Stamm2010.
73 : */
74 : template <size_t Dim>
75 1 : class Moustache : public elliptic::analytic_data::AnalyticSolution {
76 : public:
77 0 : using options = tmpl::list<>;
78 0 : static constexpr Options::String help{
79 : "A solution with a discontinuous first derivative of its source at 1/2 "
80 : "that also happens to look like a moustache. It vanishes at zero and one "
81 : "in each dimension"};
82 :
83 0 : Moustache() = default;
84 0 : Moustache(const Moustache&) = default;
85 0 : Moustache& operator=(const Moustache&) = default;
86 0 : Moustache(Moustache&&) = default;
87 0 : Moustache& operator=(Moustache&&) = default;
88 0 : ~Moustache() override = default;
89 0 : std::unique_ptr<elliptic::analytic_data::AnalyticSolution> get_clone()
90 : const override {
91 : return std::make_unique<Moustache>(*this);
92 : }
93 :
94 : /// \cond
95 : explicit Moustache(CkMigrateMessage* m)
96 : : elliptic::analytic_data::AnalyticSolution(m) {}
97 : using PUP::able::register_constructor;
98 : WRAPPED_PUPable_decl_template(Moustache); // NOLINT
99 : /// \endcond
100 :
101 : template <typename DataType, typename... RequestedTags>
102 0 : tuples::TaggedTuple<RequestedTags...> variables(
103 : const tnsr::I<DataType, Dim>& x,
104 : tmpl::list<RequestedTags...> /*meta*/) const {
105 : using VarsComputer = detail::MoustacheVariables<DataType, Dim>;
106 : typename VarsComputer::Cache cache{get_size(*x.begin())};
107 : const VarsComputer computer{x};
108 : return {cache.get_var(computer, RequestedTags{})...};
109 : }
110 :
111 : template <typename... RequestedTags>
112 0 : tuples::TaggedTuple<RequestedTags...> variables(
113 : const tnsr::I<DataVector, Dim>& x, const Mesh<Dim>& /*mesh*/,
114 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
115 : Frame::Inertial>& /*inv_jacobian*/,
116 : tmpl::list<RequestedTags...> meta) const {
117 : return variables(x, meta);
118 : }
119 : };
120 :
121 : /// \cond
122 : template <size_t Dim>
123 : PUP::able::PUP_ID Moustache<Dim>::my_PUP_ID = 0; // NOLINT
124 : /// \endcond
125 :
126 : template <size_t Dim>
127 0 : constexpr bool operator==(const Moustache<Dim>& /*lhs*/,
128 : const Moustache<Dim>& /*rhs*/) {
129 : return true;
130 : }
131 :
132 : template <size_t Dim>
133 0 : constexpr bool operator!=(const Moustache<Dim>& lhs,
134 : const Moustache<Dim>& rhs) {
135 : return not(lhs == rhs);
136 : }
137 :
138 : } // namespace Poisson::Solutions
|
