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                 :       4548 :   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                 :       4548 :   const auto& boxmassfrac = b.template get< tag::mass_fractions >();
      59                 :       4548 :   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         [ +  - ]:       9096 :   auto Ys = boxmassfrac;
      72                 :       4548 :   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                 :       4548 :   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