Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/CompFlow/Problem/BoxInitialization.cpp
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 User-defined box initialization
9 : : \details This file defines functions for initializing solutions for
10 : : compressible single-material equations inside the user-defined box.
11 : : */
12 : : // *****************************************************************************
13 : :
14 : : #include "BoxInitialization.hpp"
15 : : #include "EoS/EoS.hpp"
16 : : #include "ContainerUtil.hpp"
17 : : #include "Control/Inciter/Types.hpp"
18 : :
19 : : namespace inciter {
20 : :
21 : 759 : void initializeBox( std::size_t system,
22 : : tk::real VRatio,
23 : : tk::real t,
24 : : const inciter::ctr::box& b,
25 : : tk::real bgpreic,
26 : : tk::real cv,
27 : : std::vector< tk::real >& s )
28 : : // *****************************************************************************
29 : : // Set the solution in the user-defined IC box
30 : : //! \param[in] system Equation system index
31 : : //! \param[in] VRatio Ratio of exact box volume to discrete box volume
32 : : //! \param[in] t Physical time
33 : : //! \param[in] b IC box configuration to use
34 : : //! \param[in] bgpreic Background pressure user input
35 : : //! \param[in] cv Specific heats ratio user input
36 : : //! \param[in,out] s Solution vector that is set to box ICs
37 : : //! \details This function sets the fluid density and total specific energy
38 : : //! within a box initial condition, configured by the user. If the user
39 : : //! is specified a box where mass is specified, we also assume here that
40 : : //! internal energy content (energy per unit volume) is also
41 : : //! specified. Specific internal energy (energy per unit mass) is then
42 : : //! computed here (and added to the kinetic energy) from the internal
43 : : //! energy per unit volume by multiplying it with the total box volume
44 : : //! and dividing it by the total mass of the material in the box.
45 : : //! Example (SI) units of the quantities involved:
46 : : //! * internal energy content (energy per unit volume): J/m^3
47 : : //! * specific energy (internal energy per unit mass): J/kg
48 : : // *****************************************************************************
49 : : {
50 : : std::vector< tk::real > box
51 : 759 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
52 : 759 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
53 : 759 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
54 : :
55 : : const auto& initiate = b.template get< tag::initiate >();
56 : 759 : auto inittype = initiate.template get< tag::init >();
57 : :
58 : 759 : auto boxrho = b.template get< tag::density >();
59 : : const auto& boxvel = b.template get< tag::velocity >();
60 : 759 : auto boxpre = b.template get< tag::pressure >();
61 : 759 : auto boxene = b.template get< tag::energy >();
62 : 759 : auto boxtem = b.template get< tag::temperature >();
63 : 759 : auto boxmas = b.template get< tag::mass >();
64 : 759 : auto boxenc = b.template get< tag::energy_content >();
65 : :
66 : : tk::real rho = 0.0, ru = 0.0, rv = 0.0, rw = 0.0, re = 0.0, spi = 0.0;
67 : : bool boxmassic = false;
68 [ + + ]: 759 : if (boxmas > 0.0) {
69 : : Assert( boxenc > 0.0, "Box energy content must be nonzero" );
70 : 268 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
71 : 268 : rho = boxmas / V_ex;
72 : 268 : spi = boxenc * VRatio / rho;
73 : : boxmassic = true;
74 : : } else {
75 [ + - ]: 491 : if (boxrho > 0.0) rho = boxrho;
76 [ - + ]: 491 : if (boxvel.size() == 3) {
77 : 0 : ru = rho * boxvel[0];
78 : 0 : rv = rho * boxvel[1];
79 : 0 : rw = rho * boxvel[2];
80 : : }
81 [ + - ]: 491 : if (boxpre > 0.0) {
82 : : re = eos_totalenergy< tag::compflow >
83 : 491 : ( system, rho, ru/rho, rv/rho, rw/rho, boxpre );
84 : : }
85 [ - + ]: 491 : if (boxene > 0.0) {
86 : 0 : auto ux = ru/rho, uy = rv/rho, uz = rw/rho;
87 : 0 : auto ke = 0.5*(ux*ux + uy*uy + uz*uz);
88 : 0 : re = rho * (boxene + ke);
89 : : }
90 [ - + ]: 491 : if (boxtem > 0.0) re = rho * boxtem * cv;
91 : : }
92 : :
93 : : // Initiate type 'impulse' simply assigns the prescribed values to all
94 : : // nodes within a box.
95 [ + + ]: 759 : if (inittype == ctr::InitiateType::IMPULSE) {
96 : : // superimpose on existing velocity field
97 : 491 : auto u = s[1]/s[0], v = s[2]/s[0], w = s[3]/s[0];
98 : 491 : auto ke = 0.5*(u*u + v*v + w*w);
99 : 491 : s[0] = rho;
100 [ - + ]: 491 : if (boxmassic) {
101 : 0 : s[1] = rho * u;
102 : 0 : s[2] = rho * v;
103 : 0 : s[3] = rho * w;
104 : 0 : s[4] = rho * (spi + ke);
105 : : } else {
106 : 491 : s[1] = ru;
107 : 491 : s[2] = rv;
108 : 491 : s[3] = rw;
109 : 491 : s[4] = re;
110 : : }
111 : : }
112 : : // Initiate type 'linear' assigns the prescribed values to all nodes
113 : : // within a box. This is followed by adding a time-dependent energy
114 : : // source term representing a planar wave-front propagating along the
115 : : // z-direction with a velocity specified in the IC linear...end block.
116 : : // The wave front is smoothed out to include a couple of mesh nodes.
117 : : // see boxSrc() for further details.
118 [ + - ][ + - ]: 268 : else if (inittype == ctr::InitiateType::LINEAR && t < 1e-12) {
119 : :
120 : : // superimpose on existing velocity field
121 : 268 : const auto u = s[1]/s[0],
122 : 268 : v = s[2]/s[0],
123 : 268 : w = s[3]/s[0];
124 : 268 : const auto ke = 0.5*(u*u + v*v + w*w);
125 : :
126 : : // The linear-propagating source initialization can be done only
127 : : // based on background pressure (not on temperature): The IC box can
128 : : // have a different density than the background, while having the
129 : : // same pressure and temperature as the background. This means, the
130 : : // material in the box has a different specific heat (Cv) than the
131 : : // background material. If such a box has to be initialized based on
132 : : // temperature, the Cv of the box will have to be specified
133 : : // separately. This is not currently supported.
134 [ + - ]: 268 : if (bgpreic > 0.0) {
135 : : // energy based on box density and background pressure
136 : 268 : spi = eos_totalenergy< tag::compflow >
137 : 268 : ( system, rho, u, v, w, bgpreic ) / rho;
138 [ - - ][ - - ]: 0 : } else Throw( "Background pressure must be specified for box-IC "
[ - - ][ - - ]
[ - - ][ - - ]
139 : : "with linear propagating source");
140 : :
141 : 268 : s[0] = rho;
142 : 268 : s[1] = rho * u;
143 : 268 : s[2] = rho * v;
144 : 268 : s[3] = rho * w;
145 : 268 : s[4] = rho * (spi + ke);
146 : :
147 [ - - ][ - - ]: 0 : } else Throw( "IC box initiate type not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
148 : 759 : }
149 : :
150 : : } //inciter::
|
