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