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 : 1150 : 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 : 1150 : auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
56 : :
57 : : const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
58 : :
59 : : const auto& initiate = b.template get< tag::initiate >();
60 : :
61 : : // get material id in box (offset by 1, since input deck uses 1-based ids)
62 : 1150 : std::size_t boxmatid = b.template get< tag::materialid >() - 1;
63 : : const auto& boxvel = b.template get< tag::velocity >();
64 : 1150 : auto boxpre = b.template get< tag::pressure >();
65 : 1150 : auto boxene = b.template get< tag::energy >();
66 : 1150 : auto boxtemp = b.template get< tag::temperature >();
67 : 1150 : auto boxmas = b.template get< tag::mass >();
68 : 1150 : auto boxenc = b.template get< tag::energy_content >();
69 : :
70 : : auto alphamin = 1.0e-12;
71 : :
72 : : // [I] Compute the states inside the IC box/block based on the type of user
73 : : // input.
74 : :
75 : : // material volume fractions
76 [ + + ]: 3450 : for (std::size_t k=0; k<nmat; ++k) {
77 [ + + ]: 2300 : if (k == boxmatid) {
78 : 1150 : s[volfracIdx(nmat,k)] = 1.0 - (static_cast< tk::real >(nmat-1))*alphamin;
79 : : }
80 : : else {
81 : 1150 : s[volfracIdx(nmat,k)] = alphamin;
82 : : }
83 : : }
84 : : // material states (density, pressure, velocity)
85 : 1150 : tk::real u = 0.0, v = 0.0, w = 0.0, spi(0.0), pr(0.0), tmp(0.0), rbulk(0.0);
86 : 1150 : std::vector< tk::real > rhok(nmat, 0.0);
87 : : // 1. User-specified mass, specific energy (J/m^3) and volume of box
88 [ + + ]: 1150 : if (boxmas > 0.0) {
89 [ - + ][ - - ]: 104 : if (boxenc <= 1e-12) Throw( "Box energy content must be nonzero" );
[ - - ][ - - ]
[ - - ][ - - ]
[ - - ]
90 : : // determine density and energy of material in the box
91 [ - + ]: 104 : rhok[boxmatid] = boxmas / V_ex;
92 : 104 : spi = boxenc / rhok[boxmatid];
93 : :
94 : : // Determine pressure and temperature
95 : 104 : auto boxmat_vf = s[volfracIdx(nmat,boxmatid)];
96 : :
97 : : // Since initiate type 'linear' assigns the background IC values to all
98 : : // nodes within a box at initialization (followed by adding a time-dependent
99 : : // energy source term representing a propagating wave-front), the pressure
100 : : // in the box needs to be set to background pressure.
101 [ - + ][ - - ]: 104 : if (initiate == ctr::InitiateType::LINEAR && t < 1e-12) {
102 [ - - ][ - - ]: 0 : if (boxmas <= 1e-12 || boxenc <= 1e-12 || bgpreic <= 1e-12 ||
[ - - ]
103 : : bgtempic <= 1e-12)
104 [ - - ][ - - ]: 0 : Throw("Box mass, energy content, background pressure and background "
[ - - ][ - - ]
[ - - ][ - - ]
105 : : "temperature must be specified for IC with linear propagating source");
106 : :
107 : 0 : pr = bgpreic;
108 [ - - ]: 0 : auto te = mat_blk[boxmatid].compute< EOS::totalenergy >(
109 : : rhok[boxmatid], u, v, w, pr);
110 [ - - ]: 0 : tmp = mat_blk[boxmatid].compute< EOS::temperature >(
111 [ - - ]: 0 : boxmat_vf*rhok[boxmatid], u, v, w, boxmat_vf*te, boxmat_vf );
112 : :
113 : : // if the above density and pressure lead to an invalid thermodynamic
114 : : // state (negative temperature/energy), change temperature to background
115 : : // temperature and use corresponding density.
116 [ - - ][ - - ]: 0 : if (tmp < 0.0 || te < 0.0) {
117 : 0 : tmp = bgtempic;
118 [ - - ]: 0 : rhok[boxmatid] = mat_blk[boxmatid].compute< EOS::density >(pr, tmp);
119 : 0 : spi = boxenc / rhok[boxmatid];
120 : : }
121 : : }
122 : : // For initiate type 'impulse', pressure and temperature are determined from
123 : : // energy content that needs to be dumped into the box at IC.
124 [ + - ]: 104 : else if (initiate == ctr::InitiateType::IMPULSE) {
125 [ + - ]: 104 : pr = mat_blk[boxmatid].compute< EOS::pressure >(
126 [ + - ]: 104 : boxmat_vf*rhok[boxmatid], u, v, w, boxmat_vf*rhok[boxmatid]*spi,
127 : : boxmat_vf, boxmatid );
128 [ + - ]: 104 : tmp = mat_blk[boxmatid].compute< EOS::temperature >(
129 [ + - ][ - - ]: 208 : boxmat_vf*rhok[boxmatid], u, v, w, boxmat_vf*rhok[boxmatid]*spi,
130 : : boxmat_vf );
131 : : }
132 [ - - ][ - - ]: 0 : else Throw( "IC box initiate type not implemented for multimat" );
[ - - ][ - - ]
[ - - ][ - - ]
133 : :
134 : : // find density of trace material quantities in the box based on pressure
135 [ + + ]: 312 : for (std::size_t k=0; k<nmat; ++k) {
136 [ + + ]: 208 : if (k != boxmatid) {
137 [ + - ]: 208 : rhok[k] = mat_blk[k].compute< EOS::density >(pr, tmp);
138 : : }
139 : : }
140 : : }
141 : : // 2. User-specified temperature, pressure and velocity in box
142 : : else {
143 [ + + ]: 3138 : for (std::size_t k=0; k<nmat; ++k) {
144 [ + - ]: 4184 : rhok[k] = mat_blk[k].compute< EOS::density >(boxpre, boxtemp);
145 : : }
146 [ + - ]: 1046 : if (boxvel.size() == 3) {
147 : 1046 : u = boxvel[0];
148 : 1046 : v = boxvel[1];
149 : 1046 : w = boxvel[2];
150 : : }
151 [ + - ]: 1046 : if (boxpre > 0.0) {
152 : 1046 : pr = boxpre;
153 : : }
154 [ - + ]: 1046 : if (boxene > 0.0) {
155 [ - - ][ - - ]: 0 : Throw("IC-box with specified energy not set up for multimat");
[ - - ][ - - ]
[ - - ][ - - ]
156 : : }
157 : : }
158 : : // bulk density
159 [ + + ]: 3450 : for (std::size_t k=0; k<nmat; ++k) rbulk += s[volfracIdx(nmat,k)]*rhok[k];
160 : :
161 : : // [II] Finally initialize the solution vector
162 [ + + ]: 3450 : for (std::size_t k=0; k<nmat; ++k) {
163 : : // partial density
164 [ + + ]: 2300 : s[densityIdx(nmat,k)] = s[volfracIdx(nmat,k)] * rhok[k];
165 : : // total specific energy
166 [ + + ][ + + ]: 2300 : if (boxmas > 0.0 && k == boxmatid &&
167 [ + - ]: 104 : initiate == ctr::InitiateType::IMPULSE) {
168 : 104 : s[energyIdx(nmat,k)] = s[volfracIdx(nmat,k)] * rhok[k] * spi;
169 : : }
170 : : else {
171 : : // TEMP: Eventually we would need to initialize gk from control file
172 : : std::array< std::array< tk::real, 3 >, 3 > gk;
173 [ + + ]: 2196 : if (solidx[k] > 0) {
174 [ + + ]: 5976 : for (std::size_t i=0; i<3; ++i) {
175 [ + + ]: 17928 : for (std::size_t j=0; j<3; ++j) {
176 [ + + ]: 13446 : if (i==j) gk[i][j] = 1.0;
177 : 8964 : else gk[i][j] = 0.0;
178 : 13446 : s[deformIdx(nmat,solidx[k],i,j)] = gk[i][j];
179 : : }
180 : : }
181 : : }
182 : : else {
183 : 702 : gk = {{}};
184 : : }
185 : 2196 : s[energyIdx(nmat,k)] = s[volfracIdx(nmat,k)] *
186 [ + - ]: 2196 : mat_blk[k].compute< EOS::totalenergy >( rhok[k], u, v, w, pr, gk );
187 : : }
188 : : }
189 : : // bulk momentum
190 [ + - ]: 1150 : s[momentumIdx(nmat,0)] = rbulk * u;
191 : 1150 : s[momentumIdx(nmat,1)] = rbulk * v;
192 [ + - ]: 1150 : s[momentumIdx(nmat,2)] = rbulk * w;
193 : 1150 : }
194 : :
195 : : } //inciter::
196 : :
197 : : #endif // MultiMatBoxInitialization_h
|
