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