Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/PDE/MultiSpecies/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-species equations inside the user-defined box.
11 : : */
12 : : // *****************************************************************************
13 : : #ifndef MultiSpeciesBoxInitialization_h
14 : : #define MultiSpeciesBoxInitialization_h
15 : :
16 : : #include "Fields.hpp"
17 : : #include "EoS/EOS.hpp"
18 : : #include "ContainerUtil.hpp"
19 : : #include "MultiSpecies/MultiSpeciesIndexing.hpp"
20 : :
21 : : namespace inciter {
22 : :
23 : : using ncomp_t = tk::ncomp_t;
24 : :
25 : : template< class B >
26 : 2274 : void initializeBox( const std::vector< EOS >& mat_blk,
27 : : tk::real /*V_ex*/,
28 : : tk::real /*t*/,
29 : : const B& b,
30 : : std::vector< tk::real >& s )
31 : : // *****************************************************************************
32 : : // Set the solution in the user-defined IC box/mesh block
33 : : //! \tparam B IC-block type to operate, ctr::box, or ctr::meshblock
34 : : // //! \param[in] V_ex Exact box volume
35 : : // //! \param[in] t Physical time
36 : : //! \param[in] b IC box configuration to use
37 : : //! \param[in,out] s Solution vector that is set to box ICs
38 : : //! \details This function sets the fluid density and total specific energy
39 : : //! within a box initial condition, configured by the user. If the user
40 : : //! is specified a box where mass is specified, we also assume here that
41 : : //! internal energy content (energy per unit volume) is also
42 : : //! specified. Specific internal energy (energy per unit mass) is then
43 : : //! computed here (and added to the kinetic energy) from the internal
44 : : //! energy per unit volume by multiplying it with the total box volume
45 : : //! and dividing it by the total mass of the material in the box.
46 : : //! Example (SI) units of the quantities involved:
47 : : //! * internal energy content (energy per unit volume): J/m^3
48 : : //! * specific energy (internal energy per unit mass): J/kg
49 : : // *****************************************************************************
50 : : {
51 : :
52 : 2274 : auto nspec = g_inputdeck.get< tag::multispecies, tag::nspec >();
53 : :
54 : : const auto& initiate = b.template get< tag::initiate >();
55 : :
56 : : // get species id in box (offset by 1, since input deck uses 1-based ids)
57 : : const auto& boxmassfrac = b.template get< tag::mass_fractions >();
58 : : const auto& boxvel = b.template get< tag::velocity >();
59 : 2274 : auto boxpre = b.template get< tag::pressure >();
60 : 2274 : auto boxene = b.template get< tag::energy >();
61 : 2274 : auto boxtemp = b.template get< tag::temperature >();
62 : 2274 : auto boxmas = b.template get< tag::mass >();
63 : :
64 : 2274 : auto alphamin = 1.0e-12;
65 : :
66 : : // [I] Compute the states inside the IC box/block based on the type of user
67 : : // input.
68 : :
69 : : // species volume fractions
70 : 2274 : auto alphas = boxmassfrac;
71 : : tk::real total_al(0.0);
72 [ + + ]: 6822 : for (std::size_t k=0; k<nspec; ++k) {
73 [ + + ]: 4548 : alphas[k] = std::max(alphas[k], alphamin);
74 : 4548 : total_al += alphas[k];
75 : : }
76 [ + + ]: 6822 : for (std::size_t k=0; k<nspec; ++k) alphas[k] /= total_al;
77 : :
78 : : // material states (density, pressure, velocity)
79 : 2274 : tk::real u = 0.0, v = 0.0, w = 0.0, spi(0.0), pr(0.0), rbulk(0.0);
80 [ + - ][ - - ]: 2274 : std::vector< tk::real > rhok(nspec, 0.0);
81 : :
82 : : // 1. User-specified mass, specific energy (J/m^3) and volume of box
83 [ + - ][ + - ]: 2274 : if (boxmas > 0.0 || initiate != ctr::InitiateType::IMPULSE) {
84 [ - - ][ - - ]: 0 : Throw( "IC-box initialization type not supported in multispecies" );
[ - - ][ - - ]
[ - - ][ - - ]
85 : : }
86 : : // 2. User-specified temperature, pressure and velocity in box
87 : : else {
88 [ + + ]: 6822 : for (std::size_t k=0; k<nspec; ++k) {
89 [ + - ]: 4548 : rhok[k] = mat_blk[0].compute< EOS::density >(boxpre, boxtemp);
90 : 4548 : rbulk += alphas[k]*rhok[k];
91 : : }
92 [ + - ]: 2274 : if (boxvel.size() == 3) {
93 : 2274 : u = boxvel[0];
94 : 2274 : v = boxvel[1];
95 : 2274 : w = boxvel[2];
96 : : }
97 [ + - ]: 2274 : if (boxpre > 0.0) {
98 : 2274 : pr = boxpre;
99 : : }
100 [ - + ]: 2274 : if (boxene > 0.0) {
101 [ - - ][ - - ]: 0 : Throw("IC-box with specified energy not set up for multispecies");
[ - - ][ - - ]
[ - - ][ - - ]
102 : : }
103 : :
104 [ + - ]: 2274 : spi = mat_blk[0].compute< EOS::totalenergy >(rbulk, u, v, w, pr) / rbulk;
105 : : }
106 : :
107 : : // [II] Finally initialize the solution vector
108 : : // partial density
109 [ + + ]: 6822 : for (std::size_t k=0; k<nspec; ++k) {
110 : 4548 : s[multispecies::densityIdx(nspec,k)] = alphas[k] * rhok[k];
111 : : }
112 : : // total specific energy
113 [ + - ]: 2274 : s[multispecies::energyIdx(nspec,0)] = rbulk * spi;
114 : : // bulk momentum
115 : 2274 : s[multispecies::momentumIdx(nspec,0)] = rbulk * u;
116 [ + - ]: 2274 : s[multispecies::momentumIdx(nspec,1)] = rbulk * v;
117 [ + - ]: 2274 : s[multispecies::momentumIdx(nspec,2)] = rbulk * w;
118 : 2274 : }
119 : :
120 : : } //inciter::
121 : :
122 : : #endif // MultiSpeciesBoxInitialization_h
|
