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                 :       2274 :   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                 :       2274 :   const auto& boxmassfrac = b.template get< tag::mass_fractions >();
      58                 :       2274 :   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         [ +  - ]:       4548 :   auto alphas = boxmassfrac;
      71                 :       2274 :   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