Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiMat/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 "Control/Inciter/Types.hpp"
15 : : #include "Inciter/InputDeck/InputDeck.hpp"
16 : : #include "BoxInitialization.hpp"
17 : : #include "EoS/EoS.hpp"
18 : : #include "ContainerUtil.hpp"
19 : : #include "MultiMat/MultiMatIndexing.hpp"
20 : :
21 : : namespace inciter {
22 : :
23 : : extern ctr::InputDeck g_inputdeck;
24 : :
25 : 0 : void initializeBox( std::size_t system,
26 : : tk::real VRatio,
27 : : tk::real,
28 : : const inciter::ctr::box& b,
29 : : std::vector< tk::real >& s )
30 : : // *****************************************************************************
31 : : // Set the solution in the user-defined IC box
32 : : //! \param[in] system Equation system index
33 : : //! \param[in] VRatio Ratio of exact box volume to discrete box volume
34 : : //! \param[in] b IC box configuration to use
35 : : //! \param[in,out] s Solution vector that is set to box ICs
36 : : //! \details This function sets the fluid density and total specific energy
37 : : //! within a box initial condition, configured by the user. If the user
38 : : //! is specified a box where mass is specified, we also assume here that
39 : : //! internal energy content (energy per unit volume) is also
40 : : //! specified. Specific internal energy (energy per unit mass) is then
41 : : //! computed here (and added to the kinetic energy) from the internal
42 : : //! energy per unit volume by multiplying it with the total box volume
43 : : //! and dividing it by the total mass of the material in the box.
44 : : //! Example (SI) units of the quantities involved:
45 : : //! * internal energy content (energy per unit volume): J/m^3
46 : : //! * specific energy (internal energy per unit mass): J/kg
47 : : // *****************************************************************************
48 : : {
49 : : auto nmat =
50 : 0 : g_inputdeck.get< tag::param, tag::multimat, tag::nmat >()[system];
51 : :
52 : : std::vector< tk::real > box
53 : 0 : { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
54 : 0 : b.template get< tag::ymin >(), b.template get< tag::ymax >(),
55 : 0 : b.template get< tag::zmin >(), b.template get< tag::zmax >() };
56 : :
57 : : const auto& initiate = b.template get< tag::initiate >();
58 : 0 : auto inittype = initiate.template get< tag::init >();
59 : :
60 : 0 : auto boxmatid = b.template get< tag::materialid >();
61 : : const auto& boxvel = b.template get< tag::velocity >();
62 : 0 : auto boxpre = b.template get< tag::pressure >();
63 : 0 : auto boxene = b.template get< tag::energy >();
64 : 0 : auto boxtemp = b.template get< tag::temperature >();
65 : 0 : auto boxmas = b.template get< tag::mass >();
66 : 0 : auto boxenc = b.template get< tag::energy_content >();
67 : :
68 : : auto alphamin = 1.0e-12;
69 : :
70 : : // initialize box material volume fraction
71 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
72 [ - - ]: 0 : if (k == boxmatid-1) {
73 : 0 : s[volfracIdx(nmat,k)] = 1.0 - (static_cast< tk::real >(nmat-1))*alphamin;
74 : : }
75 : : else {
76 : 0 : s[volfracIdx(nmat,k)] = alphamin;
77 : : }
78 : : }
79 : :
80 : : // initialize box material states
81 : : tk::real u = 0.0, v = 0.0, w = 0.0;
82 [ - - ][ - - ]: 0 : std::vector< tk::real > rhok(nmat, 0.0);
83 : :
84 : : // Initiate type 'impulse' simply assigns the prescribed values to all
85 : : // nodes within a box.
86 [ - - ]: 0 : if (inittype == ctr::InitiateType::IMPULSE) {
87 [ - - ]: 0 : if (boxmas > 0.0) {
88 : : Assert( boxenc > 0.0, "Box energy content must be nonzero" );
89 : : // determine density and energy of material in the box
90 : 0 : auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
91 [ - - ]: 0 : rhok[boxmatid-1] = boxmas / V_ex;
92 : 0 : auto spi = boxenc * VRatio / rhok[boxmatid-1];
93 : :
94 : : // based on the density and energy of the material, determine pressure
95 : : // and temperature
96 : 0 : auto boxmat_vf = s[volfracIdx(nmat,boxmatid-1)];
97 : 0 : auto pr_box = eos_pressure< tag::multimat >(system,
98 [ - - ]: 0 : boxmat_vf*rhok[boxmatid-1], u, v, w, boxmat_vf*rhok[boxmatid-1]*spi,
99 : : boxmat_vf, boxmatid-1);
100 : 0 : auto t_box = eos_temperature< tag::multimat >(system,
101 : 0 : boxmat_vf*rhok[boxmatid-1], u, v, w, boxmat_vf*rhok[boxmatid-1]*spi,
102 : : boxmat_vf, boxmatid-1);
103 : :
104 : : // find density of trace material quantities in the box based on pressure
105 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
106 [ - - ]: 0 : if (k != boxmatid-1) {
107 : 0 : rhok[k] = eos_density< tag::multimat >(system, pr_box, t_box, k);
108 : : }
109 : : }
110 : :
111 : : // initialize box based on above
112 : : auto rb(0.0);
113 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
114 : : // partial density
115 [ - - ]: 0 : s[densityIdx(nmat,k)] = s[volfracIdx(nmat,k)] * rhok[k];
116 : : // total specific energy
117 [ - - ]: 0 : if (k == boxmatid-1) {
118 : 0 : s[energyIdx(nmat,k)] = s[volfracIdx(nmat,k)] * rhok[k] * spi;
119 : : }
120 : : else {
121 : 0 : s[energyIdx(nmat,k)] = s[volfracIdx(nmat,k)] *
122 : 0 : eos_totalenergy< tag::multimat >(system, rhok[k], u, v, w, pr_box,
123 : : k);
124 : : }
125 : : // bulk density
126 : 0 : rb += s[densityIdx(nmat,k)];
127 : : }
128 : : // bulk momentum
129 : 0 : s[momentumIdx(nmat,0)] = rb * u;
130 : 0 : s[momentumIdx(nmat,1)] = rb * v;
131 : 0 : s[momentumIdx(nmat,2)] = rb * w;
132 : : } else {
133 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
134 : 0 : rhok[k] = eos_density< tag::multimat >(system, boxpre, boxtemp, k);
135 : : }
136 [ - - ]: 0 : if (boxvel.size() == 3) {
137 : 0 : u = boxvel[0];
138 : 0 : v = boxvel[1];
139 : 0 : w = boxvel[2];
140 : : }
141 [ - - ]: 0 : if (boxpre > 0.0) {
142 : : auto rb(0.0);
143 [ - - ]: 0 : for (std::size_t k=0; k<nmat; ++k) {
144 : : // partial density
145 : 0 : s[densityIdx(nmat,k)] = s[volfracIdx(nmat,k)] * rhok[k];
146 : : // total specific energy
147 : 0 : s[energyIdx(nmat,k)] = s[volfracIdx(nmat,k)] *
148 : 0 : eos_totalenergy< tag::multimat >(system, rhok[k], u, v, w, boxpre,
149 : : k);
150 : : // bulk density
151 : 0 : rb += s[densityIdx(nmat,k)];
152 : : }
153 : : // bulk momentum
154 : 0 : s[momentumIdx(nmat,0)] = rb * u;
155 : 0 : s[momentumIdx(nmat,1)] = rb * v;
156 : 0 : s[momentumIdx(nmat,2)] = rb * w;
157 : : }
158 [ - - ]: 0 : if (boxene > 0.0) {
159 [ - - ][ - - ]: 0 : Throw("IC-box with specified energy not set up for multimat");
[ - - ][ - - ]
[ - - ][ - - ]
160 : : }
161 : : }
162 : : }
163 [ - - ][ - - ]: 0 : else Throw( "IC box initiate type not implemented" );
[ - - ][ - - ]
[ - - ][ - - ]
164 : 0 : }
165 : :
166 : : } //inciter::
|
