Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/CompFlow/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 single-material equations inside the user-defined box.
      11                 :            : */
      12                 :            : // *****************************************************************************
      13                 :            : #ifndef BoxInitialization_h
      14                 :            : #define BoxInitialization_h
      15                 :            : 
      16                 :            : #include "Fields.hpp"
      17                 :            : #include "EoS/EOS.hpp"
      18                 :            : #include "ContainerUtil.hpp"
      19                 :            : 
      20                 :            : namespace inciter {
      21                 :            : 
      22                 :            : //! Set the solution in the user-defined IC box
      23                 :            : 
      24                 :            : template< class B >
      25                 :       2528 : void initializeBox( const std::vector< EOS >& mat_blk,
      26                 :            :                     tk::real VRatio,
      27                 :            :                     tk::real V_ex,
      28                 :            :                     tk::real t,
      29                 :            :                     const B& b,
      30                 :            :                     tk::real bgpreic,
      31                 :            :                     tk::real cv,
      32                 :            :                     std::vector< tk::real >& s )
      33                 :            : // *****************************************************************************
      34                 :            : // Set the solution in the user-defined IC box/block
      35                 :            : //! \tparam B IC-block type to operate, ctr::box, or ctr::meshblock
      36                 :            : //! \param[in] VRatio Ratio of exact box volume to discrete box volume
      37                 :            : //! \param[in] V_ex Exact box volume
      38                 :            : //! \param[in] t Physical time
      39                 :            : //! \param[in] b IC box configuration to use
      40                 :            : //! \param[in] bgpreic Background pressure user input
      41                 :            : //! \param[in] cv Specific heats ratio user input
      42                 :            : //! \param[in,out] s Solution vector that is set to box ICs
      43                 :            : //! \details This function sets the fluid density and total specific energy
      44                 :            : //!   within a box initial condition, configured by the user. If the user
      45                 :            : //!   is specified a box where mass is specified, we also assume here that
      46                 :            : //!   internal energy content (energy per unit volume) is also
      47                 :            : //!   specified. Specific internal energy (energy per unit mass) is then
      48                 :            : //!   computed here (and added to the kinetic energy) from the internal
      49                 :            : //!   energy per unit volume by multiplying it with the total box volume
      50                 :            : //!   and dividing it by the total mass of the material in the box.
      51                 :            : //!   Example (SI) units of the quantities involved:
      52                 :            : //!    * internal energy content (energy per unit volume): J/m^3
      53                 :            : //!    * specific energy (internal energy per unit mass): J/kg
      54                 :            : // *****************************************************************************
      55                 :            : {
      56                 :       2528 :   const auto& initiate = b.template get< tag::initiate >();
      57                 :            : 
      58                 :       2528 :   auto boxrho = b.template get< tag::density >();
      59                 :       2528 :   const auto& boxvel = b.template get< tag::velocity >();
      60                 :       2528 :   auto boxpre = b.template get< tag::pressure >();
      61                 :       2528 :   auto boxene = b.template get< tag::energy >();
      62                 :       2528 :   auto boxtem = b.template get< tag::temperature >();
      63                 :       2528 :   auto boxmas = b.template get< tag::mass >();
      64                 :       2528 :   auto boxenc = b.template get< tag::energy_content >();
      65                 :            : 
      66                 :       2528 :   tk::real rho = 0.0, ru = 0.0, rv = 0.0, rw = 0.0, re = 0.0, spi = 0.0;
      67                 :       2528 :   bool boxmassic = false;
      68         [ +  + ]:       2528 :   if (boxmas > 0.0) {
      69 [ -  + ][ -  - ]:        268 :     Assert( boxenc > 0.0, "Box energy content must be nonzero" );
         [ -  - ][ -  - ]
      70                 :        268 :     rho = boxmas / V_ex;
      71                 :        268 :     spi = boxenc * VRatio / rho;
      72                 :        268 :     boxmassic = true;
      73                 :            :   } else {
      74         [ +  - ]:       2260 :     if (boxrho > 0.0) rho = boxrho;
      75         [ +  - ]:       2260 :     if (boxvel.size() == 3) {
      76                 :       2260 :       ru = rho * boxvel[0];
      77                 :       2260 :       rv = rho * boxvel[1];
      78                 :       2260 :       rw = rho * boxvel[2];
      79                 :            :     }
      80         [ +  - ]:       2260 :     if (boxpre > 0.0) {
      81         [ +  - ]:       2260 :       re = mat_blk[0].compute< EOS::totalenergy >( rho, ru/rho, rv/rho, rw/rho,
      82                 :            :         boxpre );
      83                 :            :     }
      84         [ -  + ]:       2260 :     if (boxene > 0.0) {
      85                 :          0 :       auto ux = ru/rho, uy = rv/rho, uz = rw/rho;
      86                 :          0 :       auto ke = 0.5*(ux*ux + uy*uy + uz*uz);
      87                 :          0 :       re = rho * (boxene + ke);
      88                 :            :     }
      89         [ -  + ]:       2260 :     if (boxtem > 0.0) re = rho * boxtem * cv;
      90                 :            :   }
      91                 :            : 
      92                 :            :   // Initiate type 'impulse' simply assigns the prescribed values to all
      93                 :            :   // nodes within a box.
      94         [ +  + ]:       2528 :   if (initiate == ctr::InitiateType::IMPULSE) {
      95                 :            :     // superimpose on existing velocity field
      96                 :       2260 :     auto u = s[1]/s[0], v = s[2]/s[0], w = s[3]/s[0];
      97                 :       2260 :     auto ke = 0.5*(u*u + v*v + w*w);
      98                 :       2260 :     s[0] = rho;
      99         [ -  + ]:       2260 :     if (boxmassic) {
     100                 :          0 :       s[1] = rho * u;
     101                 :          0 :       s[2] = rho * v;
     102                 :          0 :       s[3] = rho * w;
     103                 :          0 :       s[4] = rho * (spi + ke);
     104                 :            :     } else {
     105                 :       2260 :       s[1] = ru;
     106                 :       2260 :       s[2] = rv;
     107                 :       2260 :       s[3] = rw;
     108                 :       2260 :       s[4] = re;
     109                 :            :     }
     110                 :            :   }
     111                 :            :   // Initiate type 'linear' assigns the prescribed values to all nodes
     112                 :            :   // within a box. This is followed by adding a time-dependent energy
     113                 :            :   // source term representing a planar wave-front propagating along the
     114                 :            :   // z-direction with a velocity specified in the IC linear...end block.
     115                 :            :   // The wave front is smoothed out to include a couple of mesh nodes.
     116                 :            :   // see boxSrc() for further details.
     117 [ +  - ][ +  - ]:        268 :   else if (initiate == ctr::InitiateType::LINEAR && t < 1e-12) {
     118                 :            : 
     119                 :            :     // superimpose on existing velocity field
     120                 :        268 :     const auto u = s[1]/s[0],
     121                 :        268 :                v = s[2]/s[0],
     122                 :        268 :                w = s[3]/s[0];
     123                 :        268 :     const auto ke = 0.5*(u*u + v*v + w*w);
     124                 :            : 
     125                 :            :     // The linear-propagating source initialization can be done only
     126                 :            :     // based on background pressure (not on temperature): The IC box can
     127                 :            :     // have a different density than the background, while having the
     128                 :            :     // same pressure and temperature as the background. This means, the
     129                 :            :     // material in the box has a different specific heat (Cv) than the
     130                 :            :     // background material. If such a box has to be initialized based on
     131                 :            :     // temperature, the Cv of the box will have to be specified
     132                 :            :     // separately. This is not currently supported.
     133         [ +  - ]:        268 :     if (bgpreic > 0.0) {
     134                 :            :       // energy based on box density and background pressure
     135         [ +  - ]:        268 :       spi = mat_blk[0].compute< EOS::totalenergy >(rho, u, v, w, bgpreic) / rho;
     136 [ -  - ][ -  - ]:          0 :     } else Throw( "Background pressure must be specified for box-IC "
                 [ -  - ]
     137                 :            :                   "with linear propagating source");
     138                 :            : 
     139                 :        268 :     s[0] = rho;
     140                 :        268 :     s[1] = rho * u;
     141                 :        268 :     s[2] = rho * v;
     142                 :        268 :     s[3] = rho * w;
     143                 :        268 :     s[4] = rho * (spi + ke);
     144                 :            : 
     145 [ -  - ][ -  - ]:          0 :   } else Throw( "IC box initiate type not implemented" );
                 [ -  - ]
     146                 :       2528 : }
     147                 :            : 
     148                 :            : } //inciter::
     149                 :            : 
     150                 :            : #endif // BoxInitialization_h