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