Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiMat/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 multi-material equations inside the user-defined box.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef MultiMatBoxInitialization_h
14 : : #define MultiMatBoxInitialization_h
15 : :
16 : : #include "Fields.hpp"
17 : : #include "EoS/EOS.hpp"
18 : : #include "ContainerUtil.hpp"
19 : : #include "MultiMat/MultiMatIndexing.hpp"
20 : :
21 : : namespace inciter {
22 : :
23 : : using ncomp_t = tk::ncomp_t;
24 : :
25 : : template< class B >
26 : 1964 : void initializeBox( const std::vector< EOS >& mat_blk,
27 : : tk::real V_ex,
28 : : tk::real t,
29 : : const B& b,
30 : : tk::real bgpreic,
31 : : tk::real bgtempic,
32 : : std::vector< tk::real >& s )
33 : : // *****************************************************************************
34 : : // Set the solution in the user-defined IC box/mesh block
35 : : //! \tparam B IC-block type to operate, ctr::box, or ctr::meshblock
36 : : //! \param[in] V_ex Exact box volume
37 : : //! \param[in] t Physical time
38 : : //! \param[in] b IC box configuration to use
39 : : //! \param[in] bgpreic Background pressure user input
40 : : //! \param[in] bgtempic Background temperature user input
41 : : //! \param[in,out] s Solution vector that is set to box ICs
42 : : //! \details This function sets the fluid density and total specific energy
43 : : //! within a box initial condition, configured by the user. If the user
44 : : //! is specified a box where mass is specified, we also assume here that
45 : : //! internal energy content (energy per unit volume) is also
46 : : //! specified. Specific internal energy (energy per unit mass) is then
47 : : //! computed here (and added to the kinetic energy) from the internal
48 : : //! energy per unit volume by multiplying it with the total box volume
49 : : //! and dividing it by the total mass of the material in the box.
50 : : //! Example (SI) units of the quantities involved:
51 : : //! * internal energy content (energy per unit volume): J/m^3
52 : : //! * specific energy (internal energy per unit mass): J/kg
53 : : // *****************************************************************************
54 : : {
55 : 1964 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
56 : 1964 : auto alphamin = g_inputdeck.get< tag::multimat, tag::min_volumefrac >();
57 : :
58 : 1964 : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
59 : :
60 : 1964 : const auto& initiate = b.template get< tag::initiate >();
61 : :
62 : : // get material id in box (offset by 1, since input deck uses 1-based ids)
63 : 1964 : std::size_t boxmatid = b.template get< tag::materialid >() - 1;
64 : 1964 : const auto& boxvel = b.template get< tag::velocity >();
65 : 1964 : auto boxpre = b.template get< tag::pressure >();
66 : 1964 : auto boxene = b.template get< tag::energy >();
67 : 1964 : auto boxtemp = b.template get< tag::temperature >();
68 : 1964 : auto boxmas = b.template get< tag::mass >();
69 : 1964 : auto boxenc = b.template get< tag::energy_content >();
70 : :
71 : : // [I] Compute the states inside the IC box/block based on the type of user
72 : : // input.
73 : :
74 : : // material volume fractions
75 [ + + ]: 5892 : for (std::size_t k=0; k<nmat; ++k) {
76 [ + + ]: 3928 : if (k == boxmatid) {
77 : 1964 : s[volfracIdx(nmat,k)] = 1.0 - (static_cast< tk::real >(nmat-1))*alphamin;
78 : : }
79 : : else {
80 : 1964 : s[volfracIdx(nmat,k)] = alphamin;
81 : : }
82 : : }
83 : : // material states (density, pressure, velocity)
84 : 1964 : tk::real u = 0.0, v = 0.0, w = 0.0, spi(0.0), pr(0.0), tmp(0.0), rbulk(0.0);
85 [ + - ]: 1964 : std::vector< tk::real > rhok(nmat, 0.0);
86 : : // 1. User-specified mass, specific energy (J/m^3) and volume of box
87 [ - + ]: 1964 : if (boxmas > 0.0) {
88 [ - - ][ - - ]: 0 : if (boxenc <= 1e-12) Throw( "Box energy content must be nonzero" );
[ - - ][ - - ]
89 : : // determine density and energy of material in the box
90 : 0 : rhok[boxmatid] = boxmas / V_ex;
91 : 0 : spi = boxenc / rhok[boxmatid];
92 : :
93 : : // Determine pressure and temperature
94 : 0 : auto boxmat_vf = s[volfracIdx(nmat,boxmatid)];
95 : :
96 : : // Since initiate type 'linear' assigns the background IC values to all
97 : : // nodes within a box at initialization (followed by adding a time-dependent
98 : : // energy source term representing a propagating wave-front), the pressure
99 : : // in the box needs to be set to background pressure.
100 [ - - ][ - - ]: 0 : if (initiate == ctr::InitiateType::LINEAR && t < 1e-12) {
101 [ - - ][ - - ]: 0 : if (boxmas <= 1e-12 || boxenc <= 1e-12 || bgpreic <= 1e-12 ||
[ - - ][ - - ]
102 : : bgtempic <= 1e-12)
103 [ - - ][ - - ]: 0 : Throw("Box mass, energy content, background pressure and background "
[ - - ]
104 : : "temperature must be specified for IC with linear propagating source");
105 : :
106 : 0 : pr = bgpreic;
107 : 0 : auto te = mat_blk[boxmatid].compute< EOS::totalenergy >(
108 [ - - ]: 0 : boxmat_vf*rhok[boxmatid], u, v, w, boxmat_vf*pr, boxmat_vf );
109 : 0 : tmp = mat_blk[boxmatid].compute< EOS::temperature >(
110 [ - - ]: 0 : boxmat_vf*rhok[boxmatid], u, v, w, te, boxmat_vf );
111 : :
112 : : // if the above density and pressure lead to an invalid thermodynamic
113 : : // state (negative temperature/energy), change temperature to background
114 : : // temperature and use corresponding density.
115 [ - - ][ - - ]: 0 : if (tmp < 0.0 || te < 0.0) {
116 : 0 : tmp = bgtempic;
117 [ - - ]: 0 : rhok[boxmatid] = mat_blk[boxmatid].compute< EOS::density >(pr, tmp);
118 : 0 : spi = boxenc / rhok[boxmatid];
119 : 0 : }
120 : : }
121 : : // For initiate type 'impulse', pressure and temperature are determined from
122 : : // energy content that needs to be dumped into the box at IC.
123 [ - - ]: 0 : else if (initiate == ctr::InitiateType::IMPULSE) {
124 : 0 : pr = mat_blk[boxmatid].compute< EOS::pressure >(
125 [ - - ]: 0 : boxmat_vf*rhok[boxmatid], u, v, w, boxmat_vf*rhok[boxmatid]*spi,
126 : : boxmat_vf, boxmatid );
127 : 0 : tmp = mat_blk[boxmatid].compute< EOS::temperature >(
128 [ - - ]: 0 : boxmat_vf*rhok[boxmatid], u, v, w, boxmat_vf*rhok[boxmatid]*spi,
129 : : boxmat_vf );
130 : : }
131 [ - - ][ - - ]: 0 : else Throw( "IC box initiate type not implemented for multimat" );
[ - - ]
132 : :
133 : : // find density of trace material quantities in the box based on pressure
134 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
135 [ - - ]: 0 : if (k != boxmatid) {
136 [ - - ]: 0 : rhok[k] = mat_blk[k].compute< EOS::density >(pr, tmp);
137 : : }
138 : : }
139 : : }
140 : : // 2. User-specified temperature, pressure and velocity in box
141 : : else {
142 [ + + ]: 5892 : for (std::size_t k=0; k<nmat; ++k) {
143 [ + - ]: 3928 : rhok[k] = mat_blk[k].compute< EOS::density >(boxpre, boxtemp);
144 : : }
145 [ + - ]: 1964 : if (boxvel.size() == 3) {
146 : 1964 : u = boxvel[0];
147 : 1964 : v = boxvel[1];
148 : 1964 : w = boxvel[2];
149 : : }
150 [ + - ]: 1964 : if (boxpre > 0.0) {
151 : 1964 : pr = boxpre;
152 : : }
153 [ - + ]: 1964 : if (boxene > 0.0) {
154 [ - - ][ - - ]: 0 : Throw("IC-box with specified energy not set up for multimat");
[ - - ]
155 : : }
156 : : }
157 : : // bulk density
158 [ + + ]: 5892 : for (std::size_t k=0; k<nmat; ++k) rbulk += s[volfracIdx(nmat,k)]*rhok[k];
159 : :
160 : : // [II] Finally initialize the solution vector
161 [ + + ]: 5892 : for (std::size_t k=0; k<nmat; ++k) {
162 : : // partial density
163 : 3928 : s[densityIdx(nmat,k)] = s[volfracIdx(nmat,k)] * rhok[k];
164 : : // total specific energy
165 [ - + ][ - - ]: 3928 : if (boxmas > 0.0 && k == boxmatid &&
166 [ - - ]: 0 : initiate == ctr::InitiateType::IMPULSE) {
167 : 0 : s[energyIdx(nmat,k)] = s[volfracIdx(nmat,k)] * rhok[k] * spi;
168 : : }
169 : : else {
170 : : // TEMP: Eventually we would need to initialize gk from control file
171 : : std::array< std::array< tk::real, 3 >, 3 > gk;
172 [ - + ]: 3928 : if (solidx[k] > 0) {
173 [ - - ]: 0 : for (std::size_t i=0; i<3; ++i) {
174 [ - - ]: 0 : for (std::size_t j=0; j<3; ++j) {
175 [ - - ]: 0 : if (i==j) gk[i][j] = 1.0;
176 : 0 : else gk[i][j] = 0.0;
177 : 0 : s[deformIdx(nmat,solidx[k],i,j)] = gk[i][j];
178 : : }
179 : : }
180 : : }
181 : : else {
182 : 3928 : gk = {{}};
183 : : }
184 : 3928 : s[energyIdx(nmat,k)] =
185 [ + - ]: 7856 : mat_blk[k].compute< EOS::totalenergy >( s[volfracIdx(nmat,k)]*rhok[k],
186 : 7856 : u, v, w, s[volfracIdx(nmat,k)]*pr, s[volfracIdx(nmat,k)], gk );
187 : : }
188 : : }
189 : : // bulk momentum
190 : 1964 : s[momentumIdx(nmat,0)] = rbulk * u;
191 : 1964 : s[momentumIdx(nmat,1)] = rbulk * v;
192 : 1964 : s[momentumIdx(nmat,2)] = rbulk * w;
193 : 1964 : }
194 : :
195 : : } //inciter::
196 : :
197 : : #endif // MultiMatBoxInitialization_h
|
