1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144 | // *****************************************************************************
/*!
\file src/PDE/CompFlow/Problem/RayleighTaylor.hpp
\copyright 2012-2015 J. Bakosi,
2016-2018 Los Alamos National Security, LLC.,
2019-2021 Triad National Security, LLC.
All rights reserved. See the LICENSE file for details.
\brief Problem configuration for the compressible flow equations
\details This file defines a policy class for the compressible flow
equations, defined in PDE/CompFlow/CompFlow.h. See PDE/CompFlow/Problem.h
for general requirements on Problem policy classes for CompFlow.
*/
// *****************************************************************************
#ifndef CompFlowProblemRayleighTaylor_h
#define CompFlowProblemRayleighTaylor_h
#include <string>
#include <unordered_set>
#include "Types.hpp"
#include "Fields.hpp"
#include "FunctionPrototypes.hpp"
#include "SystemComponents.hpp"
#include "Inciter/Options/Problem.hpp"
#include "Inciter/InputDeck/InputDeck.hpp"
#include "EoS/EoS.hpp"
namespace inciter {
extern ctr::InputDeck g_inputdeck;
//! CompFlow system of PDEs problem: Rayleigh-Taylor
//! \see Waltz, et. al, "Manufactured solutions for the three-dimensional Euler
//! equations with relevance to Inertial Confinement Fusion", Journal of
//! Computational Physics 267 (2014) 196-209.
class CompFlowProblemRayleighTaylor {
private:
using ncomp_t = tk::ctr::ncomp_t;
using eq = tag::compflow;
public:
//! Initialize numerical solution
static tk::InitializeFn::result_type
initialize( ncomp_t system, ncomp_t, tk::real x, tk::real y,
tk::real z, tk::real t );
//! Evaluate analytical solution at (x,y,z,t) for all components
static tk::InitializeFn::result_type
analyticSolution( ncomp_t system, ncomp_t, tk::real x, tk::real y,
tk::real z, tk::real t );
//! Compute and return source term for Rayleigh-Taylor manufactured solution
//! \param[in] system Equation system index, i.e., which compressible
//! flow equation system we operate on among the systems of PDEs
//! \param[in] x X coordinate where to evaluate the solution
//! \param[in] y Y coordinate where to evaluate the solution
//! \param[in] z Z coordinate where to evaluate the solution
//! \param[in] t Physical time at which to evaluate the source
//! \param[in,out] r Density source
//! \param[in,out] ru X momentum source
//! \param[in,out] rv Y momentum source
//! \param[in,out] rw Z momentum source
//! \param[in,out] re Specific total energy source
//! \note The function signature must follow tk::SrcFn
static tk::CompFlowSrcFn::result_type
src( ncomp_t system, tk::real x, tk::real y, tk::real z, tk::real t,
tk::real& r, tk::real& ru, tk::real& rv, tk::real& rw, tk::real& re )
{
using tag::param; using std::sin; using std::cos;
// manufactured solution parameters
auto a = g_inputdeck.get< param, eq, tag::alpha >()[system];
auto bx = g_inputdeck.get< param, eq, tag::betax >()[system];
auto by = g_inputdeck.get< param, eq, tag::betay >()[system];
auto bz = g_inputdeck.get< param, eq, tag::betaz >()[system];
auto k = g_inputdeck.get< param, eq, tag::kappa >()[system];<--- Shadow variable
auto p0 = g_inputdeck.get< param, eq, tag::p0 >()[system];
// ratio of specific heats
auto g = gamma< tag::compflow >(system);
// evaluate solution at x,y,z,t
auto s = initialize( system, 5, x, y, z, t );
// density, velocity, energy, pressure
auto rho = s[0];
auto u = s[1]/s[0];
auto v = s[2]/s[0];
auto w = s[3]/s[0];
auto E = s[4]/s[0];
auto p = p0 + a*(bx*x*x + by*y*y + bz*z*z);
// spatial gradients
std::array< tk::real, 3 > drdx{{ -2.0*bx*x, -2.0*by*y, -2.0*bz*z }};
std::array< tk::real, 3 > dpdx{{ 2.0*a*bx*x, 2.0*a*by*y, 2.0*a*bz*z }};
tk::real ft = cos(k*M_PI*t);
std::array< tk::real, 3 > dudx{{ ft*M_PI*z*cos(M_PI*x),
0.0,
ft*sin(M_PI*x) }};
std::array< tk::real, 3 > dvdx{{ 0.0,
-ft*M_PI*z*sin(M_PI*y),
ft*cos(M_PI*y) }};
std::array< tk::real, 3 > dwdx{{ ft*M_PI*0.5*M_PI*z*z*sin(M_PI*x),
ft*M_PI*0.5*M_PI*z*z*cos(M_PI*y),
-ft*M_PI*z*(cos(M_PI*x) - sin(M_PI*y)) }};
std::array< tk::real, 3 > dedx{{
dpdx[0]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[0]
+ u*dudx[0] + v*dvdx[0] + w*dwdx[0],
dpdx[1]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[1]
+ u*dudx[1] + v*dvdx[1] + w*dwdx[1],
dpdx[2]/rho/(g-1.0) - p/(g-1.0)/rho/rho*drdx[2]
+ u*dudx[2] + v*dvdx[2] + w*dwdx[2] }};
// time derivatives
auto dudt = -k*M_PI*sin(k*M_PI*t)*z*sin(M_PI*x);
auto dvdt = -k*M_PI*sin(k*M_PI*t)*z*cos(M_PI*y);
auto dwdt = k*M_PI*sin(k*M_PI*t)/2*M_PI*z*z*(cos(M_PI*x) - sin(M_PI*y));
auto dedt = u*dudt + v*dvdt + w*dwdt;
// density source
r = u*drdx[0] + v*drdx[1] + w*drdx[2];
// momentum source
ru = rho*dudt+u*r+dpdx[0] + s[1]*dudx[0]+s[2]*dudx[1]+s[3]*dudx[2];
rv = rho*dvdt+v*r+dpdx[1] + s[1]*dvdx[0]+s[2]*dvdx[1]+s[3]*dvdx[2];
rw = rho*dwdt+w*r+dpdx[2] + s[1]*dwdx[0]+s[2]*dwdx[1]+s[3]*dwdx[2];
// energy source
re = rho*dedt + E*r + s[1]*dedx[0]+s[2]*dedx[1]+s[3]*dedx[2]
+ u*dpdx[0]+v*dpdx[1]+w*dpdx[2];
}
//! Return field names to be output to file
std::vector< std::string > analyticFieldNames( ncomp_t ) const;
//! Return names of integral variables to be output to diagnostics file
std::vector< std::string > names( ncomp_t /*ncomp*/ ) const;
//! Return problem type
static ctr::ProblemType type() noexcept
{ return ctr::ProblemType::RAYLEIGH_TAYLOR; }
};
} // inciter::
#endif // CompFlowProblemRayleighTaylor_h
|