Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/Physics/Problems.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.,
7 : : 2022-2024 J. Bakosi
8 : : All rights reserved. See the LICENSE file for details.
9 : : \brief Problem-specific functions. Initial conditions, source terms.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "Problems.hpp"
14 : : #include "EOS.hpp"
15 : : #include "InciterConfig.hpp"
16 : : #include "Box.hpp"
17 : :
18 : : namespace inciter {
19 : :
20 : : extern ctr::Config g_cfg;
21 : :
22 : : } // ::inciter
23 : :
24 : : namespace problems {
25 : :
26 : : using inciter::g_cfg;
27 : :
28 : : namespace userdef {
29 : :
30 : : static std::vector< tk::real >
31 : 1320556 : ic( tk::real, tk::real, tk::real, tk::real )
32 : : // *****************************************************************************
33 : : //! Set homogeneous initial conditions for a generic user-defined problem
34 : : //! \return Values of conserved variables
35 : : // *****************************************************************************
36 : : {
37 : : // pressure-based solvers
38 : :
39 : 1320556 : const auto& solver = g_cfg.get< tag::solver >();
40 [ + + ]: 1320556 : if (solver == "chocg") {
41 : 692956 : const auto& ncomp = g_cfg.get< tag::problem_ncomp >();
42 [ + - ]: 692956 : std::vector< tk::real > u( ncomp, 0.0 );
43 [ + - ]: 692956 : auto ic_velocity = g_cfg.get< tag::ic_velocity >();
44 : 692956 : auto large = std::numeric_limits< double >::max();
45 [ + - ]: 692956 : if (std::abs(ic_velocity[0] - large) > 1.0e-12) u[0] = ic_velocity[0];
46 [ + - ]: 692956 : if (std::abs(ic_velocity[1] - large) > 1.0e-12) u[1] = ic_velocity[1];
47 [ + - ]: 692956 : if (std::abs(ic_velocity[2] - large) > 1.0e-12) u[2] = ic_velocity[2];
48 : 692956 : return u;
49 : 692956 : }
50 [ + + ]: 627600 : else if (solver == "lohcg") {
51 : 511426 : const auto& ncomp = g_cfg.get< tag::problem_ncomp >();
52 [ + - ]: 511426 : std::vector< tk::real > u( ncomp, 0.0 );
53 [ + - ]: 511426 : auto ic_velocity = g_cfg.get< tag::ic_velocity >();
54 : 511426 : auto large = std::numeric_limits< double >::max();
55 [ + - ]: 511426 : if (std::abs(ic_velocity[0] - large) > 1.0e-12) u[1] = ic_velocity[0];
56 [ + - ]: 511426 : if (std::abs(ic_velocity[1] - large) > 1.0e-12) u[2] = ic_velocity[1];
57 [ + - ]: 511426 : if (std::abs(ic_velocity[2] - large) > 1.0e-12) u[3] = ic_velocity[2];
58 : 511426 : return u;
59 : 511426 : }
60 : :
61 : : // density-based solvers
62 : :
63 : 116174 : auto ic_density = g_cfg.get< tag::ic_density >();
64 : 116174 : const auto& ic_velocity = g_cfg.get< tag::ic_velocity >();
65 [ - + ][ - - ]: 116174 : ErrChk( ic_velocity.size() == 3, "ic_velocity must have 3 components" );
[ - - ][ - - ]
66 : :
67 [ + - ]: 116174 : std::vector< tk::real > u( 5, 0.0 );
68 : :
69 : 116174 : u[0] = ic_density;
70 : 116174 : u[1] = u[0] * ic_velocity[0];
71 : 116174 : u[2] = u[0] * ic_velocity[1];
72 : 116174 : u[3] = u[0] * ic_velocity[2];
73 : :
74 : 116174 : auto ic_pressure = g_cfg.get< tag::ic_pressure >();
75 : 116174 : auto ic_energy = g_cfg.get< tag::ic_energy >();
76 : 116174 : auto ic_temperature = g_cfg.get< tag::ic_temperature >();
77 : :
78 : 116174 : auto largereal = std::numeric_limits< double >::max();
79 : :
80 [ + - ]: 116174 : if (std::abs(ic_pressure - largereal) > 1.0e-12) {
81 : :
82 : 116174 : u[4] = eos::totalenergy( u[0], u[1]/u[0], u[2]/u[0], u[3]/u[0],
83 : : ic_pressure );
84 : :
85 [ - - ]: 0 : } else if (std::abs(ic_energy - largereal) > 1.0e-12) {
86 : :
87 : 0 : u[4] = u[0] * ic_energy;
88 : :
89 [ - - ]: 0 : } else if (std::abs(ic_temperature - largereal) > 1.0e-12) {
90 : :
91 : 0 : auto cv = g_cfg.get< tag::mat_spec_heat_const_vol >();
92 [ - - ]: 0 : if (std::abs(cv - largereal) > 1.0e-12) {
93 : 0 : u[4] = u[0] * ic_temperature * cv;
94 : : }
95 : :
96 : : } else {
97 : :
98 [ - - ][ - - ]: 0 : Throw( "IC background energy cannot be computed. Must specify "
[ - - ]
99 : : "one of background pressure, energy, or velocity." );
100 : :
101 : : }
102 : :
103 : 116174 : return u;
104 : 116174 : }
105 : :
106 : : static tk::real
107 : 186 : pic( tk::real, tk::real, tk::real )
108 : : // *****************************************************************************
109 : : //! Set homogeneous initial conditions for a generic user-defined problem
110 : : //! \return Value of pressure
111 : : // *****************************************************************************
112 : : {
113 : 186 : return 0.0;
114 : : }
115 : :
116 : : } // userdef::
117 : :
118 : : namespace nonlinear_energy_growth {
119 : :
120 : : static std::vector< tk::real >
121 : 73908 : ic( tk::real x, tk::real y, tk::real z, tk::real t )
122 : : // *****************************************************************************
123 : : //! Set initial conditions prescribing nonlinear energy growth
124 : : //! \param[in] x X coordinate where to evaluate the solution
125 : : //! \param[in] y Y coordinate where to evaluate the solution
126 : : //! \param[in] z Z coordinate where to evaluate the solution
127 : : //! \param[in] t Time where to evaluate the solution
128 : : //! \return Values of conserved variables
129 : : // *****************************************************************************
130 : : {
131 : : using std::cos;
132 : :
133 : : // manufactured solution parameters
134 : 73908 : auto ce = g_cfg.get< tag::problem_ce >();
135 : 73908 : auto r0 = g_cfg.get< tag::problem_r0 >();
136 : 73908 : auto a = g_cfg.get< tag::problem_alpha >();
137 : 73908 : auto k = g_cfg.get< tag::problem_kappa >();
138 : 73908 : const auto& b = g_cfg.get< tag::problem_beta >();
139 : :
140 : 73908 : auto ec = [ ce, t ]( tk::real kappa, tk::real h, tk::real p ) {
141 : 73908 : return std::pow( -3.0*(ce + kappa*h*h*t), p );
142 : 73908 : };
143 : :
144 : 295632 : auto hx = [ x, y, z, b ]() {
145 : 73908 : return cos(b[0]*M_PI*x) * cos(b[1]*M_PI*y) * cos(b[2]*M_PI*z);
146 [ + - ]: 73908 : };
147 : :
148 : : // density
149 : 73908 : auto r = r0 + std::exp(-a*t) * (1.0 - x*x - y*y - z*z);
150 : : // energy
151 : 73908 : auto re = r * ec(k,hx(),-1.0/3.0);
152 : :
153 [ + - ]: 147816 : return { r, 0.0, 0.0, 0.0, re };
154 : 73908 : }
155 : :
156 : : static std::vector< tk::real >
157 : 178308 : src( tk::real x, tk::real y, tk::real z, tk::real t )
158 : : // *****************************************************************************
159 : : //! Compute and return source term for nonlinear energy growth
160 : : //! \param[in] x X coordinate where to evaluate the source
161 : : //! \param[in] y Y coordinate where to evaluate the source
162 : : //! \param[in] z Z coordinate where to evaluate the source
163 : : //! \param[in] t Time where to evaluate the source
164 : : //! \return Source for flow variables + transported scalars
165 : : // *****************************************************************************
166 : : {
167 : : using std::sin; using std::cos; using std::pow;
168 : :
169 : : // manufactured solution parameters
170 : 178308 : auto a = g_cfg.get< tag::problem_alpha >();
171 : 178308 : const auto& b = g_cfg.get< tag::problem_beta >();
172 : 178308 : auto ce = g_cfg.get< tag::problem_ce >();
173 : 178308 : auto kappa = g_cfg.get< tag::problem_kappa >();
174 : 178308 : auto r0 = g_cfg.get< tag::problem_r0 >();
175 : : // ratio of specific heats
176 : 178308 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
177 : : // spatial component of density field
178 : 178308 : auto gx = 1.0 - x*x - y*y - z*z;
179 : : // derivative of spatial component of density field
180 : 178308 : std::array< tk::real, 3 > dg{ -2.0*x, -2.0*y, -2.0*z };
181 : : // spatial component of energy field
182 : 178308 : auto h = cos(b[0]*M_PI*x) * cos(b[1]*M_PI*y) * cos(b[2]*M_PI*z);
183 : : // derivative of spatial component of energy field
184 : : std::array< tk::real, 3 >
185 : 178308 : dh{ -b[0]*M_PI*sin(b[0]*M_PI*x)*cos(b[1]*M_PI*y)*cos(b[2]*M_PI*z),
186 : 178308 : -b[1]*M_PI*cos(b[0]*M_PI*x)*sin(b[1]*M_PI*y)*cos(b[2]*M_PI*z),
187 : 356616 : -b[2]*M_PI*cos(b[0]*M_PI*x)*cos(b[1]*M_PI*y)*sin(b[2]*M_PI*z) };
188 : : // temporal function f and its derivative
189 : 178308 : auto ft = std::exp(-a*t);
190 : 178308 : auto dfdt = -a*ft;
191 : : // density and its derivatives
192 : 178308 : auto rho = r0 + ft*gx;
193 : 178308 : std::array< tk::real, 3 > drdx{ ft*dg[0], ft*dg[1], ft*dg[2] };
194 : 178308 : auto drdt = gx*dfdt;
195 : : // internal energy and its derivatives
196 : 178308 : auto ie = pow( -3.0*(ce + kappa*h*h*t), -1.0/3.0 );
197 : 178308 : std::array< tk::real, 3 > dedx{ 2.0 * pow(ie,4.0) * kappa * h * dh[0] * t,
198 : 356616 : 2.0 * pow(ie,4.0) * kappa * h * dh[1] * t,
199 : 178308 : 2.0 * pow(ie,4.0) * kappa * h * dh[2] * t };
200 : 178308 : const auto dedt = kappa * h * h * pow(ie,4.0);
201 : :
202 [ + - ]: 178308 : std::vector< tk::real > s( 5, 0.0 );
203 : : // density source
204 : 178308 : s[0] = drdt;
205 : : // momentum source
206 : 178308 : s[1] = (g-1.0)*(rho*dedx[0] + ie*drdx[0]);
207 : 178308 : s[2] = (g-1.0)*(rho*dedx[1] + ie*drdx[1]);
208 : 178308 : s[3] = (g-1.0)*(rho*dedx[2] + ie*drdx[2]);
209 : : // energy source
210 : 178308 : s[4] = rho*dedt + ie*drdt;
211 : :
212 : 356616 : return s;
213 : : }
214 : :
215 : : } // nonlinear_energy_growth::
216 : :
217 : : namespace rayleigh_taylor {
218 : :
219 : : static std::vector< tk::real >
220 : 1058942 : ic( tk::real x, tk::real y, tk::real z, tk::real t )
221 : : // *****************************************************************************
222 : : //! Set initial conditions prescribing a Rayleigh-Taylor flow
223 : : //! \param[in] x X coordinate where to evaluate the solution
224 : : //! \param[in] y Y coordinate where to evaluate the solution
225 : : //! \param[in] z Z coordinate where to evaluate the solution
226 : : //! \param[in] t Time where to evaluate the solution
227 : : //! \return Values of conserved variables
228 : : // *****************************************************************************
229 : : {
230 : : using std::sin; using std::cos;
231 : :
232 : : // manufactured solution parameters
233 : 1058942 : auto a = g_cfg.get< tag::problem_alpha >();
234 : 1058942 : const auto& b = g_cfg.get< tag::problem_beta >();
235 : 1058942 : auto p0 = g_cfg.get< tag::problem_p0 >();
236 : 1058942 : auto r0 = g_cfg.get< tag::problem_r0 >();
237 : 1058942 : auto k = g_cfg.get< tag::problem_kappa >();
238 : :
239 : : // spatial component of density and pressure fields
240 : 1058942 : tk::real gx = b[0]*x*x + b[1]*y*y + b[2]*z*z;
241 : : // density
242 : 1058942 : tk::real r = r0 - gx;
243 : : // velocity
244 : 1058942 : tk::real ft = cos(k*M_PI*t);
245 : 1058942 : tk::real u = ft * z * sin(M_PI*x);
246 : 1058942 : tk::real v = ft * z * cos(M_PI*y);
247 : 1058942 : tk::real w = ft * ( -0.5*M_PI*z*z*(cos(M_PI*x) - sin(M_PI*y)) );
248 : : // total specific energy
249 : 1058942 : tk::real rE = eos::totalenergy( r, u, v, w, p0 + a*gx );
250 : :
251 [ + - ]: 1058942 : return { r, r*u, r*v, r*w, rE };
252 : : }
253 : :
254 : : static std::vector< tk::real >
255 : 768460 : src( tk::real x, tk::real y, tk::real z, tk::real t )
256 : : // *****************************************************************************
257 : : //! Compute and return source term for a Rayleigh-Taylor flow
258 : : //! \param[in] x X coordinate where to evaluate the source
259 : : //! \param[in] y Y coordinate where to evaluate the source
260 : : //! \param[in] z Z coordinate where to evaluate the source
261 : : //! \param[in] t Time where to evaluate the source
262 : : //! \return Source for flow variables + transported scalars
263 : : // *****************************************************************************
264 : : {
265 : : using std::sin; using std::cos;
266 : :
267 : : // manufactured solution parameters
268 : 768460 : auto a = g_cfg.get< tag::problem_alpha >();
269 : 768460 : const auto& b = g_cfg.get< tag::problem_beta >();
270 : 768460 : auto k = g_cfg.get< tag::problem_kappa >();
271 : 768460 : auto p0 = g_cfg.get< tag::problem_p0 >();
272 : 768460 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
273 : :
274 : : // evaluate solution at x,y,z,t
275 [ + - ]: 768460 : auto U = ic( x, y, z, t );
276 : :
277 : : // density, velocity, energy, pressure
278 : 768460 : auto rho = U[0];
279 : 768460 : auto u = U[1]/U[0];
280 : 768460 : auto v = U[2]/U[0];
281 : 768460 : auto w = U[3]/U[0];
282 : 768460 : auto E = U[4]/U[0];
283 : 768460 : auto p = p0 + a*(b[0]*x*x + b[1]*y*y + b[2]*z*z);
284 : :
285 : : // spatial gradients
286 : 768460 : std::array< tk::real, 3 > drdx{{ -2.0*b[0]*x, -2.0*b[1]*y, -2.0*b[2]*z }};
287 : 768460 : std::array< tk::real, 3 > dpdx{{ 2.0*a*b[0]*x, 2.0*a*b[1]*y, 2.0*a*b[2]*z }};
288 : 768460 : tk::real ft = cos(k*M_PI*t);
289 : 768460 : std::array< tk::real, 3 > dudx{{ ft*M_PI*z*cos(M_PI*x),
290 : : 0.0,
291 : 768460 : ft*sin(M_PI*x) }};
292 : 768460 : std::array< tk::real, 3 > dvdx{{ 0.0,
293 : 768460 : -ft*M_PI*z*sin(M_PI*y),
294 : 768460 : ft*cos(M_PI*y) }};
295 : 768460 : std::array< tk::real, 3 > dwdx{{ ft*M_PI*0.5*M_PI*z*z*sin(M_PI*x),
296 : 768460 : ft*M_PI*0.5*M_PI*z*z*cos(M_PI*y),
297 : 768460 : -ft*M_PI*z*(cos(M_PI*x) - sin(M_PI*y)) }};
298 : : std::array< tk::real, 3 > dedx{{
299 : 768460 : dpdx[0]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[0]
300 : 768460 : + u*dudx[0] + v*dvdx[0] + w*dwdx[0],
301 : 768460 : dpdx[1]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[1]
302 : 768460 : + u*dudx[1] + v*dvdx[1] + w*dwdx[1],
303 : 768460 : dpdx[2]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[2]
304 : 2305380 : + u*dudx[2] + v*dvdx[2] + w*dwdx[2] }};
305 : :
306 : : // time derivatives
307 : 768460 : auto dudt = -k*M_PI*sin(k*M_PI*t)*z*sin(M_PI*x);
308 : 768460 : auto dvdt = -k*M_PI*sin(k*M_PI*t)*z*cos(M_PI*y);
309 : 768460 : auto dwdt = k*M_PI*sin(k*M_PI*t)/2*M_PI*z*z*(cos(M_PI*x) - sin(M_PI*y));
310 : 768460 : auto dedt = u*dudt + v*dvdt + w*dwdt;
311 : :
312 [ + - ]: 768460 : std::vector< tk::real > s( 5, 0.0 );
313 : : // density source
314 : 768460 : s[0] = u*drdx[0] + v*drdx[1] + w*drdx[2];
315 : : // momentum source
316 : 768460 : s[1] = rho*dudt+u*s[0]+dpdx[0] + U[1]*dudx[0]+U[2]*dudx[1]+U[3]*dudx[2];
317 : 768460 : s[2] = rho*dvdt+v*s[0]+dpdx[1] + U[1]*dvdx[0]+U[2]*dvdx[1]+U[3]*dvdx[2];
318 : 768460 : s[3] = rho*dwdt+w*s[0]+dpdx[2] + U[1]*dwdx[0]+U[2]*dwdx[1]+U[3]*dwdx[2];
319 : : // energy source
320 : 768460 : s[4] = rho*dedt + E*s[0] + U[1]*dedx[0]+U[2]*dedx[1]+U[3]*dedx[2]
321 : 768460 : + u*dpdx[0]+v*dpdx[1]+w*dpdx[2];
322 : :
323 : 1536920 : return s;
324 : 768460 : }
325 : :
326 : : } // rayleigh_taylor::
327 : :
328 : : namespace sedov {
329 : : static std::vector< tk::real >
330 : 65767 : ic( tk::real x, tk::real y, tk::real z, tk::real )
331 : : // *****************************************************************************
332 : : //! Set initial conditions prescribing the Sedov blast wave
333 : : //! \param[in] x X coordinate where to evaluate the solution
334 : : //! \param[in] y Y coordinate where to evaluate the solution
335 : : //! \param[in] z Z coordinate where to evaluate the solution
336 : : //! \return Values of conserved variables
337 : : // *****************************************************************************
338 : : {
339 : : using std::abs;
340 : :
341 : : // pressure
342 : 65767 : auto eps = std::numeric_limits< tk::real >::epsilon();
343 : : tk::real p;
344 [ + + ][ + + ]: 65767 : if (abs(x) < eps && abs(y) < eps && abs(z) < eps) {
[ + + ][ + + ]
345 : 5 : p = g_cfg.get< tag::problem_p0 >();
346 : : } else {
347 : 65762 : p = 0.67e-4;
348 : : }
349 : :
350 : : // density
351 : 65767 : tk::real r = 1.0;
352 : : // velocity
353 : 65767 : tk::real u = 0.0;
354 : 65767 : tk::real v = 0.0;
355 : 65767 : tk::real w = 0.0;
356 : : // total specific energy
357 : 65767 : tk::real rE = eos::totalenergy( r, u, v, w, p );
358 : :
359 [ + - ]: 65767 : return { r, r*u, r*v, r*w, rE };
360 : :
361 : : }
362 : : } // sedov::
363 : :
364 : : namespace sod {
365 : : static std::vector< tk::real >
366 : 13179 : ic( tk::real x, tk::real, tk::real, tk::real )
367 : : // *****************************************************************************
368 : : //! Set initial conditions prescribing the Sod shocktube
369 : : //! \param[in] x X coordinate where to evaluate the solution
370 : : //! \return Values of conserved variables
371 : : // *****************************************************************************
372 : : {
373 : : tk::real r, p, u, v, w, rE;
374 : :
375 [ + + ]: 13179 : if (x<0.5) {
376 : : // density
377 : 6662 : r = 1.0;
378 : : // pressure
379 : 6662 : p = 1.0;
380 : : }
381 : : else {
382 : : // density
383 : 6517 : r = 0.125;
384 : : // pressure
385 : 6517 : p = 0.1;
386 : : }
387 : :
388 : : // velocity
389 : 13179 : u = 0.0;
390 : 13179 : v = 0.0;
391 : 13179 : w = 0.0;
392 : :
393 : : // total specific energy
394 : 13179 : rE = eos::totalenergy( r, u, v, w, p );
395 : :
396 [ + - ]: 13179 : return { r, r*u, r*v, r*w, rE };
397 : : }
398 : : } // sod::
399 : :
400 : : namespace taylor_green {
401 : :
402 : : static std::vector< tk::real >
403 : 280356 : ic( tk::real x, tk::real y, tk::real, tk::real )
404 : : // *****************************************************************************
405 : : //! Set initial conditions prescribing the Taylor-Green vortex
406 : : //! \param[in] x X coordinate where to evaluate the solution
407 : : //! \param[in] y Y coordinate where to evaluate the solution
408 : : //! \return Values of conserved variables
409 : : // *****************************************************************************
410 : : {
411 : : // density
412 : 280356 : tk::real r = 1.0;
413 : : // pressure
414 : 280356 : tk::real p = 10.0 + r/4.0*(cos(2.0*M_PI*x) + cos(2.0*M_PI*y));
415 : : // velocity
416 : 280356 : tk::real u = sin(M_PI*x) * cos(M_PI*y);
417 : 280356 : tk::real v = -cos(M_PI*x) * sin(M_PI*y);
418 : 280356 : tk::real w = 0.0;
419 : : // total specific energy
420 : 280356 : auto rE = eos::totalenergy( r, u, v, w, p );
421 : :
422 [ + - ]: 280356 : return { r, r*u, r*v, r*w, rE };
423 : : }
424 : :
425 : : static std::vector< tk::real >
426 : 932688 : src( tk::real x, tk::real y, tk::real, tk::real )
427 : : // *****************************************************************************
428 : : //! Compute and return source term for a the Taylor-Green vortex
429 : : //! \param[in] x X coordinate where to evaluate the source
430 : : //! \param[in] y Y coordinate where to evaluate the source
431 : : //! \return Source for flow variables + transported scalars
432 : : // *****************************************************************************
433 : : {
434 : : using std::cos;
435 : :
436 [ + - ]: 932688 : std::vector< tk::real > s( 5, 0.0 );
437 : 1865376 : s[4] = 3.0*M_PI/8.0*( cos(3.0*M_PI*x)*cos(M_PI*y)
438 : 932688 : - cos(3.0*M_PI*y)*cos(M_PI*x) );
439 : :
440 : 932688 : return s;
441 : : }
442 : :
443 : : } // taylor_green::
444 : :
445 : : namespace vortical_flow {
446 : :
447 : : static std::vector< tk::real >
448 : 2350800 : ic( tk::real x, tk::real y, tk::real z, tk::real )
449 : : // *****************************************************************************
450 : : //! Set initial conditions prescribing vortical flow
451 : : //! \param[in] x X coordinate where to evaluate the solution
452 : : //! \param[in] y Y coordinate where to evaluate the solution
453 : : //! \param[in] z Z coordinate where to evaluate the solution
454 : : //! \return Values of conserved variables
455 : : // *****************************************************************************
456 : : {
457 : : // manufactured solution parameters
458 : 2350800 : tk::real a = g_cfg.get< tag::problem_alpha >();
459 : 2350800 : tk::real k = g_cfg.get< tag::problem_kappa >();
460 : 2350800 : tk::real p0 = g_cfg.get< tag::problem_p0 >();
461 : : // ratio of specific heats
462 : 2350800 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
463 : : // velocity
464 : 2350800 : tk::real ru = a*x - k*y;
465 : 2350800 : tk::real rv = k*x + a*y;
466 : 2350800 : tk::real rw = -2.0*a*z;
467 : : // total specific energy
468 : 2350800 : tk::real rE = (ru*ru + rv*rv + rw*rw)/2.0 + (p0 - 2.0*a*a*z*z) / (g - 1.0);
469 : :
470 [ + - ]: 2350800 : return { 1.0, ru, rv, rw, rE };
471 : : }
472 : :
473 : : static std::vector< tk::real >
474 : 1474909 : src( tk::real x, tk::real y, tk::real z, tk::real )
475 : : // *****************************************************************************
476 : : //! Compute and return source term for vortical flow
477 : : //! \param[in] x X coordinate where to evaluate the source
478 : : //! \param[in] y Y coordinate where to evaluate the source
479 : : //! \param[in] z Z coordinate where to evaluate the source
480 : : //! \return Source for flow variables + transported scalars
481 : : // *****************************************************************************
482 : : {
483 : : // manufactured solution parameters
484 : 1474909 : auto a = g_cfg.get< tag::problem_alpha >();
485 : 1474909 : auto k = g_cfg.get< tag::problem_kappa >();
486 : : // ratio of specific heats
487 : 1474909 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
488 : : // evaluate solution at x,y,z
489 [ + - ]: 1474909 : auto u = ic( x, y, z, 0.0 );
490 : :
491 [ + - ]: 1474909 : std::vector< tk::real > s( 5, 0.0 );
492 : : // momentum source
493 : 1474909 : s[1] = a*u[1]/u[0] - k*u[2]/u[0];
494 : 1474909 : s[2] = k*u[1]/u[0] + a*u[2]/u[0];
495 : : // energy source
496 : 1474909 : s[4] = (s[1]*u[1] + s[2]*u[2])/u[0] + 8.0*a*a*a*z*z/(g-1.0);
497 : :
498 : 2949818 : return s;
499 : 1474909 : }
500 : :
501 : : } // vortical_flow::
502 : :
503 : : namespace slot_cyl {
504 : :
505 : : static std::vector< tk::real >
506 : 4300699 : ic( tk::real x, tk::real y, tk::real, tk::real t )
507 : : // *****************************************************************************
508 : : //! Set initial conditions prescribing slotted cylinder, cone, Gauss hump
509 : : //! \param[in] x X coordinate where to evaluate the solution
510 : : //! \param[in] y Y coordinate where to evaluate the solution
511 : : //! \param[in] t Time where to evaluate the solution
512 : : //! \return Values of conserved variables
513 : : // *****************************************************************************
514 : : {
515 : : using std::sin; using std::cos; using std::sqrt;
516 : :
517 : : // manufactured solution parameters
518 : 4300699 : tk::real p0 = 1.0;
519 : :
520 [ + - ]: 4300699 : std::vector< tk::real > u( 6, 0.0 );
521 : :
522 : : // prescribed velocity: rotate in x-y plane
523 : 4300699 : u[0] = 1.0;
524 : 4300699 : u[1] = u[0] * (0.5 - y);
525 : 4300699 : u[2] = u[0] * (x - 0.5);
526 : 4300699 : u[3] = 0.0;
527 : 4300699 : u[4] = eos::totalenergy( u[0], u[1]/u[0], u[2]/u[0], u[3]/u[0], p0 );
528 : :
529 : 4300699 : const tk::real R0 = 0.15;
530 : :
531 : : // center of the cone
532 : 4300699 : tk::real x0 = 0.5;
533 : 4300699 : tk::real y0 = 0.25;
534 : 4300699 : tk::real r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
535 : 4300699 : tk::real kx = 0.5 + r*sin( t );
536 : 4300699 : tk::real ky = 0.5 - r*cos( t );
537 : :
538 : : // center of the hump
539 : 4300699 : x0 = 0.25;
540 : 4300699 : y0 = 0.5;
541 : 4300699 : r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
542 : 4300699 : tk::real hx = 0.5 + r*sin( t-M_PI/2.0 ),
543 : 4300699 : hy = 0.5 - r*cos( t-M_PI/2.0 );
544 : :
545 : : // center of the slotted cylinder
546 : 4300699 : x0 = 0.5;
547 : 4300699 : y0 = 0.75;
548 : 4300699 : r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
549 : 4300699 : tk::real cx = 0.5 + r*sin( t+M_PI ),
550 : 4300699 : cy = 0.5 - r*cos( t+M_PI );
551 : :
552 : : // end points of the cylinder slot
553 : 4300699 : tk::real i1x = 0.525, i1y = cy - r*cos( std::asin(0.025/r) ),
554 : 4300699 : i2x = 0.525, i2y = 0.8,
555 : 4300699 : i3x = 0.475, i3y = 0.8;
556 : :
557 : : // rotate end points of cylinder slot
558 : 4300699 : tk::real ri1x = 0.5 + cos(t)*(i1x-0.5) - sin(t)*(i1y-0.5),
559 : 4300699 : ri1y = 0.5 + sin(t)*(i1x-0.5) + cos(t)*(i1y-0.5),
560 : 4300699 : ri2x = 0.5 + cos(t)*(i2x-0.5) - sin(t)*(i2y-0.5),
561 : 4300699 : ri2y = 0.5 + sin(t)*(i2x-0.5) + cos(t)*(i2y-0.5),
562 : 4300699 : ri3x = 0.5 + cos(t)*(i3x-0.5) - sin(t)*(i3y-0.5),
563 : 4300699 : ri3y = 0.5 + sin(t)*(i3x-0.5) + cos(t)*(i3y-0.5);
564 : :
565 : : // direction of slot sides
566 : 4300699 : tk::real v1x = ri2x-ri1x, v1y = ri2y-ri1y,
567 : 4300699 : v2x = ri3x-ri2x, v2y = ri3y-ri2y;
568 : :
569 : : // lengths of direction of slot sides vectors
570 : 4300699 : tk::real v1 = sqrt(v1x*v1x + v1y*v1y),
571 : 4300699 : v2 = sqrt(v2x*v2x + v2y*v2y);
572 : :
573 : : // cone
574 : 4300699 : r = sqrt((x-kx)*(x-kx) + (y-ky)*(y-ky)) / R0;
575 [ + + ]: 4300699 : if (r<1.0) u[5] = 0.6*(1.0-r);
576 : :
577 : : // hump
578 : 4300699 : r = sqrt((x-hx)*(x-hx) + (y-hy)*(y-hy)) / R0;
579 [ + + ]: 4300699 : if (r<1.0) u[5] = 0.2*(1.0+cos(M_PI*std::min(r,1.0)));
580 : :
581 : : // cylinder
582 : 4300699 : r = sqrt((x-cx)*(x-cx) + (y-cy)*(y-cy)) / R0;
583 : 4300699 : const std::array< tk::real, 2 > r1{{ v1x, v1y }},
584 : 4300699 : r2{{ x-ri1x, y-ri1y }};
585 : 4300699 : const auto d1 = (r1[0]*r2[1] - r2[0]*r1[1]) / v1;
586 : 4300699 : const std::array< tk::real, 2 > r3{{ v2x, v2y }},
587 : 4300699 : r4{{ x-ri2x, y-ri2y }};
588 : 4300699 : const auto d2 = (r3[0]*r4[1] - r4[0]*r3[1]) / v2;
589 [ + + ][ + + ]: 4300699 : if (r<1.0 && (d1>0.05 || d1<0.0 || d2<0.0)) u[5] = 0.6;
[ + + ][ + + ]
590 : :
591 : 8601398 : return u;
592 : : }
593 : :
594 : : static std::vector< tk::real >
595 : 2942970 : src( tk::real x, tk::real y, tk::real z, tk::real t )
596 : : // *****************************************************************************
597 : : //! Compute and return source term for slotted cylinder, cone, Gauss hump
598 : : //! \param[in] x X coordinate where to evaluate the source
599 : : //! \param[in] y Y coordinate where to evaluate the source
600 : : //! \param[in] z Z coordinate where to evaluate the source
601 : : //! \param[in] t Time where to evaluate the source
602 : : //! \return Source for flow variables + transported scalars
603 : : // *****************************************************************************
604 : : {
605 : : // evaluate solution at x,y,z,t
606 [ + - ]: 2942970 : auto u = ic( x, y, z, t );
607 : :
608 [ + - ]: 2942970 : std::vector< tk::real > s( 6, 0.0 );
609 : : // momentum source
610 : 2942970 : s[1] = -u[2];
611 : 2942970 : s[2] = u[1];
612 : :
613 : 5885940 : return s;
614 : 2942970 : }
615 : :
616 : : } // slot_cyl::
617 : :
618 : : namespace point_src {
619 : :
620 : : static std::vector< tk::real >
621 : 1858 : ic( tk::real x, tk::real y, tk::real z, tk::real t )
622 : : // *****************************************************************************
623 : : //! Set initial conditions for point source problem
624 : : //! \param[in] x X coordinate where to evaluate initial conditions
625 : : //! \param[in] y Y coordinate where to evaluate initial conditions
626 : : //! \param[in] z Z coordinate where to evaluate initial conditions
627 : : //! \param[in] t Time where to evaluate initial conditions
628 : : //! \return Values of conserved variables
629 : : // *****************************************************************************
630 : : {
631 : 1858 : auto u = userdef::ic( x, y, z, t );
632 [ + - ]: 1858 : u.push_back( 0.0 );
633 : 1858 : return u;
634 : 0 : }
635 : :
636 : : static void
637 : 600 : src( const std::array< std::vector< tk::real >, 3 >& coord,
638 : : tk::real t,
639 : : tk::Fields& U )
640 : : // *****************************************************************************
641 : : //! Apply point-source directly to numerical solution
642 : : //! \param[in] coord Mesh node coordinates
643 : : //! \param[in] t Physical time
644 : : //! \param[in,out] U Solution vector at recent time step
645 : : //! \note This is different from other source terms, because this directly
646 : : //! modifies the solution instead of applied as a source term mathematically.
647 : : //! Hence the function signature is also different.
648 : : // *****************************************************************************
649 : : {
650 [ - + ]: 600 : if (U.nprop() == 5) return;
651 : :
652 : 600 : const auto& source = g_cfg.get< tag::problem_src >();
653 : 600 : const auto& location = source.get< tag::location >();
654 : 600 : auto radius = source.get< tag::radius >();
655 : 600 : auto release_time = source.get< tag::release_time >();
656 : 600 : auto largereal = std::numeric_limits< double >::max();
657 : :
658 [ + - ]: 1200 : if (location.size() != 3 ||
659 [ + - ][ - + ]: 1200 : std::abs(radius - largereal) < 1.0e-12 ||
[ - + ]
660 : 600 : std::abs(release_time - largereal) < 1.0e-12)
661 : : {
662 : 0 : return;
663 : : }
664 : :
665 : 600 : auto sx = location[0];
666 : 600 : auto sy = location[1];
667 : 600 : auto sz = location[2];
668 : 600 : auto sr = radius;
669 : 600 : auto st = release_time;
670 : :
671 [ - + ]: 600 : if (t < st) return;
672 : :
673 : 600 : const auto& x = coord[0];
674 : 600 : const auto& y = coord[1];
675 : 600 : const auto& z = coord[2];
676 : :
677 [ + + ]: 279300 : for (std::size_t i=0; i<U.nunk(); ++i) {
678 : 278700 : auto rx = sx - x[i];
679 : 278700 : auto ry = sy - y[i];
680 : 278700 : auto rz = sz - z[i];
681 [ + + ]: 278700 : if (rx*rx + ry*ry + rz*rz < sr*sr) U(i,5) = 1.0;
682 : : }
683 : :
684 : 600 : return;
685 : : }
686 : :
687 : : } // point_src::
688 : :
689 : : namespace gyor {
690 : :
691 : : static std::vector< tk::real >
692 : 0 : ic( tk::real, tk::real, tk::real, tk::real t )
693 : : // *****************************************************************************
694 : : //! Set initial conditions prescribing gyor
695 : : //! \param[in] t Time where to evaluate initial conditions
696 : : //! \return Values of conserved variables
697 : : // *****************************************************************************
698 : : {
699 : 0 : auto r = 1.225;
700 : 0 : auto p = 1.0e5;
701 : :
702 : 0 : auto u = 5.0 * std::cos(t);
703 : 0 : auto v = 5.0 * std::sin(t);
704 : 0 : auto w = 0.0;
705 : 0 : tk::real rE = eos::totalenergy( r, u, v, w, p );
706 : :
707 [ - - ]: 0 : return { r, r*u, r*v, r*w, rE };
708 : : }
709 : :
710 : : } // gyor::
711 : :
712 : : namespace poisson {
713 : :
714 : : static std::vector< tk::real >
715 : 2142 : ic( tk::real, tk::real, tk::real, tk::real )
716 : : // *****************************************************************************
717 : : //! Set velocity initial conditions for testing a Poisson solve only
718 : : //! \return Values for initial conditions
719 : : // *****************************************************************************
720 : : {
721 [ + - ]: 2142 : return { 0, 0, 0 };
722 : : }
723 : :
724 : : } // poisson::
725 : :
726 : : namespace poisson_const {
727 : :
728 : : static tk::real
729 : 990 : pr( tk::real, tk::real, tk::real )
730 : : // *****************************************************************************
731 : : //! Set pressure rhs for testing a Poisson solve
732 : : //! \return Value for pressure rhs
733 : : // *****************************************************************************
734 : : {
735 : 990 : return 6.0;
736 : : }
737 : :
738 : : static tk::real
739 : 3774 : ic( tk::real x, tk::real y, tk::real z )
740 : : // *****************************************************************************
741 : : //! Evaluate pressure boundary condition
742 : : //! \param[in] x X coordinate where to evaluate the BC
743 : : //! \param[in] y Y coordinate where to evaluate the BC
744 : : //! \param[in] z Z coordinate where to evaluate the BC
745 : : //! \return Value for pressure BC
746 : : // *****************************************************************************
747 : : {
748 : 3774 : return x*x + y*y + z*z;
749 : : }
750 : :
751 : : } // poisson_const::
752 : :
753 : : namespace poisson_harmonic {
754 : :
755 : : static tk::real
756 : 0 : pr( tk::real, tk::real, tk::real )
757 : : // *****************************************************************************
758 : : //! Set pressure rhs for testing a Laplace solve
759 : : //! \return Value for pressure rhs
760 : : // *****************************************************************************
761 : : {
762 : 0 : return 0.0;
763 : : }
764 : :
765 : : static tk::real
766 : 0 : ic( tk::real x, tk::real y, tk::real z )
767 : : // *****************************************************************************
768 : : //! Evaluate pressure boundary condition
769 : : //! \param[in] x X coordinate where to evaluate the BC
770 : : //! \param[in] y Y coordinate where to evaluate the BC
771 : : //! \param[in] z Z coordinate where to evaluate the BC
772 : : //! \return Value for pressure BC
773 : : // *****************************************************************************
774 : : {
775 : 0 : const auto& b = g_cfg.get< tag::problem_beta >();
776 : 0 : auto x0 = b[0];
777 : 0 : auto y0 = b[1];
778 : 0 : auto z0 = b[2];
779 : :
780 : 0 : return 1.0 / std::sqrt( (x-x0)*(x-x0) + (y-y0)*(y-y0) + (z-z0)*(z-z0) );
781 : : }
782 : :
783 : : } // poisson_harmonic::
784 : :
785 : : namespace poisson_sine {
786 : :
787 : : static tk::real
788 : 278 : pr( tk::real x, tk::real y, tk::real z )
789 : : // *****************************************************************************
790 : : //! Set pressure rhs for testing a Poisson solve
791 : : //! \return Value for pressure rhs
792 : : // *****************************************************************************
793 : : {
794 : 278 : return -M_PI * M_PI * x * y * std::sin( M_PI * z );
795 : : }
796 : :
797 : : static tk::real
798 : 1067 : ic( tk::real x, tk::real y, tk::real z )
799 : : // *****************************************************************************
800 : : //! Evaluate pressure boundary condition
801 : : //! \param[in] x X coordinate where to evaluate the BC
802 : : //! \param[in] y Y coordinate where to evaluate the BC
803 : : //! \param[in] z Z coordinate where to evaluate the BC
804 : : //! \return Value for pressure BC
805 : : // *****************************************************************************
806 : : {
807 : 1067 : return x * y * std::sin( M_PI * z );
808 : : }
809 : :
810 : : } // poisson_sine::
811 : :
812 : : namespace poisson_sine3 {
813 : :
814 : : static tk::real
815 : 596 : pr( tk::real x, tk::real y, tk::real z )
816 : : // *****************************************************************************
817 : : //! Set pressure rhs for testing a Poisson solve
818 : : //! \return Value for pressure rhs
819 : : // *****************************************************************************
820 : : {
821 : : using std::sin;
822 : :
823 : 596 : return -3.0 * M_PI * M_PI * sin(M_PI*x) * sin(M_PI*y) * sin(M_PI*z);
824 : : }
825 : :
826 : : static tk::real
827 : 2281 : ic( tk::real x, tk::real y, tk::real z )
828 : : // *****************************************************************************
829 : : //! Evaluate pressure boundary condition
830 : : //! \param[in] x X coordinate where to evaluate the BC
831 : : //! \param[in] y Y coordinate where to evaluate the BC
832 : : //! \param[in] z Z coordinate where to evaluate the BC
833 : : //! \return Value for pressure BC
834 : : // *****************************************************************************
835 : : {
836 : 2281 : return sin(M_PI*x) * sin(M_PI*y) * sin(M_PI*z);
837 : : }
838 : :
839 : : } // poisson_sine3::
840 : :
841 : : namespace poisson_neumann {
842 : :
843 : : static tk::real
844 : 278 : pr( tk::real x, tk::real y, tk::real )
845 : : // *****************************************************************************
846 : : //! Set pressure rhs for testing a Poisson solve
847 : : //! \param[in] x X coordinate where to evaluate the rhs
848 : : //! \param[in] y Y coordinate where to evaluate the rhs
849 : : //! \return Value for pressure rhs
850 : : // *****************************************************************************
851 : : {
852 : 278 : return -3.0 * std::cos(2.0*x) * std::exp(y);
853 : : }
854 : :
855 : : static std::array< tk::real, 3 >
856 : 612 : pg( tk::real x, tk::real y, tk::real )
857 : : // *****************************************************************************
858 : : //! Set pressure gradient for testing a Poisson solve
859 : : //! \param[in] x X coordinate where to evaluate the pressure gradient
860 : : //! \param[in] y Y coordinate where to evaluate the pressure gradient
861 : : //! \return Value for pressure gradient at a point
862 : : // *****************************************************************************
863 : : {
864 : 612 : return { -2.0 * std::sin( 2.0 * x ) * std::exp( y ),
865 : 612 : std::cos(2.0*x) * std::exp(y),
866 : 612 : 0.0 };
867 : : }
868 : :
869 : : static tk::real
870 : 972 : ic( tk::real x, tk::real y, tk::real )
871 : : // *****************************************************************************
872 : : //! Evaluate pressure boundary condition
873 : : //! \param[in] x X coordinate where to evaluate the IC / analytic solution
874 : : //! \param[in] y Y coordinate where to evaluate the IC / analytic solution
875 : : //! \return Value for pressure
876 : : // *****************************************************************************
877 : : {
878 : 972 : return std::cos(2.0*x) * std::exp(y);
879 : : }
880 : :
881 : : } // poisson_neumann::
882 : :
883 : :
884 : : namespace poiseuille {
885 : :
886 : : static std::vector< tk::real >
887 : 0 : ic( tk::real, tk::real y, tk::real, tk::real )
888 : : // *****************************************************************************
889 : : //! Set initial conditions prescribing the Poisuille problem
890 : : //! \param[in] y Y coordinate where to evaluate the solution
891 : : //! \return Values of conserved variables
892 : : // *****************************************************************************
893 : : {
894 : 0 : auto eps = std::numeric_limits< tk::real >::epsilon();
895 : 0 : auto nu = g_cfg.get< tag::mat_dyn_viscosity >();
896 [ - - ][ - - ]: 0 : if (nu < eps) Throw( "Poiseuille flow needs nonzero viscosity" );
[ - - ][ - - ]
897 : :
898 : 0 : auto dpdx = -0.12;
899 : 0 : auto u = -dpdx * y * (1.0 - y) / 2.0 / nu;
900 : :
901 [ - - ]: 0 : return { u, 0.0, 0.0 };
902 : : }
903 : :
904 : : static tk::real
905 : 0 : pic( tk::real, tk::real, tk::real )
906 : : // *****************************************************************************
907 : : //! Set homogeneous initial conditions for Poiseuille
908 : : //! \return Value of pressure
909 : : // *****************************************************************************
910 : : {
911 : 0 : return 0.0;
912 : : }
913 : :
914 : : } // poiseuille::
915 : :
916 : :
917 : : std::function< std::vector< tk::real >
918 : : ( tk::real, tk::real, tk::real, tk::real ) >
919 : 54822 : IC()
920 : : // *****************************************************************************
921 : : // Query user config and assign function to set initial conditions
922 : : //! \return The function to call to set initial conditions
923 : : // *****************************************************************************
924 : : {
925 : 54822 : const auto& problem = inciter::g_cfg.get< tag::problem >();
926 : :
927 : : std::function< std::vector< tk::real >
928 : 54822 : ( tk::real, tk::real, tk::real, tk::real ) > ic;
929 : :
930 [ + + ]: 54822 : if (problem == "userdef")
931 : 5341 : ic = userdef::ic;
932 [ + + ]: 49481 : else if (problem == "nonlinear_energy_growth")
933 : 1254 : ic = nonlinear_energy_growth::ic;
934 [ + + ]: 48227 : else if (problem == "rayleigh_taylor")
935 : 3624 : ic = rayleigh_taylor::ic;
936 [ + + ]: 44603 : else if (problem == "sedov")
937 : 59 : ic = sedov::ic;
938 [ + + ]: 44544 : else if (problem == "sod")
939 : 629 : ic = sod::ic;
940 [ + + ]: 43915 : else if (problem == "taylor_green")
941 : 4812 : ic = taylor_green::ic;
942 [ + + ]: 39103 : else if (problem == "vortical_flow")
943 : 29151 : ic = vortical_flow::ic;
944 [ + + ]: 9952 : else if (problem == "slot_cyl")
945 : 9916 : ic = slot_cyl::ic;
946 [ + + ]: 36 : else if (problem == "point_src")
947 : 4 : ic = point_src::ic;
948 [ - + ]: 32 : else if (problem == "gyor")
949 : 0 : ic = gyor::ic;
950 [ + - ]: 32 : else if (problem.find("poisson") != std::string::npos)
951 : 32 : ic = poisson::ic;
952 [ - - ]: 0 : else if (problem == "poiseuille")
953 : 0 : ic = poiseuille::ic;
954 : : else
955 [ - - ][ - - ]: 0 : Throw( "problem type ic not hooked up" );
[ - - ]
956 : :
957 : 54822 : return ic;
958 : 0 : }
959 : :
960 : : std::function< std::vector< tk::real >
961 : : ( tk::real, tk::real, tk::real, tk::real ) >
962 : 36184 : SOL()
963 : : // *****************************************************************************
964 : : // Query user config and assign function to query analytic solutions
965 : : //! \return The function to call to query analytic solutions
966 : : // *****************************************************************************
967 : : {
968 : 36184 : const auto& problem = inciter::g_cfg.get< tag::problem >();
969 : :
970 [ + + ]: 58170 : if (problem == "userdef" ||
971 [ + + ]: 41557 : problem == "sod" ||
972 [ + + ][ + + ]: 77741 : problem == "sedov" ||
[ + + ]
973 : 17557 : problem == "point_src")
974 : 18683 : return {};
975 : : else
976 : 17501 : return IC();
977 : : }
978 : :
979 : : void
980 : 2671 : initialize( const std::array< std::vector< tk::real >, 3 >& coord,
981 : : tk::Fields& U,
982 : : tk::real t,
983 : : const std::vector< std::unordered_set< std::size_t > >& boxnodes )
984 : : // *****************************************************************************
985 : : // Set inital conditions
986 : : //! \param[in] coord Mesh node coordinates
987 : : //! \param[in,out] U Array of unknowns
988 : : //! \param[in] t Physical time
989 : : //! \param[in] boxnodes Nodes at which box user ICs are set (for each box IC)
990 : : // *****************************************************************************
991 : : {
992 [ - + ][ - - ]: 2671 : Assert( coord[0].size() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
993 : :
994 [ + - ]: 2671 : auto ic = IC();
995 : 2671 : const auto& x = coord[0];
996 : 2671 : const auto& y = coord[1];
997 : 2671 : const auto& z = coord[2];
998 : :
999 : : // Set initial conditions dependeing on problem configured
1000 [ + + ]: 347744 : for (std::size_t i=0; i<x.size(); ++i) {
1001 : :
1002 : : // Set background ICs
1003 [ + - ]: 345073 : auto s = ic( x[i], y[i], z[i], t );
1004 [ - + ][ - - ]: 345073 : Assert( s.size() == U.nprop(), "Size mismatch" );
[ - - ][ - - ]
1005 : :
1006 : : // Initialize user-defined ICs in boxes
1007 [ + - ]: 345073 : box( i, s, boxnodes );
1008 : :
1009 : : // Set values for ICs
1010 [ + - ][ + + ]: 2104547 : for (std::size_t c=0; c<s.size(); ++c) U(i,c) = s[c];
1011 : :
1012 : 345073 : }
1013 : 2671 : }
1014 : :
1015 : : std::function< tk::real( tk::real, tk::real, tk::real ) >
1016 : 4978 : PRESSURE_RHS()
1017 : : // *****************************************************************************
1018 : : // Query user config and assign function to set pressure rhs
1019 : : //! \return The function to call to set pressure rhs
1020 : : // *****************************************************************************
1021 : : {
1022 : 4978 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1023 : :
1024 : 4978 : std::function< tk::real( tk::real, tk::real, tk::real ) > pr;
1025 : :
1026 [ + + ]: 4978 : if (problem == "poisson_const")
1027 : 22 : pr = poisson_const::pr;
1028 [ - + ]: 4956 : else if (problem == "poisson_harmonic")
1029 : 0 : pr = poisson_harmonic::pr;
1030 [ + + ]: 4956 : else if (problem == "poisson_sine")
1031 : 2 : pr = poisson_sine::pr;
1032 [ + + ]: 4954 : else if (problem == "poisson_sine3")
1033 : 6 : pr = poisson_sine3::pr;
1034 [ + + ]: 4948 : else if (problem == "poisson_neumann")
1035 : 2 : pr = poisson_neumann::pr;
1036 : :
1037 : 4978 : return pr;
1038 : : }
1039 : :
1040 : : std::function< tk::real( tk::real, tk::real, tk::real ) >
1041 : 2432 : PRESSURE_IC()
1042 : : // *****************************************************************************
1043 : : // Query user config and assign function to set pressure initial conditions
1044 : : //! \return The function to call to set pressure initial conditions
1045 : : // *****************************************************************************
1046 : : {
1047 : 2432 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1048 : :
1049 : 2432 : std::function< tk::real( tk::real, tk::real, tk::real ) > ic;
1050 : :
1051 [ + + ]: 2432 : if (problem == "userdef")
1052 : 2290 : ic = userdef::pic;
1053 [ + + ]: 142 : else if (problem == "poisson_const")
1054 : 94 : ic = poisson_const::ic;
1055 [ - + ]: 48 : else if (problem == "poisson_harmonic")
1056 : 0 : ic = poisson_harmonic::ic;
1057 [ + + ]: 48 : else if (problem == "poisson_sine")
1058 : 10 : ic = poisson_sine::ic;
1059 [ + + ]: 38 : else if (problem == "poisson_sine3")
1060 : 28 : ic = poisson_sine3::ic;
1061 [ + - ]: 10 : else if (problem == "poisson_neumann")
1062 : 10 : ic = poisson_neumann::ic;
1063 [ - - ]: 0 : else if (problem == "poiseuille")
1064 : 0 : ic = poiseuille::pic;
1065 : : else
1066 [ - - ][ - - ]: 0 : Throw( "problem type not hooked up" );
[ - - ]
1067 : :
1068 : 2432 : return ic;
1069 : 0 : }
1070 : :
1071 : : std::function< tk::real( tk::real, tk::real, tk::real ) >
1072 : 4430 : PRESSURE_SOL()
1073 : : // *****************************************************************************
1074 : : // Query user config and assign function to query analytic pressure solutions
1075 : : //! \return The function to call to query analytic solutions
1076 : : // *****************************************************************************
1077 : : {
1078 : 4430 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1079 : :
1080 [ + + ][ - + ]: 4430 : if (problem == "userdef" || problem == "poiseuille")
[ + + ]
1081 : 4320 : return {};
1082 : : else
1083 : 110 : return PRESSURE_IC();
1084 : : }
1085 : :
1086 : : std::function< std::array< tk::real, 3 >( tk::real, tk::real, tk::real ) >
1087 : 4978 : PRESSURE_GRAD()
1088 : : // *****************************************************************************
1089 : : // Assign function to query pressure gradient at a point
1090 : : //! \return The function to call to query the pressure gradient
1091 : : // *****************************************************************************
1092 : : {
1093 : 4978 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1094 : :
1095 [ + + ]: 4978 : if (problem == "poisson_neumann")
1096 : 2 : return poisson_neumann::pg;
1097 : :
1098 : 4976 : return {};
1099 : : }
1100 : :
1101 : : tk::real
1102 : 0 : initialize( tk::real x, tk::real y, tk::real z )
1103 : : // *****************************************************************************
1104 : : // Evaluate initial condition for pressure
1105 : : //! \param[in] x X coordinate where to evaluate the pressure initial condition
1106 : : //! \param[in] y Y coordinate where to evaluate the pressure initial condition
1107 : : //! \param[in] z Z coordinate where to evaluate the pressure initial condition
1108 : : //! \return Pressure initial condition
1109 : : // *****************************************************************************
1110 : : {
1111 : 0 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1112 : :
1113 : 0 : std::function< tk::real( tk::real, tk::real, tk::real ) > ic;
1114 : :
1115 [ - - ]: 0 : if (problem == "poisson_const")
1116 : 0 : ic = poisson_const::ic;
1117 [ - - ]: 0 : else if (problem == "poisson_harmonic")
1118 : 0 : ic = poisson_harmonic::ic;
1119 [ - - ]: 0 : else if (problem == "poisson_sine")
1120 : 0 : ic = poisson_sine::ic;
1121 [ - - ]: 0 : else if (problem == "poisson_sine3")
1122 : 0 : ic = poisson_sine3::ic;
1123 : : else
1124 [ - - ][ - - ]: 0 : Throw( "problem type not hooked up" );
[ - - ]
1125 : :
1126 [ - - ]: 0 : return ic( x, y, z );
1127 : 0 : }
1128 : :
1129 : : std::function< std::vector< tk::real >
1130 : : ( tk::real, tk::real, tk::real, tk::real ) >
1131 : 62611 : SRC()
1132 : : // *****************************************************************************
1133 : : // Query user config and assign function to add a source term
1134 : : //! \return The function to call to evaluate a problem-sepcific source term
1135 : : // *****************************************************************************
1136 : : {
1137 : 62611 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1138 : :
1139 : : std::function<
1140 : 62611 : std::vector< tk::real >( tk::real, tk::real, tk::real, tk::real ) > src;
1141 : :
1142 [ + + ]: 62611 : if (problem == "nonlinear_energy_growth")
1143 : 676 : src = nonlinear_energy_growth::src;
1144 [ + + ]: 61935 : else if (problem == "rayleigh_taylor")
1145 : 2000 : src = rayleigh_taylor::src;
1146 [ + + ]: 59935 : else if (problem == "taylor_green")
1147 : 3536 : src = taylor_green::src;
1148 [ + + ]: 56399 : else if (problem == "vortical_flow")
1149 : 18547 : src = vortical_flow::src;
1150 [ + + ]: 37852 : else if (problem == "slot_cyl")
1151 : 5497 : src = slot_cyl::src;
1152 : :
1153 : 62611 : return src;
1154 : : }
1155 : :
1156 : : std::function< void( const std::array< std::vector< tk::real >, 3 >&,
1157 : : tk::real,
1158 : : tk::Fields& ) >
1159 : 71021 : PHYS_SRC()
1160 : : // *****************************************************************************
1161 : : // Query user config and assign function to apply source to numerical solution
1162 : : //! \return The function to call to evaluate a problem-sepcific source term
1163 : : // *****************************************************************************
1164 : : {
1165 : 71021 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1166 : :
1167 : : std::function< void( const std::array< std::vector< tk::real >, 3 >&,
1168 : : tk::real,
1169 : 71021 : tk::Fields& ) > src;
1170 : :
1171 [ + + ]: 71021 : if (problem == "point_src") {
1172 : 600 : src = point_src::src;
1173 : : }
1174 : :
1175 : 71021 : return src;
1176 : : }
1177 : :
1178 : : } // problems::
|
