``` 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 #include #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 ```