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-2025 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 : 1276019 : 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 : 1276019 : const auto& solver = g_cfg.get< tag::solver >();
40 [ + + ]: 1276019 : if (solver == "chocg") {
41 : 659594 : const auto& ncomp = g_cfg.get< tag::problem_ncomp >();
42 [ + - ]: 659594 : std::vector< tk::real > u( ncomp, 0.0 );
43 [ + - ]: 659594 : auto ic_velocity = g_cfg.get< tag::ic_velocity >();
44 : 659594 : auto large = std::numeric_limits< double >::max();
45 [ + - ]: 659594 : if (std::abs(ic_velocity[0] - large) > 1.0e-12) u[0] = ic_velocity[0];
46 [ + - ]: 659594 : if (std::abs(ic_velocity[1] - large) > 1.0e-12) u[1] = ic_velocity[1];
47 [ + - ]: 659594 : if (std::abs(ic_velocity[2] - large) > 1.0e-12) u[2] = ic_velocity[2];
48 : 659594 : return u;
49 : 659594 : }
50 [ + + ]: 616425 : else if (solver == "lohcg") {
51 : 475985 : const auto& ncomp = g_cfg.get< tag::problem_ncomp >();
52 [ + - ]: 475985 : std::vector< tk::real > u( ncomp, 0.0 );
53 [ + - ]: 475985 : auto ic_velocity = g_cfg.get< tag::ic_velocity >();
54 : 475985 : auto large = std::numeric_limits< double >::max();
55 [ + - ]: 475985 : if (std::abs(ic_velocity[0] - large) > 1.0e-12) u[1] = ic_velocity[0];
56 [ + - ]: 475985 : if (std::abs(ic_velocity[1] - large) > 1.0e-12) u[2] = ic_velocity[1];
57 [ + - ]: 475985 : if (std::abs(ic_velocity[2] - large) > 1.0e-12) u[3] = ic_velocity[2];
58 : 475985 : return u;
59 : 475985 : }
60 : :
61 : : // density-based solvers
62 : :
63 : 140440 : auto ic_density = g_cfg.get< tag::ic_density >();
64 : 140440 : const auto& ic_velocity = g_cfg.get< tag::ic_velocity >();
65 [ - + ][ - - ]: 140440 : ErrChk( ic_velocity.size() == 3, "ic_velocity must have 3 components" );
[ - - ][ - - ]
66 : :
67 [ + - ]: 140440 : std::vector< tk::real > u( 5, 0.0 );
68 : :
69 : 140440 : u[0] = ic_density;
70 : 140440 : u[1] = u[0] * ic_velocity[0];
71 : 140440 : u[2] = u[0] * ic_velocity[1];
72 : 140440 : u[3] = u[0] * ic_velocity[2];
73 : :
74 : 140440 : auto ic_pressure = g_cfg.get< tag::ic_pressure >();
75 : 140440 : auto ic_energy = g_cfg.get< tag::ic_energy >();
76 : 140440 : auto ic_temperature = g_cfg.get< tag::ic_temperature >();
77 : :
78 : 140440 : auto largereal = std::numeric_limits< double >::max();
79 : :
80 [ + - ]: 140440 : if (std::abs(ic_pressure - largereal) > 1.0e-12) {
81 : :
82 : 140440 : 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 : 140440 : return u;
104 : 140440 : }
105 : :
106 : : static tk::real
107 : 8084 : 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 : 8084 : 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 : 66650 : 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 : 66650 : auto eps = std::numeric_limits< tk::real >::epsilon();
343 : : tk::real p;
344 [ + + ][ + + ]: 66650 : if (abs(x) < eps && abs(y) < eps && abs(z) < eps) {
[ + + ][ + + ]
345 : 5 : p = g_cfg.get< tag::problem_p0 >();
346 : : } else {
347 : 66645 : p = 0.67e-4;
348 : : }
349 : :
350 : : // density
351 : 66650 : tk::real r = 1.0;
352 : : // velocity
353 : 66650 : tk::real u = 0.0;
354 : 66650 : tk::real v = 0.0;
355 : 66650 : tk::real w = 0.0;
356 : : // total specific energy
357 : 66650 : tk::real rE = eos::totalenergy( r, u, v, w, p );
358 : :
359 [ + - ]: 66650 : 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 : 6179899 : 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 : 6179899 : tk::real p0 = 1.0;
519 : :
520 : : // configure number of scalar components
521 : 6179899 : std::size_t ncomp = 6;
522 : 6179899 : const auto& solver = g_cfg.get< tag::solver >();
523 [ + + ]: 6179899 : if (solver == "chocg") {
524 : 997120 : ncomp = 4;
525 : : }
526 [ + + ]: 5182779 : else if (solver == "lohcg") {
527 : 882080 : ncomp = 5;
528 : : }
529 : :
530 [ + - ]: 6179899 : std::vector< tk::real > u( ncomp, 0.0 );
531 : :
532 : : // prescribed velocity: rotate in x-y plane
533 : 6179899 : std::size_t sc = 5;
534 [ + + ]: 6179899 : if (solver == "chocg") {
535 : 997120 : u[0] = 0.5 - y;
536 : 997120 : u[1] = x - 0.5;
537 : 997120 : u[2] = 0.0;
538 : 997120 : sc = 3;
539 : : }
540 [ + + ]: 5182779 : else if (solver == "lohcg") {
541 : 882080 : u[0] = 0.0;
542 : 882080 : u[1] = 0.5 - y;
543 : 882080 : u[2] = x - 0.5;
544 : 882080 : u[3] = 0.0;
545 : 882080 : sc = 4;
546 : : }
547 : : else {
548 : 4300699 : u[0] = 1.0;
549 : 4300699 : u[1] = u[0] * (0.5 - y);
550 : 4300699 : u[2] = u[0] * (x - 0.5);
551 : 4300699 : u[3] = 0.0;
552 : 4300699 : u[4] = eos::totalenergy( u[0], u[1]/u[0], u[2]/u[0], u[3]/u[0], p0 );
553 : : }
554 : :
555 : 6179899 : const tk::real R0 = 0.15;
556 : :
557 : : // center of the cone
558 : 6179899 : tk::real x0 = 0.5;
559 : 6179899 : tk::real y0 = 0.25;
560 : 6179899 : tk::real r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
561 : 6179899 : tk::real kx = 0.5 + r*sin( t );
562 : 6179899 : tk::real ky = 0.5 - r*cos( t );
563 : :
564 : : // center of the hump
565 : 6179899 : x0 = 0.25;
566 : 6179899 : y0 = 0.5;
567 : 6179899 : r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
568 : 6179899 : tk::real hx = 0.5 + r*sin( t-M_PI/2.0 ),
569 : 6179899 : hy = 0.5 - r*cos( t-M_PI/2.0 );
570 : :
571 : : // center of the slotted cylinder
572 : 6179899 : x0 = 0.5;
573 : 6179899 : y0 = 0.75;
574 : 6179899 : r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
575 : 6179899 : tk::real cx = 0.5 + r*sin( t+M_PI ),
576 : 6179899 : cy = 0.5 - r*cos( t+M_PI );
577 : :
578 : : // end points of the cylinder slot
579 : 6179899 : tk::real i1x = 0.525, i1y = cy - r*cos( std::asin(0.025/r) ),
580 : 6179899 : i2x = 0.525, i2y = 0.8,
581 : 6179899 : i3x = 0.475, i3y = 0.8;
582 : :
583 : : // rotate end points of cylinder slot
584 : 6179899 : tk::real ri1x = 0.5 + cos(t)*(i1x-0.5) - sin(t)*(i1y-0.5),
585 : 6179899 : ri1y = 0.5 + sin(t)*(i1x-0.5) + cos(t)*(i1y-0.5),
586 : 6179899 : ri2x = 0.5 + cos(t)*(i2x-0.5) - sin(t)*(i2y-0.5),
587 : 6179899 : ri2y = 0.5 + sin(t)*(i2x-0.5) + cos(t)*(i2y-0.5),
588 : 6179899 : ri3x = 0.5 + cos(t)*(i3x-0.5) - sin(t)*(i3y-0.5),
589 : 6179899 : ri3y = 0.5 + sin(t)*(i3x-0.5) + cos(t)*(i3y-0.5);
590 : :
591 : : // direction of slot sides
592 : 6179899 : tk::real v1x = ri2x-ri1x, v1y = ri2y-ri1y,
593 : 6179899 : v2x = ri3x-ri2x, v2y = ri3y-ri2y;
594 : :
595 : : // lengths of direction of slot sides vectors
596 : 6179899 : tk::real v1 = sqrt(v1x*v1x + v1y*v1y),
597 : 6179899 : v2 = sqrt(v2x*v2x + v2y*v2y);
598 : :
599 : : // cone
600 : 6179899 : r = sqrt((x-kx)*(x-kx) + (y-ky)*(y-ky)) / R0;
601 [ + + ]: 6179899 : if (r<1.0) u[sc] = 0.6*(1.0-r);
602 : :
603 : : // hump
604 : 6179899 : r = sqrt((x-hx)*(x-hx) + (y-hy)*(y-hy)) / R0;
605 [ + + ]: 6179899 : if (r<1.0) u[sc] = 0.2*(1.0+cos(M_PI*std::min(r,1.0)));
606 : :
607 : : // cylinder
608 : 6179899 : r = sqrt((x-cx)*(x-cx) + (y-cy)*(y-cy)) / R0;
609 : 6179899 : const std::array< tk::real, 2 > r1{{ v1x, v1y }},
610 : 6179899 : r2{{ x-ri1x, y-ri1y }};
611 : 6179899 : const auto d1 = (r1[0]*r2[1] - r2[0]*r1[1]) / v1;
612 : 6179899 : const std::array< tk::real, 2 > r3{{ v2x, v2y }},
613 : 6179899 : r4{{ x-ri2x, y-ri2y }};
614 : 6179899 : const auto d2 = (r3[0]*r4[1] - r4[0]*r3[1]) / v2;
615 [ + + ][ + + ]: 6179899 : if (r<1.0 && (d1>0.05 || d1<0.0 || d2<0.0)) u[sc] = 0.6;
[ + + ][ + + ]
616 : :
617 : 12359798 : return u;
618 : : }
619 : :
620 : : static std::vector< tk::real >
621 : 3681250 : src( tk::real x, tk::real y, tk::real z, tk::real t )
622 : : // *****************************************************************************
623 : : //! Compute and return source term for slotted cylinder, cone, Gauss hump
624 : : //! \param[in] x X coordinate where to evaluate the source
625 : : //! \param[in] y Y coordinate where to evaluate the source
626 : : //! \param[in] z Z coordinate where to evaluate the source
627 : : //! \param[in] t Time where to evaluate the source
628 : : //! \return Source for flow variables + transported scalars
629 : : // *****************************************************************************
630 : : {
631 : : // evaluate solution at x,y,z,t
632 [ + - ]: 3681250 : auto u = ic( x, y, z, t );
633 : :
634 : : // configure number of scalar components
635 : 3681250 : std::size_t ncomp = 6;
636 : 3681250 : const auto& solver = g_cfg.get< tag::solver >();
637 [ + + ]: 3681250 : if (solver == "chocg") {
638 : 369140 : ncomp = 4;
639 : : }
640 [ + + ]: 3312110 : else if (solver == "lohcg") {
641 : 369140 : ncomp = 5;
642 : : }
643 : :
644 [ + - ]: 3681250 : std::vector< tk::real > s( ncomp, 0.0 );
645 : :
646 : : // momentum source
647 [ + + ]: 3681250 : if (solver == "chocg") {
648 : 369140 : s[0] = -u[1];
649 : 369140 : s[1] = u[0];
650 : : }
651 [ + + ]: 3312110 : else if (solver == "lohcg") {
652 : 369140 : s[1] = -u[2];
653 : 369140 : s[2] = u[1];
654 : : }
655 : : else {
656 : 2942970 : s[1] = -u[2];
657 : 2942970 : s[2] = u[1];
658 : : }
659 : :
660 : 7362500 : return s;
661 : 3681250 : }
662 : :
663 : : } // slot_cyl::
664 : :
665 : : namespace sheardiff {
666 : :
667 : : static std::vector< tk::real >
668 : 0 : ic( tk::real x, tk::real y, tk::real, tk::real t )
669 : : // *****************************************************************************
670 : : //! Set initial conditions prescribing shear-diffusion in 2D
671 : : //! \param[in] x X coordinate where to evaluate the solution
672 : : //! \param[in] y Y coordinate where to evaluate the solution
673 : : //! \param[in] t Time where to evaluate the solution
674 : : //! \return Values of conserved variables
675 : : //! \see A. Okubo, M.J. Karweit, Diffusion from a continuous source in a uniform
676 : : //! shear flow, Limnology and Oceanography, 14, 1969,
677 : : //! https://doi.org/10.4319/lo.1969.14.4.0514.
678 : : // *****************************************************************************
679 : : {
680 : : using std::exp;
681 : : using std::pow;
682 : : using std::sqrt;
683 : :
684 : : // manufactured solution parameters
685 : 0 : auto V0 = g_cfg.get< tag::problem_p0 >(); // translation velocity in x
686 : 0 : auto L = g_cfg.get< tag::problem_alpha >();// shear velocity
687 : 0 : auto t0 = g_cfg.get< tag::t0 >(); // initial time
688 : 0 : auto dif = g_cfg.get< tag::mat_dyn_diffusivity >(); // diffusivity
689 : 0 : auto eps = std::numeric_limits< tk::real >::epsilon();
690 [ - - ][ - - ]: 0 : if (dif < eps) Throw( "Diffusivity must be positive" );
[ - - ][ - - ]
691 : :
692 : : // configure number of scalar components
693 : 0 : std::size_t ncomp = 4;
694 : 0 : const auto& solver = g_cfg.get< tag::solver >();
695 [ - - ]: 0 : if (solver == "chocg") {
696 : : } else
697 [ - - ]: 0 : if (solver == "lohcg") {
698 : 0 : ncomp = 5;
699 : : }
700 [ - - ][ - - ]: 0 : else Throw( "Shear-diff IC not setup for this solver" );
[ - - ]
701 : :
702 : : // prescribed velocity
703 [ - - ]: 0 : std::vector< tk::real > u( ncomp, 0.0 );
704 : 0 : std::size_t sc = 3;
705 [ - - ]: 0 : if (solver == "chocg") {
706 : 0 : u[0] = V0 + L*y;
707 : 0 : u[1] = u[2] = 0.0;
708 : : }
709 [ - - ]: 0 : else if (solver == "lohcg") {
710 : 0 : u[1] = V0 + L*y;
711 : 0 : u[2] = u[3] = 0.0;
712 : 0 : sc = 4;
713 : : }
714 : :
715 : 0 : auto M = [=](tk::real T){ return 4.0*M_PI*T*sqrt(1.0 + L*L*T*T/12.0); };
716 : :
717 : 0 : u[sc] = M(t0) / M(t) *
718 : 0 : exp( -pow( x - V0*t - 0.5*L*y*t, 2.0 ) /
719 : 0 : (4.0*dif*t*(1.0 + L*L*t*t/12.0)) - y*y / (4.0*dif*t) );
720 : :
721 : 0 : return u;
722 : : }
723 : :
724 : : static std::vector< tk::real >
725 : 0 : src( tk::real, tk::real, tk::real, tk::real )
726 : : // *****************************************************************************
727 : : //! Compute and return source term for shear-diffusion in 2D
728 : : //! \return Source for flow variables + transported scalars
729 : : //! \see A. Okubo, M.J. Karweit, Diffusion from a continuous source in a uniform
730 : : //! shear flow, Limnology and Oceanography, 14, 1969,
731 : : //! https://doi.org/10.4319/lo.1969.14.4.0514.
732 : : // *****************************************************************************
733 : : {
734 : : // configure number of scalar components
735 : 0 : std::size_t ncomp = 4;
736 : 0 : const auto& solver = g_cfg.get< tag::solver >();
737 [ - - ]: 0 : if (solver == "lohcg") {
738 : 0 : ncomp = 5;
739 : : }
740 : :
741 : : // manufactured solution parameters
742 : 0 : auto L = g_cfg.get< tag::problem_alpha >();// shear velocity
743 : :
744 : : // source
745 [ - - ]: 0 : std::vector< tk::real > s( ncomp, 0.0 );
746 [ - - ]: 0 : if (solver == "lohcg") {
747 : 0 : s[0] = -L;
748 : : }
749 : :
750 : 0 : return s;
751 : : }
752 : :
753 : : } // sheardiff::
754 : :
755 : : namespace point_src {
756 : :
757 : : static std::vector< tk::real >
758 : 1858 : ic( tk::real x, tk::real y, tk::real z, tk::real t )
759 : : // *****************************************************************************
760 : : //! Set initial conditions for point source problem
761 : : //! \param[in] x X coordinate where to evaluate initial conditions
762 : : //! \param[in] y Y coordinate where to evaluate initial conditions
763 : : //! \param[in] z Z coordinate where to evaluate initial conditions
764 : : //! \param[in] t Time where to evaluate initial conditions
765 : : //! \return Values of conserved variables
766 : : // *****************************************************************************
767 : : {
768 : 1858 : auto u = userdef::ic( x, y, z, t );
769 [ + - ]: 1858 : u.push_back( 0.0 );
770 : 1858 : return u;
771 : 0 : }
772 : :
773 : : static void
774 : 600 : src( const std::array< std::vector< tk::real >, 3 >& coord,
775 : : tk::real t,
776 : : tk::Fields& U )
777 : : // *****************************************************************************
778 : : //! Apply point-source directly to numerical solution
779 : : //! \param[in] coord Mesh node coordinates
780 : : //! \param[in] t Physical time
781 : : //! \param[in,out] U Solution vector at recent time step
782 : : //! \note This is different from other source terms, because this directly
783 : : //! modifies the solution instead of applied as a source term mathematically.
784 : : //! Hence the function signature is also different.
785 : : // *****************************************************************************
786 : : {
787 [ - + ]: 600 : if (U.nprop() == 5) return;
788 : :
789 : 600 : const auto& source = g_cfg.get< tag::problem_src >();
790 : 600 : const auto& location = source.get< tag::location >();
791 : 600 : auto radius = source.get< tag::radius >();
792 : 600 : auto release_time = source.get< tag::release_time >();
793 : 600 : auto largereal = std::numeric_limits< double >::max();
794 : :
795 [ + - ]: 1200 : if (location.size() != 3 ||
796 [ + - ][ - + ]: 1200 : std::abs(radius - largereal) < 1.0e-12 ||
[ - + ]
797 : 600 : std::abs(release_time - largereal) < 1.0e-12)
798 : : {
799 : 0 : return;
800 : : }
801 : :
802 : 600 : auto sx = location[0];
803 : 600 : auto sy = location[1];
804 : 600 : auto sz = location[2];
805 : 600 : auto sr = radius;
806 : 600 : auto st = release_time;
807 : :
808 [ - + ]: 600 : if (t < st) return;
809 : :
810 : 600 : const auto& x = coord[0];
811 : 600 : const auto& y = coord[1];
812 : 600 : const auto& z = coord[2];
813 : :
814 [ + + ]: 279300 : for (std::size_t i=0; i<U.nunk(); ++i) {
815 : 278700 : auto rx = sx - x[i];
816 : 278700 : auto ry = sy - y[i];
817 : 278700 : auto rz = sz - z[i];
818 [ + + ]: 278700 : if (rx*rx + ry*ry + rz*rz < sr*sr) U(i,5) = 1.0;
819 : : }
820 : :
821 : 600 : return;
822 : : }
823 : :
824 : : } // point_src::
825 : :
826 : : namespace poisson {
827 : :
828 : : static std::vector< tk::real >
829 : 8568 : ic( tk::real, tk::real, tk::real, tk::real )
830 : : // *****************************************************************************
831 : : //! Set velocity initial conditions for testing a Poisson solve only
832 : : //! \return Values for initial conditions
833 : : // *****************************************************************************
834 : : {
835 [ + - ]: 8568 : return { 0, 0, 0 };
836 : : }
837 : :
838 : : } // poisson::
839 : :
840 : : namespace poisson_const {
841 : :
842 : : static tk::real
843 : 990 : pr( tk::real, tk::real, tk::real )
844 : : // *****************************************************************************
845 : : //! Set pressure rhs for testing a Poisson solve
846 : : //! \return Value for pressure rhs
847 : : // *****************************************************************************
848 : : {
849 : 990 : return 6.0;
850 : : }
851 : :
852 : : static tk::real
853 : 3781 : ic( tk::real x, tk::real y, tk::real z )
854 : : // *****************************************************************************
855 : : //! Evaluate pressure boundary condition
856 : : //! \param[in] x X coordinate where to evaluate the BC
857 : : //! \param[in] y Y coordinate where to evaluate the BC
858 : : //! \param[in] z Z coordinate where to evaluate the BC
859 : : //! \return Value for pressure BC
860 : : // *****************************************************************************
861 : : {
862 : 3781 : return x*x + y*y + z*z;
863 : : }
864 : :
865 : : } // poisson_const::
866 : :
867 : : namespace poisson_harmonic {
868 : :
869 : : static tk::real
870 : 0 : pr( tk::real, tk::real, tk::real )
871 : : // *****************************************************************************
872 : : //! Set pressure rhs for testing a Laplace solve
873 : : //! \return Value for pressure rhs
874 : : // *****************************************************************************
875 : : {
876 : 0 : return 0.0;
877 : : }
878 : :
879 : : static tk::real
880 : 0 : ic( tk::real x, tk::real y, tk::real z )
881 : : // *****************************************************************************
882 : : //! Evaluate pressure boundary condition
883 : : //! \param[in] x X coordinate where to evaluate the BC
884 : : //! \param[in] y Y coordinate where to evaluate the BC
885 : : //! \param[in] z Z coordinate where to evaluate the BC
886 : : //! \return Value for pressure BC
887 : : // *****************************************************************************
888 : : {
889 : 0 : const auto& b = g_cfg.get< tag::problem_beta >();
890 : 0 : auto x0 = b[0];
891 : 0 : auto y0 = b[1];
892 : 0 : auto z0 = b[2];
893 : :
894 : 0 : return 1.0 / std::sqrt( (x-x0)*(x-x0) + (y-y0)*(y-y0) + (z-z0)*(z-z0) );
895 : : }
896 : :
897 : : } // poisson_harmonic::
898 : :
899 : : namespace poisson_sine {
900 : :
901 : : static tk::real
902 : 278 : pr( tk::real x, tk::real y, tk::real z )
903 : : // *****************************************************************************
904 : : //! Set pressure rhs for testing a Poisson solve
905 : : //! \return Value for pressure rhs
906 : : // *****************************************************************************
907 : : {
908 : 278 : return -M_PI * M_PI * x * y * std::sin( M_PI * z );
909 : : }
910 : :
911 : : static tk::real
912 : 1067 : ic( tk::real x, tk::real y, tk::real z )
913 : : // *****************************************************************************
914 : : //! Evaluate pressure boundary condition
915 : : //! \param[in] x X coordinate where to evaluate the BC
916 : : //! \param[in] y Y coordinate where to evaluate the BC
917 : : //! \param[in] z Z coordinate where to evaluate the BC
918 : : //! \return Value for pressure BC
919 : : // *****************************************************************************
920 : : {
921 : 1067 : return x * y * std::sin( M_PI * z );
922 : : }
923 : :
924 : : } // poisson_sine::
925 : :
926 : : namespace poisson_sine3 {
927 : :
928 : : static tk::real
929 : 596 : pr( tk::real x, tk::real y, tk::real z )
930 : : // *****************************************************************************
931 : : //! Set pressure rhs for testing a Poisson solve
932 : : //! \return Value for pressure rhs
933 : : // *****************************************************************************
934 : : {
935 : : using std::sin;
936 : :
937 : 596 : return -3.0 * M_PI * M_PI * sin(M_PI*x) * sin(M_PI*y) * sin(M_PI*z);
938 : : }
939 : :
940 : : static tk::real
941 : 2281 : ic( tk::real x, tk::real y, tk::real z )
942 : : // *****************************************************************************
943 : : //! Evaluate pressure boundary condition
944 : : //! \param[in] x X coordinate where to evaluate the BC
945 : : //! \param[in] y Y coordinate where to evaluate the BC
946 : : //! \param[in] z Z coordinate where to evaluate the BC
947 : : //! \return Value for pressure BC
948 : : // *****************************************************************************
949 : : {
950 : 2281 : return sin(M_PI*x) * sin(M_PI*y) * sin(M_PI*z);
951 : : }
952 : :
953 : : } // poisson_sine3::
954 : :
955 : : namespace poisson_neumann {
956 : :
957 : : static tk::real
958 : 278 : pr( tk::real x, tk::real y, tk::real )
959 : : // *****************************************************************************
960 : : //! Set pressure rhs for testing a Poisson solve
961 : : //! \param[in] x X coordinate where to evaluate the rhs
962 : : //! \param[in] y Y coordinate where to evaluate the rhs
963 : : //! \return Value for pressure rhs
964 : : // *****************************************************************************
965 : : {
966 : 278 : return -3.0 * std::cos(2.0*x) * std::exp(y);
967 : : }
968 : :
969 : : static std::array< tk::real, 3 >
970 : 612 : pg( tk::real x, tk::real y, tk::real )
971 : : // *****************************************************************************
972 : : //! Set pressure gradient for testing a Poisson solve
973 : : //! \param[in] x X coordinate where to evaluate the pressure gradient
974 : : //! \param[in] y Y coordinate where to evaluate the pressure gradient
975 : : //! \return Value for pressure gradient at a point
976 : : // *****************************************************************************
977 : : {
978 : 612 : return { -2.0 * std::sin( 2.0 * x ) * std::exp( y ),
979 : 612 : std::cos(2.0*x) * std::exp(y),
980 : 612 : 0.0 };
981 : : }
982 : :
983 : : static tk::real
984 : 972 : ic( tk::real x, tk::real y, tk::real )
985 : : // *****************************************************************************
986 : : //! Evaluate pressure boundary condition
987 : : //! \param[in] x X coordinate where to evaluate the IC / analytic solution
988 : : //! \param[in] y Y coordinate where to evaluate the IC / analytic solution
989 : : //! \return Value for pressure
990 : : // *****************************************************************************
991 : : {
992 : 972 : return std::cos(2.0*x) * std::exp(y);
993 : : }
994 : :
995 : : } // poisson_neumann::
996 : :
997 : :
998 : : namespace poiseuille {
999 : :
1000 : : static std::vector< tk::real >
1001 : 0 : ic( tk::real, tk::real y, tk::real, tk::real )
1002 : : // *****************************************************************************
1003 : : //! Set initial conditions prescribing the Poisuille problem
1004 : : //! \param[in] y Y coordinate where to evaluate the solution
1005 : : //! \return Values of physics variables
1006 : : // *****************************************************************************
1007 : : {
1008 : 0 : auto eps = std::numeric_limits< tk::real >::epsilon();
1009 : 0 : auto nu = g_cfg.get< tag::mat_dyn_viscosity >();
1010 [ - - ][ - - ]: 0 : if (nu < eps) Throw( "Poiseuille flow needs nonzero viscosity" );
[ - - ][ - - ]
1011 : :
1012 : 0 : auto dpdx = -0.12;
1013 : 0 : auto u = -dpdx * y * (1.0 - y) / 2.0 / nu;
1014 : :
1015 : 0 : const auto& solver = g_cfg.get< tag::solver >();
1016 [ - - ]: 0 : if (solver == "chocg") {
1017 [ - - ]: 0 : return { u, 0.0, 0.0 };
1018 : : } else
1019 [ - - ]: 0 : if (solver == "lohcg") {
1020 [ - - ]: 0 : return { 0.0, 0.0, 0.0, 0.0 };
1021 : : }
1022 [ - - ][ - - ]: 0 : else Throw( "Poiseuille IC not setup for this solver" );
[ - - ]
1023 : : }
1024 : :
1025 : : static std::vector< tk::real >
1026 : 0 : sol( tk::real, tk::real y, tk::real, tk::real )
1027 : : // *****************************************************************************
1028 : : //! Set analytic solution of the Poisuille problem
1029 : : //! \param[in] y Y coordinate where to evaluate the solution
1030 : : //! \return Values of physics variables
1031 : : // *****************************************************************************
1032 : : {
1033 : 0 : auto eps = std::numeric_limits< tk::real >::epsilon();
1034 : 0 : auto nu = g_cfg.get< tag::mat_dyn_viscosity >();
1035 [ - - ][ - - ]: 0 : if (nu < eps) Throw( "Poiseuille flow needs nonzero viscosity" );
[ - - ][ - - ]
1036 : :
1037 : 0 : auto dpdx = -0.12;
1038 : 0 : auto u = -dpdx * y * (1.0 - y) / 2.0 / nu;
1039 : :
1040 : 0 : const auto& solver = g_cfg.get< tag::solver >();
1041 [ - - ]: 0 : if (solver == "chocg") {
1042 [ - - ]: 0 : return { u, 0.0, 0.0 };
1043 : : } else
1044 [ - - ]: 0 : if (solver == "lohcg") {
1045 [ - - ]: 0 : return { 0.0, u, 0.0, 0.0 };
1046 : : }
1047 [ - - ][ - - ]: 0 : else Throw( "Poiseuille IC not setup for this solver" );
[ - - ]
1048 : : }
1049 : :
1050 : : } // poiseuille::
1051 : :
1052 : : std::function< std::vector< tk::real >
1053 : : ( tk::real, tk::real, tk::real, tk::real ) >
1054 : 65648 : IC()
1055 : : // *****************************************************************************
1056 : : // Query user config and assign function to set initial conditions
1057 : : //! \return The function to call to set initial conditions
1058 : : // *****************************************************************************
1059 : : {
1060 : 65648 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1061 : :
1062 : : std::function< std::vector< tk::real >
1063 : 65648 : ( tk::real, tk::real, tk::real, tk::real ) > ic;
1064 : :
1065 [ + + ]: 65648 : if (problem == "userdef")
1066 : 6356 : ic = userdef::ic;
1067 [ + + ]: 59292 : else if (problem == "nonlinear_energy_growth")
1068 : 1254 : ic = nonlinear_energy_growth::ic;
1069 [ + + ]: 58038 : else if (problem == "rayleigh_taylor")
1070 : 3624 : ic = rayleigh_taylor::ic;
1071 [ + + ]: 54414 : else if (problem == "sedov")
1072 : 59 : ic = sedov::ic;
1073 [ + + ]: 54355 : else if (problem == "sod")
1074 : 629 : ic = sod::ic;
1075 [ + + ]: 53726 : else if (problem == "taylor_green")
1076 : 4812 : ic = taylor_green::ic;
1077 [ + + ]: 48914 : else if (problem == "vortical_flow")
1078 : 29151 : ic = vortical_flow::ic;
1079 [ + + ]: 19763 : else if (problem == "slot_cyl")
1080 : 19624 : ic = slot_cyl::ic;
1081 [ - + ]: 139 : else if (problem == "sheardiff")
1082 : 0 : ic = sheardiff::ic;
1083 [ + + ]: 139 : else if (problem == "point_src")
1084 : 4 : ic = point_src::ic;
1085 [ + - ]: 135 : else if (problem.find("poisson") != std::string::npos)
1086 : 135 : ic = poisson::ic;
1087 [ - - ]: 0 : else if (problem == "poiseuille")
1088 : 0 : ic = poiseuille::ic;
1089 : : else
1090 [ - - ][ - - ]: 0 : Throw( "problem type ic not hooked up" );
[ - - ]
1091 : :
1092 : 65648 : return ic;
1093 : 0 : }
1094 : :
1095 : : std::function< std::vector< tk::real >
1096 : : ( tk::real, tk::real, tk::real, tk::real ) >
1097 : 47235 : SOL()
1098 : : // *****************************************************************************
1099 : : // Query user config and assign function to query analytic solutions
1100 : : //! \return The function to call to query analytic solutions
1101 : : // *****************************************************************************
1102 : : {
1103 : 47235 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1104 : :
1105 [ + + ]: 71692 : if (problem == "userdef" or
1106 [ + + ]: 46499 : problem == "sod" or
1107 [ + + ][ + + ]: 93734 : problem == "sedov" or
[ + + ]
1108 : 20028 : problem == "point_src")
1109 : 27263 : return {};
1110 [ - + ]: 19972 : else if (problem == "poiseuille")
1111 : 0 : return poiseuille::sol;
1112 : : else
1113 : 19972 : return IC();
1114 : : }
1115 : :
1116 : : void
1117 : 2786 : initialize( const std::array< std::vector< tk::real >, 3 >& coord,
1118 : : tk::Fields& U,
1119 : : tk::real t,
1120 : : const std::vector< std::unordered_set< std::size_t > >& boxnodes )
1121 : : // *****************************************************************************
1122 : : // Set inital conditions
1123 : : //! \param[in] coord Mesh node coordinates
1124 : : //! \param[in,out] U Array of unknowns
1125 : : //! \param[in] t Physical time
1126 : : //! \param[in] boxnodes Nodes at which box user ICs are set (for each box IC)
1127 : : // *****************************************************************************
1128 : : {
1129 [ - + ][ - - ]: 2786 : Assert( coord[0].size() == U.nunk(), "Size mismatch" );
[ - - ][ - - ]
1130 : :
1131 [ + - ]: 2786 : auto ic = IC();
1132 : 2786 : const auto& x = coord[0];
1133 : 2786 : const auto& y = coord[1];
1134 : 2786 : const auto& z = coord[2];
1135 : :
1136 : : // Set initial conditions dependeing on problem configured
1137 [ + + ]: 385551 : for (std::size_t i=0; i<x.size(); ++i) {
1138 : :
1139 : : // Set background ICs
1140 [ + - ]: 382765 : auto s = ic( x[i], y[i], z[i], t );
1141 [ - + ][ - - ]: 382765 : Assert( s.size() == U.nprop(), "Size mismatch" );
[ - - ][ - - ]
1142 : :
1143 : : // Initialize user-defined ICs in boxes
1144 [ + - ]: 382765 : box( i, s, boxnodes );
1145 : :
1146 : : // Set values for ICs
1147 [ + - ][ + + ]: 2324142 : for (std::size_t c=0; c<s.size(); ++c) U(i,c) = s[c];
1148 : :
1149 : 382765 : }
1150 : 2786 : }
1151 : :
1152 : : std::function< tk::real( tk::real, tk::real, tk::real ) >
1153 : 6198 : PRESSURE_RHS()
1154 : : // *****************************************************************************
1155 : : // Query user config and assign function to set pressure rhs
1156 : : //! \return The function to call to set pressure rhs
1157 : : // *****************************************************************************
1158 : : {
1159 : 6198 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1160 : :
1161 : 6198 : std::function< tk::real( tk::real, tk::real, tk::real ) > pr;
1162 : :
1163 [ + + ]: 6198 : if (problem == "poisson_const")
1164 : 22 : pr = poisson_const::pr;
1165 [ - + ]: 6176 : else if (problem == "poisson_harmonic")
1166 : 0 : pr = poisson_harmonic::pr;
1167 [ + + ]: 6176 : else if (problem == "poisson_sine")
1168 : 2 : pr = poisson_sine::pr;
1169 [ + + ]: 6174 : else if (problem == "poisson_sine3")
1170 : 6 : pr = poisson_sine3::pr;
1171 [ + + ]: 6168 : else if (problem == "poisson_neumann")
1172 : 2 : pr = poisson_neumann::pr;
1173 : :
1174 : 6198 : return pr;
1175 : : }
1176 : :
1177 : : std::function< tk::real( tk::real, tk::real, tk::real ) >
1178 : 3412 : PRESSURE_IC()
1179 : : // *****************************************************************************
1180 : : // Query user config and assign function to set pressure initial conditions
1181 : : //! \return The function to call to set pressure initial conditions
1182 : : // *****************************************************************************
1183 : : {
1184 : 3412 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1185 : :
1186 : 3412 : std::function< tk::real( tk::real, tk::real, tk::real ) > ic;
1187 : :
1188 [ + + ]: 3650 : if ( problem == "userdef" or
1189 [ + - ]: 380 : problem == "slot_cyl" or
1190 [ + + ][ - + ]: 3792 : problem == "sheardiff" or
[ + + ]
1191 : 142 : problem == "poiseuille" )
1192 : 3270 : ic = userdef::pic;
1193 [ + + ]: 142 : else if (problem == "poisson_const")
1194 : 94 : ic = poisson_const::ic;
1195 [ - + ]: 48 : else if (problem == "poisson_harmonic")
1196 : 0 : ic = poisson_harmonic::ic;
1197 [ + + ]: 48 : else if (problem == "poisson_sine")
1198 : 10 : ic = poisson_sine::ic;
1199 [ + + ]: 38 : else if (problem == "poisson_sine3")
1200 : 28 : ic = poisson_sine3::ic;
1201 [ + - ]: 10 : else if (problem == "poisson_neumann")
1202 : 10 : ic = poisson_neumann::ic;
1203 : : else
1204 [ - - ][ - - ]: 0 : Throw( "problem type not hooked up" );
[ - - ]
1205 : :
1206 : 3412 : return ic;
1207 : 0 : }
1208 : :
1209 : : std::function< tk::real( tk::real, tk::real, tk::real ) >
1210 : 5612 : PRESSURE_SOL()
1211 : : // *****************************************************************************
1212 : : // Query user config and assign function to query analytic pressure solutions
1213 : : //! \return The function to call to query analytic solutions
1214 : : // *****************************************************************************
1215 : : {
1216 : 5612 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1217 : :
1218 [ + + ]: 6906 : if ( problem == "userdef" or
1219 [ + - ]: 1404 : problem == "slot_cyl" or
1220 [ + + ][ - + ]: 7016 : problem == "poiseuille" or
[ + + ]
1221 : 110 : problem == "sheardiff" )
1222 : : {
1223 : 5502 : return {};
1224 : : }
1225 : : else {
1226 : 110 : return PRESSURE_IC();
1227 : : }
1228 : : }
1229 : :
1230 : : std::function< std::array< tk::real, 3 >( tk::real, tk::real, tk::real ) >
1231 : 6198 : PRESSURE_GRAD()
1232 : : // *****************************************************************************
1233 : : // Assign function to query pressure gradient at a point
1234 : : //! \return The function to call to query the pressure gradient
1235 : : // *****************************************************************************
1236 : : {
1237 : 6198 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1238 : :
1239 [ + + ]: 6198 : if (problem == "poisson_neumann")
1240 : 2 : return poisson_neumann::pg;
1241 : :
1242 : 6196 : return {};
1243 : : }
1244 : :
1245 : : tk::real
1246 : 0 : initialize( tk::real x, tk::real y, tk::real z )
1247 : : // *****************************************************************************
1248 : : // Evaluate initial condition for pressure
1249 : : //! \param[in] x X coordinate where to evaluate the pressure initial condition
1250 : : //! \param[in] y Y coordinate where to evaluate the pressure initial condition
1251 : : //! \param[in] z Z coordinate where to evaluate the pressure initial condition
1252 : : //! \return Pressure initial condition
1253 : : // *****************************************************************************
1254 : : {
1255 : 0 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1256 : :
1257 : 0 : std::function< tk::real( tk::real, tk::real, tk::real ) > ic;
1258 : :
1259 [ - - ]: 0 : if (problem == "poisson_const")
1260 : 0 : ic = poisson_const::ic;
1261 [ - - ]: 0 : else if (problem == "poisson_harmonic")
1262 : 0 : ic = poisson_harmonic::ic;
1263 [ - - ]: 0 : else if (problem == "poisson_sine")
1264 : 0 : ic = poisson_sine::ic;
1265 [ - - ]: 0 : else if (problem == "poisson_sine3")
1266 : 0 : ic = poisson_sine3::ic;
1267 : : else
1268 [ - - ][ - - ]: 0 : Throw( "problem type not hooked up" );
[ - - ]
1269 : :
1270 [ - - ]: 0 : return ic( x, y, z );
1271 : 0 : }
1272 : :
1273 : : std::function< std::vector< tk::real >
1274 : : ( tk::real, tk::real, tk::real, tk::real ) >
1275 : 78156 : SRC()
1276 : : // *****************************************************************************
1277 : : // Query user config and assign function to add a source term
1278 : : //! \return The function to call to evaluate a problem-sepcific source term
1279 : : // *****************************************************************************
1280 : : {
1281 : 78156 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1282 : :
1283 : : std::function<
1284 : 78156 : std::vector< tk::real >( tk::real, tk::real, tk::real, tk::real ) > src;
1285 : :
1286 [ + + ]: 78156 : if (problem == "nonlinear_energy_growth")
1287 : 676 : src = nonlinear_energy_growth::src;
1288 [ + + ]: 77480 : else if (problem == "rayleigh_taylor")
1289 : 2000 : src = rayleigh_taylor::src;
1290 [ + + ]: 75480 : else if (problem == "taylor_green")
1291 : 3536 : src = taylor_green::src;
1292 [ + + ]: 71944 : else if (problem == "vortical_flow")
1293 : 18547 : src = vortical_flow::src;
1294 [ + + ]: 53397 : else if (problem == "slot_cyl")
1295 : 11537 : src = slot_cyl::src;
1296 [ - + ]: 41860 : else if (problem == "sheardiff")
1297 : 0 : src = sheardiff::src;
1298 : :
1299 : 78156 : return src;
1300 : : }
1301 : :
1302 : : std::function< void( const std::array< std::vector< tk::real >, 3 >&,
1303 : : tk::real,
1304 : : tk::Fields& ) >
1305 : 78156 : PHYS_SRC()
1306 : : // *****************************************************************************
1307 : : // Query user config and assign function to apply source to numerical solution
1308 : : //! \return The function to call to evaluate a problem-sepcific source term
1309 : : // *****************************************************************************
1310 : : {
1311 : 78156 : const auto& problem = inciter::g_cfg.get< tag::problem >();
1312 : :
1313 : : std::function< void( const std::array< std::vector< tk::real >, 3 >&,
1314 : : tk::real,
1315 : 78156 : tk::Fields& ) > src;
1316 : :
1317 [ + + ]: 78156 : if (problem == "point_src") {
1318 : 600 : src = point_src::src;
1319 : : }
1320 : :
1321 : 78156 : return src;
1322 : : }
1323 : :
1324 : : } // problems::
|
