Branch data     Line data    Source code

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