1 #include <iostream>
2 #include <stdexcept>
3
4 #include "metaphysicl_config.h"
5
6 // If we have MASA we test ourselves against an MMS solution; if not
7 // we just test that this compiles.
8 #ifdef METAPHYSICL_HAVE_MASA
9 # include <masa.h>
10 using namespace MASA;
11 #else
12 # define masa_get_param(arg) 1;
13 #endif // METAPHYSICL_HAVE_MASA
14
15 using namespace MetaPhysicL;
16
17 template <typename Vector>
18 static double evaluate_q (const Vector& xyz, const int);
19
20 using namespace MetaPhysicL;
21
main(void)22 int main(void)
23 {
24 int N = 10; // mesh pts. in x and y
25 double s2u,s2v,s2e,s2p;
26
27 #ifdef METAPHYSICL_HAVE_MASA
28 int err = 0;
29 double su,sv,sp,se;
30 double pnorm, unorm, vnorm, enorm;
31 double pnorm_max = 0., unorm_max = 0., vnorm_max = 0., enorm_max = 0.;
32 double prnorm_max = 0., urnorm_max = 0., vrnorm_max = 0., ernorm_max = 0.;
33 #endif
34
35 typedef VectorUnitVector<NDIM,0,RawScalar>::type XVector;
36 XVector xvec = VectorUnitVector<NDIM,0,RawScalar>::value();
37
38 typedef VectorUnitVector<NDIM,1,RawScalar>::type YVector;
39 YVector yvec = VectorUnitVector<NDIM,1,RawScalar>::value();
40
41 typedef DualNumber<RawScalar, XVector> XFirstDerivType;
42 typedef DualNumber<RawScalar, YVector> YFirstDerivType;
43 // typedef DualNumber<XFirstDerivType, XVector::rebind<XFirstDerivType>::other > XSecondDerivType;
44 // typedef DualNumber<YFirstDerivType, YVector::rebind<YFirstDerivType>::other > YSecondDerivType;
45
46 typedef XFirstDerivType XADType;
47 typedef YFirstDerivType YADType;
48
49 typedef VectorOf<NDIM, 0, XADType, 1, YADType>::type Vector;
50
51 #ifdef METAPHYSICL_HAVE_MASA
52 // initialize the problem in MASA
53 err += masa_init("euler-maple","euler_2d");
54
55 // call the sanity check routine
56 // (tests that all variables have been initialized)
57 err += masa_sanity_check();
58 //err += masa_printid<double>();
59 #endif // METAPHYSICL_HAVE_MASA
60
61 Vector xy;
62 xy.insert<0>() = XADType(1., xvec);
63 xy.insert<1>() = YADType(1., yvec);
64
65 // the input argument xyz is another NumberVector
66 // a vector just like Q_rho_u, a spatial location rather
67 // than a vector-valued forcing function.
68 double h = 1.0/N;
69 for (int i=0; i != N+1; ++i)
70 {
71 xy.get<0>() = XADType(i*h, xvec);
72
73 for (int j=0; j != N+1; ++j)
74 {
75 xy.get<1>() = YADType(j*h, yvec);
76
77 // AD source terms
78 s2u = evaluate_q(xy,1);
79 s2v = evaluate_q(xy,2);
80 s2p = evaluate_q(xy,3);
81 s2e = evaluate_q(xy,4);
82
83 #ifdef METAPHYSICL_HAVE_MASA
84 // evaluate masa source terms
85 su = masa_eval_source_rho_u<double>(i*h,j*h);
86 sv = masa_eval_source_rho_v<double>(i*h,j*h);
87 sp = masa_eval_source_rho <double>(i*h,j*h);
88 se = masa_eval_source_rho_e<double>(i*h,j*h);
89
90 unorm = fabs(su-s2u);
91 vnorm = fabs(sv-s2v);
92 pnorm = fabs(sp-s2p);
93 enorm = fabs(se-s2e);
94
95 double urnorm = fabs(su-s2u)/std::max(su,s2u);
96 double vrnorm = fabs(sv-s2v)/std::max(sv,s2v);
97 double prnorm = fabs(sp-s2p)/std::max(sp,s2p);
98 double ernorm = fabs(se-s2e)/std::max(se,s2e);
99
100 unorm_max = std::max(unorm, unorm_max);
101 vnorm_max = std::max(vnorm, vnorm_max);
102 pnorm_max = std::max(pnorm, pnorm_max);
103 enorm_max = std::max(enorm, enorm_max);
104
105 urnorm_max = std::max(urnorm, urnorm_max);
106 vrnorm_max = std::max(vrnorm, vrnorm_max);
107 prnorm_max = std::max(prnorm, prnorm_max);
108 ernorm_max = std::max(ernorm, ernorm_max);
109 #else
110 // Avoid "set but not used" variable warnings;
111 (void) s2u;
112 (void) s2v;
113 (void) s2p;
114 (void) s2e;
115 #endif // METAPHYSICL_HAVE_MASA
116
117 }
118 }
119
120 #ifdef METAPHYSICL_HAVE_MASA
121 std::cout << "max error in u : " << unorm_max << std::endl;
122 std::cout << "max error in v : " << vnorm_max << std::endl;
123 std::cout << "max error in density: " << pnorm_max << std::endl;
124 std::cout << "max error in energy : " << enorm_max << std::endl;
125
126 std::cout << "max relative error in u : " << urnorm_max << std::endl;
127 std::cout << "max relative error in v : " << vrnorm_max << std::endl;
128 std::cout << "max relative error in density: " << prnorm_max << std::endl;
129 std::cout << "max relative error in energy : " << ernorm_max << std::endl;
130 #endif // METAPHYSICL_HAVE_MASA
131
132 // steady as she goes...
133 return 0;
134
135 }
136
137 // Note: ADScalar needs to be a FirstDerivType or better since we have
138 // first derivatives here. Adding diffusion will require a
139 // SecondDerivType or better
140
141 template <typename Vector>
evaluate_q(const Vector & xyz,const int ret)142 double evaluate_q (const Vector& xyz, const int ret)
143 {
144 typedef typename Vector::value_type ADScalar;
145
146 typedef typename RawType<ADScalar>::value_type Scalar;
147
148 typedef typename Vector::template rebind<Scalar>::other RawVector;
149
150 typedef typename Vector::template rebind<ADScalar>::other FullVector;
151
152 const Scalar PI = std::acos(Scalar(-1));
153
154 const Scalar u_0 = 200.23;
155 const Scalar u_x = 1.1;
156 const Scalar u_y = 1.08;
157 const Scalar v_0 = 1.2;
158 const Scalar v_x = 1.6;
159 const Scalar v_y = .47;
160 const Scalar rho_0 = 100.02;
161 const Scalar rho_x = 2.22;
162 const Scalar rho_y = 0.8;
163 const Scalar p_0 = 150.2;
164 const Scalar p_x = .91;
165 const Scalar p_y = .623;
166 const Scalar a_px = .165;
167 const Scalar a_py = .612;
168 const Scalar a_rhox = 1.0;
169 const Scalar a_rhoy = 1.0;
170 const Scalar a_ux = .1987;
171 const Scalar a_uy = 1.189;
172 const Scalar a_vx = 1.91;
173 const Scalar a_vy = 1.0;
174 const Scalar Gamma = 1.01;
175 const Scalar L = 3.02;
176
177 const typename Vector::template entry_type<0>::type& x =
178 xyz.template get<0>();
179 const typename Vector::template entry_type<1>::type& y =
180 xyz.template get<1>();
181
182 // Treat velocity as a vector
183 FullVector U;
184
185 // Arbitrary manufactured solution
186 U.template insert<0>() = u_0 + u_x * std::sin(a_ux * PI * x / L) + u_y * std::cos(a_uy * PI * y / L);
187 U.template insert<1>() = v_0 + v_x * std::cos(a_vx * PI * x / L) + v_y * std::sin(a_vy * PI * y / L);
188 ADScalar RHO = rho_0 + rho_x * std::sin(a_rhox * PI * x / L) + rho_y * std::cos(a_rhoy * PI * y / L);
189 ADScalar P = p_0 + p_x * std::cos(a_px * PI * x / L) + p_y * std::sin(a_py * PI * y / L);
190
191 P/RHO;
192
193
194
195
196
197
198 // Perfect gas energies
199 ADScalar E = 1./(Gamma-1.)*P/RHO;
200 ADScalar ET = E + .5 * U.dot(U);
201
202 // Euler equation residuals
203 Scalar Q_rho = raw_value(divergence(RHO*U));
204 RawVector Q_rho_u = raw_value(divergence(RHO*U.outerproduct(U)) + P.derivatives());
205
206 // energy equation
207 Scalar Q_rho_e = raw_value(divergence((RHO*ET+P)*U));
208
209 switch(ret)
210 {
211
212 // u
213 case 1:
214 return Q_rho_u.template get<0>();
215 break;
216
217 // v
218 case 2:
219 return Q_rho_u.template get<1>();
220 break;
221
222 // rho
223 case 3:
224 return Q_rho;
225 break;
226
227 // energy
228 case 4:
229 return Q_rho_e;
230 break;
231
232 default:
233 throw std::domain_error("Bad evaluate_q input request");
234 }
235 return 0.;
236 }
237