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 : 1724801 : ic( tk::real, tk::real, tk::real, tk::real, std::size_t meshid )
32 : : // *****************************************************************************
33 : : //! Set homogeneous initial conditions for a generic user-defined problem
34 : : //! \param[in] meshid Mesh id to use
35 : : //! \return Values of conserved variables
36 : : // *****************************************************************************
37 : : {
38 : : bool multi = g_cfg.get< tag::input >().size() > 1;
39 : : const auto& ncomp = g_cfg.get< tag::problem_ncomp >();
40 : :
41 : : // pressure-based solvers
42 : :
43 : : const auto& solver = g_cfg.get< tag::solver >();
44 [ + + ]: 1724801 : if (solver == "chocg") {
45 [ + - ]: 662459 : std::vector< tk::real > u( ncomp, 0.0 );
46 : : const auto& ic_velocity = g_cfg.get< tag::ic, tag::velocity >();
47 : : auto large = std::numeric_limits< double >::max();
48 [ + - ]: 662459 : if (std::abs(ic_velocity[0] - large) > 1.0e-12) u[0] = ic_velocity[0];
49 [ + - ]: 662459 : if (std::abs(ic_velocity[1] - large) > 1.0e-12) u[1] = ic_velocity[1];
50 [ + - ]: 662459 : if (std::abs(ic_velocity[2] - large) > 1.0e-12) u[2] = ic_velocity[2];
51 : : return u;
52 : : }
53 [ + + ]: 1062342 : else if (solver == "lohcg") {
54 : 921902 : std::vector< tk::real > u( ncomp, 0.0 );
55 : :
56 : : const auto& ic_velocity =
57 : : multi ? g_cfg.get< tag::ic_ >()[ meshid ].get< tag::velocity >()
58 [ - + ]: 921902 : : g_cfg.get< tag::ic, tag::velocity >();
59 : :
60 : : auto large = std::numeric_limits< double >::max();
61 [ + - ]: 921902 : if (std::abs(ic_velocity[0] - large) > 1.0e-12) u[1] = ic_velocity[0];
62 [ + - ]: 921902 : if (std::abs(ic_velocity[1] - large) > 1.0e-12) u[2] = ic_velocity[1];
63 [ + - ]: 921902 : if (std::abs(ic_velocity[2] - large) > 1.0e-12) u[3] = ic_velocity[2];
64 : : return u;
65 : : }
66 : :
67 : : // density-based solvers
68 : : const auto& tic = g_cfg.get< tag::ic >();
69 [ - + ]: 140440 : auto ic_density = tic.get< tag::density >();
70 : : const auto& ic_velocity = tic.get< tag::velocity >();
71 [ - + ][ - - ]: 140440 : ErrChk( ic_velocity.size() == 3, "ic_velocity must have 3 components" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
72 : :
73 [ + - ]: 140440 : std::vector< tk::real > u( ncomp, 0.0 );
74 : :
75 [ + - ]: 140440 : u[0] = ic_density;
76 : 140440 : u[1] = u[0] * ic_velocity[0];
77 : 140440 : u[2] = u[0] * ic_velocity[1];
78 : 140440 : u[3] = u[0] * ic_velocity[2];
79 : :
80 : 140440 : auto ic_pressure = tic.get< tag::pressure >();
81 : 140440 : auto ic_energy = tic.get< tag::energy >();
82 : 140440 : auto ic_temperature = tic.get< tag::temperature >();
83 : :
84 : : auto largereal = std::numeric_limits< double >::max();
85 : :
86 [ + - ]: 140440 : if (std::abs(ic_pressure - largereal) > 1.0e-12) {
87 : :
88 : 140440 : u[4] = eos::totalenergy( u[0], u[1]/u[0], u[2]/u[0], u[3]/u[0],
89 : : ic_pressure );
90 : :
91 [ - - ]: 0 : } else if (std::abs(ic_energy - largereal) > 1.0e-12) {
92 : :
93 : 0 : u[4] = u[0] * ic_energy;
94 : :
95 [ - - ]: 0 : } else if (std::abs(ic_temperature - largereal) > 1.0e-12) {
96 : :
97 : 0 : auto cv = g_cfg.get< tag::mat_spec_heat_const_vol >();
98 [ - - ]: 0 : if (std::abs(cv - largereal) > 1.0e-12) {
99 : 0 : u[4] = u[0] * ic_temperature * cv;
100 : : }
101 : :
102 : : } else {
103 : :
104 [ - - ][ - - ]: 0 : Throw( "IC background energy cannot be computed. Must specify "
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
105 : : "one of background pressure, energy, or velocity." );
106 : :
107 : : }
108 : :
109 : : return u;
110 : : }
111 : :
112 : : static tk::real
113 : 16746 : pic( tk::real, tk::real, tk::real, std::size_t )
114 : : // *****************************************************************************
115 : : //! Set homogeneous initial conditions for a generic user-defined problem
116 : : //! \return Value of pressure
117 : : // *****************************************************************************
118 : : {
119 : 16746 : return 0.0;
120 : : }
121 : :
122 : : } // userdef::
123 : :
124 : : namespace nonlinear_energy_growth {
125 : :
126 : : static std::vector< tk::real >
127 : 73908 : ic( tk::real x, tk::real y, tk::real z, tk::real t, std::size_t )
128 : : // *****************************************************************************
129 : : //! Set initial conditions prescribing nonlinear energy growth
130 : : //! \param[in] x X coordinate where to evaluate the solution
131 : : //! \param[in] y Y coordinate where to evaluate the solution
132 : : //! \param[in] z Z coordinate where to evaluate the solution
133 : : //! \param[in] t Time where to evaluate the solution
134 : : //! \return Values of conserved variables
135 : : // *****************************************************************************
136 : : {
137 : : using std::cos;
138 : :
139 : : // manufactured solution parameters
140 : 73908 : auto ce = g_cfg.get< tag::problem_ce >();
141 : 73908 : auto r0 = g_cfg.get< tag::problem_r0 >();
142 : 73908 : auto a = g_cfg.get< tag::problem_alpha >();
143 : 73908 : auto k = g_cfg.get< tag::problem_kappa >();
144 : : const auto& b = g_cfg.get< tag::problem_beta >();
145 : :
146 : : auto ec = [ ce, t ]( tk::real kappa, tk::real h, tk::real p ) {
147 : 73908 : return std::pow( -3.0*(ce + kappa*h*h*t), p );
148 : : };
149 : :
150 : 147816 : auto hx = [ x, y, z, b ]() {
151 : 73908 : return cos(b[0]*M_PI*x) * cos(b[1]*M_PI*y) * cos(b[2]*M_PI*z);
152 : 73908 : };
153 : :
154 : : // density
155 : 73908 : auto r = r0 + std::exp(-a*t) * (1.0 - x*x - y*y - z*z);
156 : : // energy
157 : 73908 : auto re = r * ec(k,hx(),-1.0/3.0);
158 : :
159 [ + - ][ + - ]: 147816 : return { r, 0.0, 0.0, 0.0, re };
[ - - ]
160 : : }
161 : :
162 : : static std::vector< tk::real >
163 : 178308 : src( tk::real x, tk::real y, tk::real z, tk::real t, std::size_t )
164 : : // *****************************************************************************
165 : : //! Compute and return source term for nonlinear energy growth
166 : : //! \param[in] x X coordinate where to evaluate the source
167 : : //! \param[in] y Y coordinate where to evaluate the source
168 : : //! \param[in] z Z coordinate where to evaluate the source
169 : : //! \param[in] t Time where to evaluate the source
170 : : //! \return Source for flow variables + transported scalars
171 : : // *****************************************************************************
172 : : {
173 : : using std::sin; using std::cos; using std::pow;
174 : :
175 : : // manufactured solution parameters
176 : 178308 : auto a = g_cfg.get< tag::problem_alpha >();
177 : : const auto& b = g_cfg.get< tag::problem_beta >();
178 : 178308 : auto ce = g_cfg.get< tag::problem_ce >();
179 : 178308 : auto kappa = g_cfg.get< tag::problem_kappa >();
180 : 178308 : auto r0 = g_cfg.get< tag::problem_r0 >();
181 : : // ratio of specific heats
182 : 178308 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
183 : : // spatial component of density field
184 : 178308 : auto gx = 1.0 - x*x - y*y - z*z;
185 : : // derivative of spatial component of density field
186 : 178308 : std::array< tk::real, 3 > dg{ -2.0*x, -2.0*y, -2.0*z };
187 : : // spatial component of energy field
188 : 178308 : auto h = cos(b[0]*M_PI*x) * cos(b[1]*M_PI*y) * cos(b[2]*M_PI*z);
189 : : // derivative of spatial component of energy field
190 : : std::array< tk::real, 3 >
191 : 178308 : dh{ -b[0]*M_PI*sin(b[0]*M_PI*x)*cos(b[1]*M_PI*y)*cos(b[2]*M_PI*z),
192 : 178308 : -b[1]*M_PI*cos(b[0]*M_PI*x)*sin(b[1]*M_PI*y)*cos(b[2]*M_PI*z),
193 : 178308 : -b[2]*M_PI*cos(b[0]*M_PI*x)*cos(b[1]*M_PI*y)*sin(b[2]*M_PI*z) };
194 : : // temporal function f and its derivative
195 : 178308 : auto ft = std::exp(-a*t);
196 : 178308 : auto dfdt = -a*ft;
197 : : // density and its derivatives
198 : 178308 : auto rho = r0 + ft*gx;
199 : 178308 : std::array< tk::real, 3 > drdx{ ft*dg[0], ft*dg[1], ft*dg[2] };
200 : 178308 : auto drdt = gx*dfdt;
201 : : // internal energy and its derivatives
202 : 178308 : auto ie = pow( -3.0*(ce + kappa*h*h*t), -1.0/3.0 );
203 : 178308 : std::array< tk::real, 3 > dedx{ 2.0 * pow(ie,4.0) * kappa * h * dh[0] * t,
204 : 178308 : 2.0 * pow(ie,4.0) * kappa * h * dh[1] * t,
205 : 178308 : 2.0 * pow(ie,4.0) * kappa * h * dh[2] * t };
206 : 178308 : const auto dedt = kappa * h * h * pow(ie,4.0);
207 : :
208 : 178308 : std::vector< tk::real > s( 5, 0.0 );
209 : : // density source
210 : 178308 : s[0] = drdt;
211 : : // momentum source
212 : 178308 : s[1] = (g-1.0)*(rho*dedx[0] + ie*drdx[0]);
213 : 178308 : s[2] = (g-1.0)*(rho*dedx[1] + ie*drdx[1]);
214 : 178308 : s[3] = (g-1.0)*(rho*dedx[2] + ie*drdx[2]);
215 : : // energy source
216 : 178308 : s[4] = rho*dedt + ie*drdt;
217 : :
218 : 178308 : return s;
219 : : }
220 : :
221 : : } // nonlinear_energy_growth::
222 : :
223 : : namespace rayleigh_taylor {
224 : :
225 : : static std::vector< tk::real >
226 : 1058942 : ic( tk::real x, tk::real y, tk::real z, tk::real t, std::size_t )
227 : : // *****************************************************************************
228 : : //! Set initial conditions prescribing a Rayleigh-Taylor flow
229 : : //! \param[in] x X coordinate where to evaluate the solution
230 : : //! \param[in] y Y coordinate where to evaluate the solution
231 : : //! \param[in] z Z coordinate where to evaluate the solution
232 : : //! \param[in] t Time where to evaluate the solution
233 : : //! \return Values of conserved variables
234 : : // *****************************************************************************
235 : : {
236 : : using std::sin; using std::cos;
237 : :
238 : : // manufactured solution parameters
239 : 1058942 : auto a = g_cfg.get< tag::problem_alpha >();
240 : : const auto& b = g_cfg.get< tag::problem_beta >();
241 : 1058942 : auto p0 = g_cfg.get< tag::problem_p0 >();
242 : 1058942 : auto r0 = g_cfg.get< tag::problem_r0 >();
243 : 1058942 : auto k = g_cfg.get< tag::problem_kappa >();
244 : :
245 : : // spatial component of density and pressure fields
246 : 1058942 : tk::real gx = b[0]*x*x + b[1]*y*y + b[2]*z*z;
247 : : // density
248 : 1058942 : tk::real r = r0 - gx;
249 : : // velocity
250 : 1058942 : tk::real ft = cos(k*M_PI*t);
251 : 1058942 : tk::real u = ft * z * sin(M_PI*x);
252 : 1058942 : tk::real v = ft * z * cos(M_PI*y);
253 : 1058942 : tk::real w = ft * ( -0.5*M_PI*z*z*(cos(M_PI*x) - sin(M_PI*y)) );
254 : : // total specific energy
255 : 1058942 : tk::real rE = eos::totalenergy( r, u, v, w, p0 + a*gx );
256 : :
257 : 1058942 : return { r, r*u, r*v, r*w, rE };
258 : : }
259 : :
260 : : static std::vector< tk::real >
261 : 768460 : src( tk::real x, tk::real y, tk::real z, tk::real t, std::size_t meshid )
262 : : // *****************************************************************************
263 : : //! Compute and return source term for a Rayleigh-Taylor flow
264 : : //! \param[in] x X coordinate where to evaluate the source
265 : : //! \param[in] y Y coordinate where to evaluate the source
266 : : //! \param[in] z Z coordinate where to evaluate the source
267 : : //! \param[in] t Time where to evaluate the source
268 : : //! \param[in] meshid Mesh id to use
269 : : //! \return Source for flow variables + transported scalars
270 : : // *****************************************************************************
271 : : {
272 : : using std::sin; using std::cos;
273 : :
274 : : // manufactured solution parameters
275 : 768460 : auto a = g_cfg.get< tag::problem_alpha >();
276 : : const auto& b = g_cfg.get< tag::problem_beta >();
277 : 768460 : auto k = g_cfg.get< tag::problem_kappa >();
278 : 768460 : auto p0 = g_cfg.get< tag::problem_p0 >();
279 : 768460 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
280 : :
281 : : // evaluate solution at x,y,z,t
282 : 768460 : auto U = ic( x, y, z, t, meshid );
283 : :
284 : : // density, velocity, energy, pressure
285 : 768460 : auto rho = U[0];
286 : 768460 : auto u = U[1]/U[0];
287 : 768460 : auto v = U[2]/U[0];
288 : 768460 : auto w = U[3]/U[0];
289 [ + - ]: 768460 : auto E = U[4]/U[0];
290 : 768460 : auto p = p0 + a*(b[0]*x*x + b[1]*y*y + b[2]*z*z);
291 : :
292 : : // spatial gradients
293 : 768460 : std::array< tk::real, 3 > drdx{{ -2.0*b[0]*x, -2.0*b[1]*y, -2.0*b[2]*z }};
294 : 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 }};
295 : 768460 : tk::real ft = cos(k*M_PI*t);
296 : 768460 : std::array< tk::real, 3 > dudx{{ ft*M_PI*z*cos(M_PI*x),
297 : : 0.0,
298 : 768460 : ft*sin(M_PI*x) }};
299 : : std::array< tk::real, 3 > dvdx{{ 0.0,
300 : 768460 : -ft*M_PI*z*sin(M_PI*y),
301 : 768460 : ft*cos(M_PI*y) }};
302 : 768460 : std::array< tk::real, 3 > dwdx{{ ft*M_PI*0.5*M_PI*z*z*sin(M_PI*x),
303 : 768460 : ft*M_PI*0.5*M_PI*z*z*cos(M_PI*y),
304 : 768460 : -ft*M_PI*z*(cos(M_PI*x) - sin(M_PI*y)) }};
305 : : std::array< tk::real, 3 > dedx{{
306 : 768460 : dpdx[0]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[0]
307 : 768460 : + u*dudx[0] + v*dvdx[0] + w*dwdx[0],
308 : 768460 : dpdx[1]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[1]
309 : 768460 : + u*dudx[1] + v*dvdx[1] + w*dwdx[1],
310 : 768460 : dpdx[2]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[2]
311 : 768460 : + u*dudx[2] + v*dvdx[2] + w*dwdx[2] }};
312 : :
313 : : // time derivatives
314 : 768460 : auto dudt = -k*M_PI*sin(k*M_PI*t)*z*sin(M_PI*x);
315 : 768460 : auto dvdt = -k*M_PI*sin(k*M_PI*t)*z*cos(M_PI*y);
316 : 768460 : auto dwdt = k*M_PI*sin(k*M_PI*t)/2*M_PI*z*z*(cos(M_PI*x) - sin(M_PI*y));
317 : 768460 : auto dedt = u*dudt + v*dvdt + w*dwdt;
318 : :
319 [ + - ]: 768460 : std::vector< tk::real > s( 5, 0.0 );
320 : : // density source
321 : 768460 : s[0] = u*drdx[0] + v*drdx[1] + w*drdx[2];
322 : : // momentum source
323 : 768460 : s[1] = rho*dudt+u*s[0]+dpdx[0] + U[1]*dudx[0]+U[2]*dudx[1]+U[3]*dudx[2];
324 : 768460 : s[2] = rho*dvdt+v*s[0]+dpdx[1] + U[1]*dvdx[0]+U[2]*dvdx[1]+U[3]*dvdx[2];
325 : 768460 : s[3] = rho*dwdt+w*s[0]+dpdx[2] + U[1]*dwdx[0]+U[2]*dwdx[1]+U[3]*dwdx[2];
326 : : // energy source
327 : 768460 : s[4] = rho*dedt + E*s[0] + U[1]*dedx[0]+U[2]*dedx[1]+U[3]*dedx[2]
328 : 768460 : + u*dpdx[0]+v*dpdx[1]+w*dpdx[2];
329 : :
330 : 768460 : return s;
331 : : }
332 : :
333 : : } // rayleigh_taylor::
334 : :
335 : : namespace sedov {
336 : : static std::vector< tk::real >
337 [ + + ]: 66564 : ic( tk::real x, tk::real y, tk::real z, tk::real, std::size_t )
338 : : // *****************************************************************************
339 : : //! Set initial conditions prescribing the Sedov blast wave
340 : : //! \param[in] x X coordinate where to evaluate the solution
341 : : //! \param[in] y Y coordinate where to evaluate the solution
342 : : //! \param[in] z Z coordinate where to evaluate the solution
343 : : //! \return Values of conserved variables
344 : : // *****************************************************************************
345 : : {
346 : : using std::abs;
347 : :
348 : : // pressure
349 : : auto eps = std::numeric_limits< tk::real >::epsilon();
350 : : tk::real p;
351 [ + + ][ + + ]: 66564 : if (abs(x) < eps && abs(y) < eps && abs(z) < eps) {
[ + + ]
352 : 5 : p = g_cfg.get< tag::problem_p0 >();
353 : : } else {
354 : : p = 0.67e-4;
355 : : }
356 : :
357 : : // density
358 : : tk::real r = 1.0;
359 : : // velocity
360 : : tk::real u = 0.0;
361 : : tk::real v = 0.0;
362 : : tk::real w = 0.0;
363 : : // total specific energy
364 : : tk::real rE = eos::totalenergy( r, u, v, w, p );
365 : :
366 : 66564 : return { r, r*u, r*v, r*w, rE };
367 : :
368 : : }
369 : : } // sedov::
370 : :
371 : : namespace sod {
372 : : static std::vector< tk::real >
373 : 13179 : ic( tk::real x, tk::real, tk::real, tk::real, std::size_t )
374 : : // *****************************************************************************
375 : : //! Set initial conditions prescribing the Sod shocktube
376 : : //! \param[in] x X coordinate where to evaluate the solution
377 : : //! \return Values of conserved variables
378 : : // *****************************************************************************
379 : : {
380 : : tk::real r, p, u, v, w, rE;
381 : :
382 [ + + ]: 13179 : if (x<0.5) {
383 : : // density
384 : : r = 1.0;
385 : : // pressure
386 : : p = 1.0;
387 : : }
388 : : else {
389 : : // density
390 : : r = 0.125;
391 : : // pressure
392 : : p = 0.1;
393 : : }
394 : :
395 : : // velocity
396 : : u = 0.0;
397 : : v = 0.0;
398 : : w = 0.0;
399 : :
400 : : // total specific energy
401 : : rE = eos::totalenergy( r, u, v, w, p );
402 : :
403 : 13179 : return { r, r*u, r*v, r*w, rE };
404 : : }
405 : : } // sod::
406 : :
407 : : namespace taylor_green {
408 : :
409 : : static std::vector< tk::real >
410 : 280356 : ic( tk::real x, tk::real y, tk::real, tk::real, std::size_t )
411 : : // *****************************************************************************
412 : : //! Set initial conditions prescribing the Taylor-Green vortex
413 : : //! \param[in] x X coordinate where to evaluate the solution
414 : : //! \param[in] y Y coordinate where to evaluate the solution
415 : : //! \return Values of conserved variables
416 : : // *****************************************************************************
417 : : {
418 : : // density
419 : : tk::real r = 1.0;
420 : : // pressure
421 : 280356 : tk::real p = 10.0 + r/4.0*(cos(2.0*M_PI*x) + cos(2.0*M_PI*y));
422 : : // velocity
423 : 280356 : tk::real u = sin(M_PI*x) * cos(M_PI*y);
424 : 280356 : tk::real v = -cos(M_PI*x) * sin(M_PI*y);
425 : : tk::real w = 0.0;
426 : : // total specific energy
427 : : auto rE = eos::totalenergy( r, u, v, w, p );
428 : :
429 : 280356 : return { r, r*u, r*v, r*w, rE };
430 : : }
431 : :
432 : : static std::vector< tk::real >
433 : 932688 : src( tk::real x, tk::real y, tk::real, tk::real, std::size_t )
434 : : // *****************************************************************************
435 : : //! Compute and return source term for a the Taylor-Green vortex
436 : : //! \param[in] x X coordinate where to evaluate the source
437 : : //! \param[in] y Y coordinate where to evaluate the source
438 : : //! \return Source for flow variables + transported scalars
439 : : // *****************************************************************************
440 : : {
441 : : using std::cos;
442 : :
443 : 932688 : std::vector< tk::real > s( 5, 0.0 );
444 : 932688 : s[4] = 3.0*M_PI/8.0*( cos(3.0*M_PI*x)*cos(M_PI*y)
445 : 932688 : - cos(3.0*M_PI*y)*cos(M_PI*x) );
446 : :
447 : 932688 : return s;
448 : : }
449 : :
450 : : } // taylor_green::
451 : :
452 : : namespace vortical_flow {
453 : :
454 : : static std::vector< tk::real >
455 : 2350800 : ic( tk::real x, tk::real y, tk::real z, tk::real, std::size_t )
456 : : // *****************************************************************************
457 : : //! Set initial conditions prescribing vortical flow
458 : : //! \param[in] x X coordinate where to evaluate the solution
459 : : //! \param[in] y Y coordinate where to evaluate the solution
460 : : //! \param[in] z Z coordinate where to evaluate the solution
461 : : //! \return Values of conserved variables
462 : : // *****************************************************************************
463 : : {
464 : : // manufactured solution parameters
465 : 2350800 : tk::real a = g_cfg.get< tag::problem_alpha >();
466 : 2350800 : tk::real k = g_cfg.get< tag::problem_kappa >();
467 : 2350800 : tk::real p0 = g_cfg.get< tag::problem_p0 >();
468 : : // ratio of specific heats
469 : 2350800 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
470 : : // velocity
471 : 2350800 : tk::real ru = a*x - k*y;
472 : 2350800 : tk::real rv = k*x + a*y;
473 : 2350800 : tk::real rw = -2.0*a*z;
474 : : // total specific energy
475 : 2350800 : tk::real rE = (ru*ru + rv*rv + rw*rw)/2.0 + (p0 - 2.0*a*a*z*z) / (g - 1.0);
476 : :
477 : 2350800 : return { 1.0, ru, rv, rw, rE };
478 : : }
479 : :
480 : : static std::vector< tk::real >
481 : 1474909 : src( tk::real x, tk::real y, tk::real z, tk::real, std::size_t meshid )
482 : : // *****************************************************************************
483 : : //! Compute and return source term for vortical flow
484 : : //! \param[in] x X coordinate where to evaluate the source
485 : : //! \param[in] y Y coordinate where to evaluate the source
486 : : //! \param[in] z Z coordinate where to evaluate the source
487 : : //! \param[in] meshid Mesh id to use
488 : : //! \return Source for flow variables + transported scalars
489 : : // *****************************************************************************
490 : : {
491 : : // manufactured solution parameters
492 : 1474909 : auto a = g_cfg.get< tag::problem_alpha >();
493 : 1474909 : auto k = g_cfg.get< tag::problem_kappa >();
494 : : // ratio of specific heats
495 : 1474909 : auto g = g_cfg.get< tag::mat_spec_heat_ratio >();
496 : : // evaluate solution at x,y,z
497 : 1474909 : auto u = ic( x, y, z, 0.0, meshid );
498 : :
499 [ + - ][ - - ]: 1474909 : std::vector< tk::real > s( 5, 0.0 );
500 : : // momentum source
501 : 1474909 : s[1] = a*u[1]/u[0] - k*u[2]/u[0];
502 : 1474909 : s[2] = k*u[1]/u[0] + a*u[2]/u[0];
503 : : // energy source
504 : 1474909 : s[4] = (s[1]*u[1] + s[2]*u[2])/u[0] + 8.0*a*a*a*z*z/(g-1.0);
505 : :
506 : 1474909 : return s;
507 : : }
508 : :
509 : : } // vortical_flow::
510 : :
511 : : namespace slot_cyl {
512 : :
513 : : static std::vector< tk::real >
514 : 6795244 : ic( tk::real x, tk::real y, tk::real, tk::real t, std::size_t )
515 : : // *****************************************************************************
516 : : //! Set initial conditions prescribing slotted cylinder, cone, Gauss hump
517 : : //! \param[in] x X coordinate where to evaluate the solution
518 : : //! \param[in] y Y coordinate where to evaluate the solution
519 : : //! \param[in] t Time where to evaluate the solution
520 : : //! \return Values of conserved variables
521 : : // *****************************************************************************
522 : : {
523 : : using std::sin; using std::cos; using std::sqrt;
524 : :
525 : : // manufactured solution parameters
526 : : tk::real p0 = 1.0;
527 : :
528 : : // configure number of scalar components
529 : : std::size_t ncomp = 6;
530 : : const auto& solver = g_cfg.get< tag::solver >();
531 [ + + ]: 6795244 : if (solver == "chocg") {
532 : : ncomp = 4;
533 : : }
534 [ + + ]: 5551919 : else if (solver == "lohcg") {
535 : : ncomp = 5;
536 : : }
537 : :
538 : 6795244 : std::vector< tk::real > u( ncomp, 0.0 );
539 : :
540 : : // prescribed velocity: rotate in x-y plane
541 : : std::size_t sc = 5;
542 [ + + ]: 6795244 : if (solver == "chocg") {
543 : 1243325 : u[0] = 0.5 - y;
544 : 1243325 : u[1] = x - 0.5;
545 : 1243325 : u[2] = 0.0;
546 : : sc = 3;
547 : : }
548 [ + + ]: 5551919 : else if (solver == "lohcg") {
549 : 1251220 : u[0] = 0.0;
550 : 1251220 : u[1] = 0.5 - y;
551 : 1251220 : u[2] = x - 0.5;
552 : 1251220 : u[3] = 0.0;
553 : : sc = 4;
554 : : }
555 : : else {
556 : 4300699 : u[0] = 1.0;
557 : 4300699 : u[1] = u[0] * (0.5 - y);
558 : 4300699 : u[2] = u[0] * (x - 0.5);
559 : 4300699 : u[3] = 0.0;
560 : 4300699 : u[4] = eos::totalenergy( u[0], u[1]/u[0], u[2]/u[0], u[3]/u[0], p0 );
561 : : }
562 : :
563 : : const tk::real R0 = 0.15;
564 : :
565 : : // center of the cone
566 : : tk::real x0 = 0.5;
567 : : tk::real y0 = 0.25;
568 : : tk::real r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
569 : 6795244 : tk::real kx = 0.5 + r*sin( t );
570 : 6795244 : tk::real ky = 0.5 - r*cos( t );
571 : :
572 : : // center of the hump
573 : : x0 = 0.25;
574 : : y0 = 0.5;
575 : : r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
576 : 6795244 : tk::real hx = 0.5 + r*sin( t-M_PI/2.0 ),
577 : 6795244 : hy = 0.5 - r*cos( t-M_PI/2.0 );
578 : :
579 : : // center of the slotted cylinder
580 : : x0 = 0.5;
581 : : y0 = 0.75;
582 : : r = sqrt((x0-0.5)*(x0-0.5) + (y0-0.5)*(y0-0.5));
583 : 6795244 : tk::real cx = 0.5 + r*sin( t+M_PI ),
584 : 6795244 : cy = 0.5 - r*cos( t+M_PI );
585 : :
586 : : // end points of the cylinder slot
587 : 6795244 : tk::real i1x = 0.525, i1y = cy - r*cos( std::asin(0.025/r) ),
588 : : i2x = 0.525, i2y = 0.8,
589 : : i3x = 0.475, i3y = 0.8;
590 : :
591 : : // rotate end points of cylinder slot
592 : 6795244 : tk::real ri1x = 0.5 + cos(t)*(i1x-0.5) - sin(t)*(i1y-0.5),
593 : 6795244 : ri1y = 0.5 + sin(t)*(i1x-0.5) + cos(t)*(i1y-0.5),
594 : 6795244 : ri2x = 0.5 + cos(t)*(i2x-0.5) - sin(t)*(i2y-0.5),
595 : 6795244 : ri2y = 0.5 + sin(t)*(i2x-0.5) + cos(t)*(i2y-0.5),
596 : 6795244 : ri3x = 0.5 + cos(t)*(i3x-0.5) - sin(t)*(i3y-0.5),
597 : 6795244 : ri3y = 0.5 + sin(t)*(i3x-0.5) + cos(t)*(i3y-0.5);
598 : :
599 : : // direction of slot sides
600 : 6795244 : tk::real v1x = ri2x-ri1x, v1y = ri2y-ri1y,
601 : 6795244 : v2x = ri3x-ri2x, v2y = ri3y-ri2y;
602 : :
603 : : // lengths of direction of slot sides vectors
604 : 6795244 : tk::real v1 = sqrt(v1x*v1x + v1y*v1y),
605 : 6795244 : v2 = sqrt(v2x*v2x + v2y*v2y);
606 : :
607 : : // cone
608 : 6795244 : r = sqrt((x-kx)*(x-kx) + (y-ky)*(y-ky)) / R0;
609 [ + + ]: 6795244 : if (r<1.0) u[sc] = 0.6*(1.0-r);
610 : :
611 : : // hump
612 : 6795244 : r = sqrt((x-hx)*(x-hx) + (y-hy)*(y-hy)) / R0;
613 [ + + ][ - + ]: 6795244 : if (r<1.0) u[sc] = 0.2*(1.0+cos(M_PI*std::min(r,1.0)));
614 : :
615 : : // cylinder
616 : 6795244 : r = sqrt((x-cx)*(x-cx) + (y-cy)*(y-cy)) / R0;
617 : : const std::array< tk::real, 2 > r1{{ v1x, v1y }},
618 : 6795244 : r2{{ x-ri1x, y-ri1y }};
619 : 6795244 : const auto d1 = (r1[0]*r2[1] - r2[0]*r1[1]) / v1;
620 : : const std::array< tk::real, 2 > r3{{ v2x, v2y }},
621 : 6795244 : r4{{ x-ri2x, y-ri2y }};
622 : 6795244 : const auto d2 = (r3[0]*r4[1] - r4[0]*r3[1]) / v2;
623 [ + + ][ + + ]: 6795244 : if (r<1.0 && (d1>0.05 || d1<0.0 || d2<0.0)) u[sc] = 0.6;
[ + + ][ + + ]
624 : :
625 : 6795244 : return u;
626 : : }
627 : :
628 : : static std::vector< tk::real >
629 : 3777330 : src( tk::real x, tk::real y, tk::real z, tk::real t, std::size_t meshid )
630 : : // *****************************************************************************
631 : : //! Compute and return source term for slotted cylinder, cone, Gauss hump
632 : : //! \param[in] x X coordinate where to evaluate the source
633 : : //! \param[in] y Y coordinate where to evaluate the source
634 : : //! \param[in] z Z coordinate where to evaluate the source
635 : : //! \param[in] t Time where to evaluate the source
636 : : //! \param[in] meshid Mesh id to use
637 : : //! \return Source for flow variables + transported scalars
638 : : // *****************************************************************************
639 : : {
640 : : // evaluate solution at x,y,z,t
641 : 3777330 : auto u = ic( x, y, z, t, meshid );
642 : :
643 : : // configure number of scalar components
644 : : std::size_t ncomp = 6;
645 : : const auto& solver = g_cfg.get< tag::solver >();
646 [ + + ]: 3777330 : if (solver == "chocg") {
647 : : ncomp = 4;
648 : : }
649 [ + + ]: 3312110 : else if (solver == "lohcg") {
650 : : ncomp = 5;
651 : : }
652 : :
653 [ + - ][ - - ]: 3777330 : std::vector< tk::real > s( ncomp, 0.0 );
654 : :
655 : : // momentum source
656 [ + + ]: 3777330 : if (solver == "chocg") {
657 : 465220 : s[0] = -u[1];
658 : 465220 : s[1] = u[0];
659 : : }
660 [ + + ]: 3312110 : else if (solver == "lohcg") {
661 : 369140 : s[1] = -u[2];
662 : 369140 : s[2] = u[1];
663 : : }
664 : : else {
665 : 2942970 : s[1] = -u[2];
666 : 2942970 : s[2] = u[1];
667 : : }
668 : :
669 : 3777330 : return s;
670 : : }
671 : :
672 : : } // slot_cyl::
673 : :
674 : : namespace sheardiff {
675 : :
676 : : static std::vector< tk::real >
677 : 0 : ic( tk::real x, tk::real y, tk::real, tk::real t, std::size_t )
678 : : // *****************************************************************************
679 : : //! Set initial conditions prescribing shear-diffusion in 2D
680 : : //! \param[in] x X coordinate where to evaluate the solution
681 : : //! \param[in] y Y coordinate where to evaluate the solution
682 : : //! \param[in] t Time where to evaluate the solution
683 : : //! \return Values of conserved variables
684 : : //! \see A. Okubo, M.J. Karweit, Diffusion from a continuous source in a uniform
685 : : //! shear flow, Limnology and Oceanography, 14, 1969,
686 : : //! https://doi.org/10.4319/lo.1969.14.4.0514.
687 : : // *****************************************************************************
688 : : {
689 : : using std::exp;
690 : : using std::pow;
691 : : using std::sqrt;
692 : :
693 : : // manufactured solution parameters
694 : 0 : auto V0 = g_cfg.get< tag::problem_p0 >(); // translation velocity in x
695 : 0 : auto L = g_cfg.get< tag::problem_alpha >();// shear velocity
696 : 0 : auto t0 = g_cfg.get< tag::t0 >(); // initial time
697 : 0 : auto dif = g_cfg.get< tag::mat_dyn_diffusivity >(); // diffusivity
698 : : auto eps = std::numeric_limits< tk::real >::epsilon();
699 [ - - ][ - - ]: 0 : if (dif < eps) Throw( "Diffusivity must be positive" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
700 : :
701 : : // configure number of scalar components
702 : : std::size_t ncomp = 4;
703 : : const auto& solver = g_cfg.get< tag::solver >();
704 [ - - ]: 0 : if (solver == "chocg") {
705 : : } else
706 [ - - ]: 0 : if (solver == "lohcg") {
707 : : ncomp = 5;
708 : : }
709 [ - - ][ - - ]: 0 : else Throw( "Shear-diff IC not setup for this solver" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
710 : :
711 : : // prescribed velocity
712 : 0 : std::vector< tk::real > u( ncomp, 0.0 );
713 : : std::size_t sc = 3;
714 [ - - ]: 0 : if (solver == "chocg") {
715 : 0 : u[0] = V0 + L*y;
716 : 0 : u[1] = u[2] = 0.0;
717 : : }
718 [ - - ]: 0 : else if (solver == "lohcg") {
719 : 0 : u[1] = V0 + L*y;
720 : 0 : u[2] = u[3] = 0.0;
721 : : sc = 4;
722 : : }
723 : :
724 : 0 : auto M = [=](tk::real T){ return 4.0*M_PI*T*sqrt(1.0 + L*L*T*T/12.0); };
725 : :
726 : 0 : u[sc] = M(t0) / M(t) *
727 : 0 : exp( -pow( x - V0*t - 0.5*L*y*t, 2.0 ) /
728 : 0 : (4.0*dif*t*(1.0 + L*L*t*t/12.0)) - y*y / (4.0*dif*t) );
729 : :
730 : 0 : return u;
731 : : }
732 : :
733 : : static std::vector< tk::real >
734 : 0 : src( tk::real, tk::real, tk::real, tk::real, std::size_t )
735 : : // *****************************************************************************
736 : : //! Compute and return source term for shear-diffusion in 2D
737 : : //! \return Source for flow variables + transported scalars
738 : : //! \see A. Okubo, M.J. Karweit, Diffusion from a continuous source in a uniform
739 : : //! shear flow, Limnology and Oceanography, 14, 1969,
740 : : //! https://doi.org/10.4319/lo.1969.14.4.0514.
741 : : // *****************************************************************************
742 : : {
743 : : // configure number of scalar components
744 : : std::size_t ncomp = 4;
745 : : const auto& solver = g_cfg.get< tag::solver >();
746 [ - - ]: 0 : if (solver == "lohcg") {
747 : : ncomp = 5;
748 : : }
749 : :
750 : : // manufactured solution parameters
751 : 0 : auto L = g_cfg.get< tag::problem_alpha >();// shear velocity
752 : :
753 : : // source
754 : 0 : std::vector< tk::real > s( ncomp, 0.0 );
755 [ - - ]: 0 : if (solver == "lohcg") {
756 : 0 : s[0] = -L;
757 : : }
758 : :
759 : 0 : return s;
760 : : }
761 : :
762 : : } // sheardiff::
763 : :
764 : : namespace point_src {
765 : :
766 : : static void
767 [ + - ]: 620 : src( const std::array< std::vector< tk::real >, 3 >& coord,
768 : : tk::real t,
769 : : tk::Fields& U )
770 : : // *****************************************************************************
771 : : //! Apply point-source directly to numerical solution
772 : : //! \param[in] coord Mesh node coordinates
773 : : //! \param[in] t Physical time
774 : : //! \param[in,out] U Solution vector at recent time step
775 : : //! \note This is different from other source terms, because this directly
776 : : //! modifies the solution instead of applied as a source term mathematically.
777 : : //! Hence the function signature is also different.
778 : : // *****************************************************************************
779 : : {
780 [ + - ]: 620 : if (U.nprop() == 5) return;
781 : :
782 : : const auto& source = g_cfg.get< tag::problem_src >();
783 : : const auto& location = source.get< tag::location >();
784 : 620 : auto radius = source.get< tag::radius >();
785 [ + - ]: 620 : auto release_time = source.get< tag::release_time >();
786 : : auto largereal = std::numeric_limits< double >::max();
787 : :
788 [ + - ]: 620 : if (location.size() != 3 ||
789 [ + - ][ + - ]: 620 : std::abs(radius - largereal) < 1.0e-12 ||
[ + - ]
790 [ + - ]: 620 : std::abs(release_time - largereal) < 1.0e-12)
791 : : {
792 : : return;
793 : : }
794 : :
795 : 620 : auto sx = location[0];
796 : 620 : auto sy = location[1];
797 : 620 : auto sz = location[2];
798 : : auto sr = radius;
799 : : auto st = release_time;
800 : :
801 [ + - ]: 620 : if (t < st) return;
802 : :
803 : : //configure scalar location
804 : : const auto& solver = g_cfg.get< tag::solver >();
805 : : std::size_t sc = 5;
806 [ + + ]: 620 : if(solver == "chocg"){
807 : : sc = 3;
808 [ - + ]: 600 : } else if(solver =="lohcg"){
809 : : sc = 4;
810 : : }
811 : : const auto& x = coord[0];
812 : : const auto& y = coord[1];
813 : : const auto& z = coord[2];
814 : :
815 [ + + ]: 299660 : for (std::size_t i=0; i<U.nunk(); ++i) {
816 : 299040 : auto rx = sx - x[i];
817 : 299040 : auto ry = sy - y[i];
818 : 299040 : auto rz = sz - z[i];
819 [ + + ]: 299040 : if (rx*rx + ry*ry + rz*rz < sr*sr) U(i,sc) = 1.0;
820 : : }
821 : :
822 : : return;
823 : : }
824 : :
825 : : } // point_src::
826 : :
827 : : namespace poisson {
828 : :
829 : : static std::vector< tk::real >
830 : 8524 : ic( tk::real, tk::real, tk::real, tk::real, std::size_t )
831 : : // *****************************************************************************
832 : : //! Set velocity initial conditions for testing a Poisson solve only
833 : : //! \return Values for initial conditions
834 : : // *****************************************************************************
835 : : {
836 : 8524 : return { 0, 0, 0 };
837 : : }
838 : :
839 : : } // poisson::
840 : :
841 : : namespace poisson_const {
842 : :
843 : : static tk::real
844 : 979 : pr( tk::real, tk::real, tk::real )
845 : : // *****************************************************************************
846 : : //! Set pressure rhs for testing a Poisson solve
847 : : //! \return Value for pressure rhs
848 : : // *****************************************************************************
849 : : {
850 : 979 : return 6.0;
851 : : }
852 : :
853 : : static tk::real
854 : 3739 : ic( tk::real x, tk::real y, tk::real z, std::size_t )
855 : : // *****************************************************************************
856 : : //! Evaluate pressure boundary condition
857 : : //! \param[in] x X coordinate where to evaluate the BC
858 : : //! \param[in] y Y coordinate where to evaluate the BC
859 : : //! \param[in] z Z coordinate where to evaluate the BC
860 : : //! \return Value for pressure BC
861 : : // *****************************************************************************
862 : : {
863 : 3739 : return x*x + y*y + z*z;
864 : : }
865 : :
866 : : } // poisson_const::
867 : :
868 : : namespace poisson_harmonic {
869 : :
870 : : static tk::real
871 : 0 : pr( tk::real, tk::real, tk::real )
872 : : // *****************************************************************************
873 : : //! Set pressure rhs for testing a Laplace solve
874 : : //! \return Value for pressure rhs
875 : : // *****************************************************************************
876 : : {
877 : 0 : return 0.0;
878 : : }
879 : :
880 : : static tk::real
881 : 0 : ic( tk::real x, tk::real y, tk::real z, std::size_t )
882 : : // *****************************************************************************
883 : : //! Evaluate pressure boundary condition
884 : : //! \param[in] x X coordinate where to evaluate the BC
885 : : //! \param[in] y Y coordinate where to evaluate the BC
886 : : //! \param[in] z Z coordinate where to evaluate the BC
887 : : //! \return Value for pressure BC
888 : : // *****************************************************************************
889 : : {
890 : : const auto& b = g_cfg.get< tag::problem_beta >();
891 : 0 : auto x0 = b[0];
892 : 0 : auto y0 = b[1];
893 : 0 : auto z0 = b[2];
894 : :
895 : 0 : return 1.0 / std::sqrt( (x-x0)*(x-x0) + (y-y0)*(y-y0) + (z-z0)*(z-z0) );
896 : : }
897 : :
898 : : } // poisson_harmonic::
899 : :
900 : : namespace poisson_sine {
901 : :
902 : : static tk::real
903 : 278 : pr( tk::real x, tk::real y, tk::real z )
904 : : // *****************************************************************************
905 : : //! Set pressure rhs for testing a Poisson solve
906 : : //! \return Value for pressure rhs
907 : : // *****************************************************************************
908 : : {
909 : 278 : return -M_PI * M_PI * x * y * std::sin( M_PI * z );
910 : : }
911 : :
912 : : static tk::real
913 : 1067 : ic( tk::real x, tk::real y, tk::real z, std::size_t )
914 : : // *****************************************************************************
915 : : //! Evaluate pressure boundary condition
916 : : //! \param[in] x X coordinate where to evaluate the BC
917 : : //! \param[in] y Y coordinate where to evaluate the BC
918 : : //! \param[in] z Z coordinate where to evaluate the BC
919 : : //! \return Value for pressure BC
920 : : // *****************************************************************************
921 : : {
922 : 1067 : return x * y * std::sin( M_PI * z );
923 : : }
924 : :
925 : : } // poisson_sine::
926 : :
927 : : namespace poisson_sine3 {
928 : :
929 : : static tk::real
930 : 596 : pr( tk::real x, tk::real y, tk::real z )
931 : : // *****************************************************************************
932 : : //! Set pressure rhs for testing a Poisson solve
933 : : //! \return Value for pressure rhs
934 : : // *****************************************************************************
935 : : {
936 : : using std::sin;
937 : :
938 : 596 : return -3.0 * M_PI * M_PI * sin(M_PI*x) * sin(M_PI*y) * sin(M_PI*z);
939 : : }
940 : :
941 : : static tk::real
942 : 2281 : ic( tk::real x, tk::real y, tk::real z, std::size_t )
943 : : // *****************************************************************************
944 : : //! Evaluate pressure boundary condition
945 : : //! \param[in] x X coordinate where to evaluate the BC
946 : : //! \param[in] y Y coordinate where to evaluate the BC
947 : : //! \param[in] z Z coordinate where to evaluate the BC
948 : : //! \return Value for pressure BC
949 : : // *****************************************************************************
950 : : {
951 : 2281 : return sin(M_PI*x) * sin(M_PI*y) * sin(M_PI*z);
952 : : }
953 : :
954 : : } // poisson_sine3::
955 : :
956 : : namespace poisson_neumann {
957 : :
958 : : static tk::real
959 : 278 : pr( tk::real x, tk::real y, tk::real )
960 : : // *****************************************************************************
961 : : //! Set pressure rhs for testing a Poisson solve
962 : : //! \param[in] x X coordinate where to evaluate the rhs
963 : : //! \param[in] y Y coordinate where to evaluate the rhs
964 : : //! \return Value for pressure rhs
965 : : // *****************************************************************************
966 : : {
967 : 278 : return -3.0 * std::cos(2.0*x) * std::exp(y);
968 : : }
969 : :
970 : : static std::array< tk::real, 3 >
971 : 612 : pg( tk::real x, tk::real y, tk::real )
972 : : // *****************************************************************************
973 : : //! Set pressure gradient for testing a Poisson solve
974 : : //! \param[in] x X coordinate where to evaluate the pressure gradient
975 : : //! \param[in] y Y coordinate where to evaluate the pressure gradient
976 : : //! \return Value for pressure gradient at a point
977 : : // *****************************************************************************
978 : : {
979 : 612 : return { -2.0 * std::sin( 2.0 * x ) * std::exp( y ),
980 : 612 : std::cos(2.0*x) * std::exp(y),
981 : 612 : 0.0 };
982 : : }
983 : :
984 : : static tk::real
985 : 972 : ic( tk::real x, tk::real y, tk::real, std::size_t )
986 : : // *****************************************************************************
987 : : //! Evaluate pressure boundary condition
988 : : //! \param[in] x X coordinate where to evaluate the IC / analytic solution
989 : : //! \param[in] y Y coordinate where to evaluate the IC / analytic solution
990 : : //! \return Value for pressure
991 : : // *****************************************************************************
992 : : {
993 : 972 : return std::cos(2.0*x) * std::exp(y);
994 : : }
995 : :
996 : : } // poisson_neumann::
997 : :
998 : :
999 : : namespace poiseuille {
1000 : :
1001 : : static std::vector< tk::real >
1002 : 0 : ic( tk::real, tk::real y, tk::real, tk::real, std::size_t )
1003 : : // *****************************************************************************
1004 : : //! Set initial conditions prescribing the Poisuille problem
1005 : : //! \param[in] y Y coordinate where to evaluate the solution
1006 : : //! \return Values of physics variables
1007 : : // *****************************************************************************
1008 : : {
1009 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1010 : 0 : auto nu = g_cfg.get< tag::mat_dyn_viscosity >();
1011 [ - - ][ - - ]: 0 : if (nu < eps) Throw( "Poiseuille flow needs nonzero viscosity" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
1012 : :
1013 : : auto dpdx = -0.12;
1014 : 0 : auto u = -dpdx * y * (1.0 - y) / 2.0 / nu;
1015 : :
1016 : : const auto& solver = g_cfg.get< tag::solver >();
1017 [ - - ]: 0 : if (solver == "chocg") {
1018 : 0 : return { u, 0.0, 0.0 };
1019 : : } else
1020 [ - - ]: 0 : if (solver == "lohcg") {
1021 : 0 : return { 0.0, 0.0, 0.0, 0.0 };
1022 : : }
1023 [ - - ][ - - ]: 0 : else Throw( "Poiseuille IC not setup for this solver" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1024 : : }
1025 : :
1026 : : static std::vector< tk::real >
1027 : 0 : sol( tk::real, tk::real y, tk::real, tk::real, std::size_t )
1028 : : // *****************************************************************************
1029 : : //! Set analytic solution of the Poisuille problem
1030 : : //! \param[in] y Y coordinate where to evaluate the solution
1031 : : //! \return Values of physics variables
1032 : : // *****************************************************************************
1033 : : {
1034 : : auto eps = std::numeric_limits< tk::real >::epsilon();
1035 : 0 : auto nu = g_cfg.get< tag::mat_dyn_viscosity >();
1036 [ - - ][ - - ]: 0 : if (nu < eps) Throw( "Poiseuille flow needs nonzero viscosity" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ][ - - ]
1037 : :
1038 : : auto dpdx = -0.12;
1039 : 0 : auto u = -dpdx * y * (1.0 - y) / 2.0 / nu;
1040 : :
1041 : : const auto& solver = g_cfg.get< tag::solver >();
1042 [ - - ]: 0 : if (solver == "chocg") {
1043 : 0 : return { u, 0.0, 0.0 };
1044 : : } else
1045 [ - - ]: 0 : if (solver == "lohcg") {
1046 : 0 : return { 0.0, u, 0.0, 0.0 };
1047 : : }
1048 [ - - ][ - - ]: 0 : else Throw( "Poiseuille IC not setup for this solver" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1049 : : }
1050 : :
1051 : : } // poiseuille::
1052 : :
1053 : : namespace overset_linear {
1054 : :
1055 : : static std::vector< tk::real >
1056 : 0 : ic( tk::real x, tk::real, tk::real, tk::real, std::size_t meshid )
1057 : : // *****************************************************************************
1058 : : //! Set initial conditions prescribing linear a field to test overset
1059 : : //! \param[in] x X coordinate where to set the solution
1060 : : //! \param[in] meshid Mesh id to use
1061 : : //! \return Values of physics variables
1062 : : // *****************************************************************************
1063 : : {
1064 [ - - ]: 0 : return { 0.0, meshid == 0 ? x : 0.0, 0.0, 0.0 };
1065 : : }
1066 : :
1067 : : } // overset_linear::
1068 : :
1069 : : std::function< std::vector< tk::real >
1070 : : ( tk::real, tk::real, tk::real, tk::real, std::size_t ) >
1071 : 117707 : IC()
1072 : : // *****************************************************************************
1073 : : // Query user config and assign function to set initial conditions
1074 : : //! \return The function to call to set initial conditions
1075 : : // *****************************************************************************
1076 : : {
1077 : : const auto& problem = inciter::g_cfg.get< tag::problem >();
1078 : :
1079 : : std::function< std::vector< tk::real >
1080 : : ( tk::real, tk::real, tk::real, tk::real, std::size_t ) > ic;
1081 : :
1082 [ + + ][ + + ]: 186014 : if (problem == "userdef" or problem == "point_src")
1083 : 50051 : ic = userdef::ic;
1084 [ + + ]: 67656 : else if (problem == "nonlinear_energy_growth")
1085 : 1254 : ic = nonlinear_energy_growth::ic;
1086 [ + + ]: 66402 : else if (problem == "rayleigh_taylor")
1087 : 3624 : ic = rayleigh_taylor::ic;
1088 [ + + ]: 62778 : else if (problem == "sedov")
1089 : 1358 : ic = sedov::ic;
1090 [ + + ]: 61420 : else if (problem == "sod")
1091 : 4468 : ic = sod::ic;
1092 [ + + ]: 56952 : else if (problem == "taylor_green")
1093 : 4812 : ic = taylor_green::ic;
1094 [ + + ]: 52140 : else if (problem == "vortical_flow")
1095 : 29151 : ic = vortical_flow::ic;
1096 [ + + ]: 22989 : else if (problem == "slot_cyl")
1097 : 22790 : ic = slot_cyl::ic;
1098 [ - + ]: 199 : else if (problem == "sheardiff")
1099 : 0 : ic = sheardiff::ic;
1100 [ + - ]: 199 : else if (problem.find("poisson") != std::string::npos)
1101 : 199 : ic = poisson::ic;
1102 [ - - ]: 0 : else if (problem == "poiseuille")
1103 : 0 : ic = poiseuille::ic;
1104 [ - - ]: 0 : else if (problem == "overset_linear")
1105 : 0 : ic = overset_linear::ic;
1106 : : else
1107 [ - - ][ - - ]: 0 : Throw( "problem type ic not hooked up" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1108 : :
1109 : 117707 : return ic;
1110 : : }
1111 : :
1112 : : std::function< std::vector< tk::real >
1113 : : ( tk::real, tk::real, tk::real, tk::real, std::size_t ) >
1114 : 47292 : SOL()
1115 : : // *****************************************************************************
1116 : : // Query user config and assign function to query analytic solutions
1117 : : //! \return The function to call to query analytic solutions
1118 : : // *****************************************************************************
1119 : : {
1120 : : const auto& problem = inciter::g_cfg.get< tag::problem >();
1121 : :
1122 [ + + ]: 24514 : if (problem == "userdef" or
1123 [ + + ]: 22099 : problem == "sod" or
1124 [ + + ][ + + ]: 67377 : problem == "sedov" or
1125 : : problem == "point_src")
1126 : : return {};
1127 [ - + ]: 20015 : else if (problem == "poiseuille")
1128 : : return poiseuille::sol;
1129 : : else
1130 : 20015 : return IC();
1131 : : }
1132 : :
1133 : : void
1134 : 2788 : initialize( const std::array< std::vector< tk::real >, 3 >& coord,
1135 : : tk::Fields& U,
1136 : : tk::real t,
1137 : : std::size_t meshid,
1138 : : const std::vector< std::unordered_set< std::size_t > >& boxnodes )
1139 : : // *****************************************************************************
1140 : : // Set inital conditions
1141 : : //! \param[in] coord Mesh node coordinates
1142 : : //! \param[in,out] U Array of unknowns
1143 : : //! \param[in] t Physical time
1144 : : //! \param[in] meshid Mesh id to use
1145 : : //! \param[in] boxnodes Nodes at which box user ICs are set (for each box IC)
1146 : : // *****************************************************************************
1147 : : {
1148 : : Assert( coord[0].size() == U.nunk(), "Size mismatch" );
1149 : :
1150 : 2788 : auto ic = IC();
1151 : : const auto& x = coord[0];
1152 : : const auto& y = coord[1];
1153 : : const auto& z = coord[2];
1154 : :
1155 : : // Set initial conditions dependeing on problem configured
1156 [ + + ]: 387695 : for (std::size_t i=0; i<x.size(); ++i) {
1157 : :
1158 : : // Set background ICs
1159 [ - + ]: 384907 : auto s = ic( x[i], y[i], z[i], t, meshid );
1160 : : Assert( s.size() == U.nprop(), "Size mismatch" );
1161 : :
1162 : : // Initialize user-defined ICs in boxes
1163 [ + - ]: 384907 : box( i, s, boxnodes );
1164 : :
1165 : : // Set values for ICs
1166 [ + + ]: 2334777 : for (std::size_t c=0; c<s.size(); ++c) U(i,c) = s[c];
1167 : :
1168 : : }
1169 : 2788 : }
1170 : :
1171 : : std::function< tk::real( tk::real, tk::real, tk::real ) >
1172 : 6242 : PRESSURE_RHS()
1173 : : // *****************************************************************************
1174 : : // Query user config and assign function to set pressure rhs
1175 : : //! \return The function to call to set pressure rhs
1176 : : // *****************************************************************************
1177 : : {
1178 : : const auto& problem = inciter::g_cfg.get< tag::problem >();
1179 : :
1180 : : std::function< tk::real( tk::real, tk::real, tk::real ) > pr;
1181 : :
1182 [ + + ]: 6242 : if (problem == "poisson_const")
1183 : 22 : pr = poisson_const::pr;
1184 [ - + ]: 6220 : else if (problem == "poisson_harmonic")
1185 : 0 : pr = poisson_harmonic::pr;
1186 [ + + ]: 6220 : else if (problem == "poisson_sine")
1187 : 2 : pr = poisson_sine::pr;
1188 [ + + ]: 6218 : else if (problem == "poisson_sine3")
1189 : 6 : pr = poisson_sine3::pr;
1190 [ + + ]: 6212 : else if (problem == "poisson_neumann")
1191 : 2 : pr = poisson_neumann::pr;
1192 : :
1193 : 6242 : return pr;
1194 : : }
1195 : :
1196 : : std::function< tk::real( tk::real, tk::real, tk::real, std::size_t ) >
1197 : 18316 : PRESSURE_IC()
1198 : : // *****************************************************************************
1199 : : // Query user config and assign function to set pressure initial conditions
1200 : : //! \return The function to call to set pressure initial conditions
1201 : : // *****************************************************************************
1202 : : {
1203 : : const auto& problem = inciter::g_cfg.get< tag::problem >();
1204 : :
1205 : : std::function< tk::real( tk::real, tk::real, tk::real, std::size_t ) > ic;
1206 : :
1207 [ + + ]: 6322 : if ( problem == "userdef" or
1208 [ + - ]: 164 : problem == "slot_cyl" or
1209 [ + - ]: 164 : problem == "sheardiff" or
1210 [ + + ][ + + ]: 18480 : problem == "poiseuille" or
1211 : : problem == "point_src"
1212 : : )
1213 : :
1214 : 18174 : ic = userdef::pic;
1215 [ + + ]: 142 : else if (problem == "poisson_const")
1216 : 94 : ic = poisson_const::ic;
1217 [ - + ]: 48 : else if (problem == "poisson_harmonic")
1218 : 0 : ic = poisson_harmonic::ic;
1219 [ + + ]: 48 : else if (problem == "poisson_sine")
1220 : 10 : ic = poisson_sine::ic;
1221 [ + + ]: 38 : else if (problem == "poisson_sine3")
1222 : 28 : ic = poisson_sine3::ic;
1223 [ + - ]: 10 : else if (problem == "poisson_neumann")
1224 : 10 : ic = poisson_neumann::ic;
1225 : : else
1226 [ - - ][ - - ]: 0 : Throw( "problem type not hooked up" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1227 : :
1228 : 18316 : return ic;
1229 : : }
1230 : :
1231 : : std::function< tk::real( tk::real, tk::real, tk::real, std::size_t ) >
1232 : 5669 : PRESSURE_SOL()
1233 : : // *****************************************************************************
1234 : : // Query user config and assign function to query analytic pressure solutions
1235 : : //! \return The function to call to query analytic solutions
1236 : : // *****************************************************************************
1237 : : {
1238 : : const auto& problem = inciter::g_cfg.get< tag::problem >();
1239 : :
1240 [ + + ]: 1351 : if ( problem == "userdef" or
1241 [ + - ]: 124 : problem == "slot_cyl" or
1242 [ + + ]: 124 : problem == "poiseuille" or
1243 [ + + ][ - + ]: 5779 : problem == "point_src" or
1244 : : problem == "sheardiff" )
1245 : : {
1246 : : return {};
1247 : : }
1248 : : else {
1249 : 110 : return PRESSURE_IC();
1250 : : }
1251 : : }
1252 : :
1253 : : std::function< std::array< tk::real, 3 >( tk::real, tk::real, tk::real ) >
1254 : 6242 : PRESSURE_GRAD()
1255 : : // *****************************************************************************
1256 : : // Assign function to query pressure gradient at a point
1257 : : //! \return The function to call to query the pressure gradient
1258 : : // *****************************************************************************
1259 : : {
1260 : : const auto& problem = inciter::g_cfg.get< tag::problem >();
1261 : :
1262 [ + + ]: 6242 : if (problem == "poisson_neumann")
1263 : : return poisson_neumann::pg;
1264 : :
1265 : : return {};
1266 : : }
1267 : :
1268 : : tk::real
1269 : 0 : initialize( tk::real x, tk::real y, tk::real z, std::size_t meshid )
1270 : : // *****************************************************************************
1271 : : // Evaluate initial condition for pressure
1272 : : //! \param[in] x X coordinate where to evaluate the pressure initial condition
1273 : : //! \param[in] y Y coordinate where to evaluate the pressure initial condition
1274 : : //! \param[in] z Z coordinate where to evaluate the pressure initial condition
1275 : : //! \param[in] meshid Mesh id to use
1276 : : //! \return Pressure initial condition
1277 : : // *****************************************************************************
1278 : : {
1279 : : const auto& problem = inciter::g_cfg.get< tag::problem >();
1280 : :
1281 : : std::function< tk::real( tk::real, tk::real, tk::real, std::size_t ) > ic;
1282 : :
1283 [ - - ]: 0 : if (problem == "poisson_const")
1284 : 0 : ic = poisson_const::ic;
1285 [ - - ]: 0 : else if (problem == "poisson_harmonic")
1286 : 0 : ic = poisson_harmonic::ic;
1287 [ - - ]: 0 : else if (problem == "poisson_sine")
1288 : 0 : ic = poisson_sine::ic;
1289 [ - - ]: 0 : else if (problem == "poisson_sine3")
1290 : 0 : ic = poisson_sine3::ic;
1291 : : else
1292 [ - - ][ - - ]: 0 : Throw( "problem type not hooked up" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
1293 : :
1294 [ - - ]: 0 : return ic( x, y, z, meshid );
1295 : : }
1296 : :
1297 : : std::function< std::vector< tk::real >
1298 : : ( tk::real, tk::real, tk::real, tk::real, std::size_t ) >
1299 : 78256 : SRC()
1300 : : // *****************************************************************************
1301 : : // Query user config and assign function to add a source term
1302 : : //! \return The function to call to evaluate a problem-sepcific source term
1303 : : // *****************************************************************************
1304 : : {
1305 : : const auto& problem = inciter::g_cfg.get< tag::problem >();
1306 : :
1307 : : std::function<
1308 : : std::vector< tk::real >
1309 : : ( tk::real, tk::real, tk::real, tk::real, std::size_t ) > src;
1310 : :
1311 [ + + ]: 78256 : if (problem == "nonlinear_energy_growth")
1312 : 676 : src = nonlinear_energy_growth::src;
1313 [ + + ]: 77580 : else if (problem == "rayleigh_taylor")
1314 : 2000 : src = rayleigh_taylor::src;
1315 [ + + ]: 75580 : else if (problem == "taylor_green")
1316 : 3536 : src = taylor_green::src;
1317 [ + + ]: 72044 : else if (problem == "vortical_flow")
1318 : 18547 : src = vortical_flow::src;
1319 [ + + ]: 53497 : else if (problem == "slot_cyl")
1320 : 11617 : src = slot_cyl::src;
1321 [ - + ]: 41880 : else if (problem == "sheardiff")
1322 : 0 : src = sheardiff::src;
1323 : :
1324 : 78256 : return src;
1325 : : }
1326 : :
1327 : : std::function< void( const std::array< std::vector< tk::real >, 3 >&,
1328 : : tk::real,
1329 : : tk::Fields& ) >
1330 : 78256 : PHYS_SRC()
1331 : : // *****************************************************************************
1332 : : // Query user config and assign function to apply source to numerical solution
1333 : : //! \return The function to call to evaluate a problem-sepcific source term
1334 : : // *****************************************************************************
1335 : : {
1336 : : const auto& problem = inciter::g_cfg.get< tag::problem >();
1337 : :
1338 : : std::function< void( const std::array< std::vector< tk::real >, 3 >&,
1339 : : tk::real,
1340 : : tk::Fields& ) > src;
1341 : :
1342 [ + + ]: 78256 : if (problem == "point_src") {
1343 : 620 : src = point_src::src;
1344 : : }
1345 : :
1346 : 78256 : return src;
1347 : : }
1348 : :
1349 : : } // problems::
|
