Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/MultiMat/Physics/FVEnergyPill.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     Physics policy for the Euler equation governing multi-material flow
       9                 :            :     using a finite volume method
      10                 :            :   \details   This file defines a Physics policy class for the compressible
      11                 :            :     flow equations class fv::MultiMat, defined in PDE/MultiMat/FVMultiMat.h.
      12                 :            :     This specific class allows energy pill initialization of a user defined
      13                 :            :     box/block for multi-material flow. See PDE/MultiMat/Physics/FV.h for general
      14                 :            :     requirements on Physics policy classes for fv::MultiMat.
      15                 :            : */
      16                 :            : // *****************************************************************************
      17                 :            : 
      18                 :            : #include "FVEnergyPill.hpp"
      19                 :            : #include "Inciter/InputDeck/InputDeck.hpp"
      20                 :            : #include "ContainerUtil.hpp"
      21                 :            : #include "MultiMat/MultiMatIndexing.hpp"
      22                 :            : #include "NoWarning/charm++.hpp"
      23                 :            : 
      24                 :            : namespace inciter {
      25                 :            : 
      26                 :            : extern ctr::InputDeck g_inputdeck;
      27                 :            : 
      28                 :            : } // inciter::
      29                 :            : 
      30                 :            : using inciter::fv::MultiMatPhysicsEnergyPill;
      31                 :            : 
      32                 :            : tk::real
      33                 :        200 : MultiMatPhysicsEnergyPill::dtRestriction(
      34                 :            :   const tk::Fields& geoElem,
      35                 :            :   std::size_t nelem,
      36                 :            :   const std::vector< int >& engSrcAd ) const
      37                 :            : // *****************************************************************************
      38                 :            : //  Compute the time step size restriction based on this physics
      39                 :            : //! \param[in] geoElem Element geometry array
      40                 :            : //! \param[in] nelem Number of elements
      41                 :            : //! \param[in] engSrcAd Whether the energy source was added
      42                 :            : //! \return Maximum allowable time step based on front propagation speed
      43                 :            : // *****************************************************************************
      44                 :            : {
      45                 :        200 :   auto mindt = std::numeric_limits< tk::real >::max();
      46                 :            :   // determine front propagation speed
      47                 :            :   const auto& icmbk = g_inputdeck.get< tag::ic, tag::meshblock >();
      48                 :        200 :   tk::real v_front(0.0);
      49         [ +  + ]:        400 :   for (const auto& b : icmbk) { // for all blocks
      50                 :        200 :     auto inittype = b.template get< tag::initiate >();
      51         [ -  + ]:        200 :     if (inittype == ctr::InitiateType::LINEAR) {
      52                 :          0 :       v_front = std::max(v_front,
      53                 :            :         b.template get< tag::front_speed >());
      54                 :            :     }
      55                 :            :   }
      56                 :            : 
      57         [ +  + ]:      46400 :   for (std::size_t e=0; e<nelem; ++e)
      58                 :            :   {
      59                 :            :     // characteristic length (radius of insphere)
      60         [ -  + ]:      46200 :     auto dx = std::min(std::cbrt(geoElem(e,0)), geoElem(e,4))
      61                 :      46200 :       /std::sqrt(24.0);
      62                 :            : 
      63                 :            :     // element dt restriction if relevant energy sources were added
      64 [ -  + ][ -  - ]:      46200 :     if (engSrcAd[e] > 0 && std::abs(v_front) > 1e-8)
      65         [ -  - ]:          0 :       mindt = std::min(mindt, dx/v_front);
      66                 :            :   }
      67                 :            : 
      68                 :        200 :   return mindt;
      69                 :            : }
      70                 :            : 
      71                 :        400 : void MultiMatPhysicsEnergyPill::
      72                 :            : physSrc(
      73                 :            :   std::size_t nmat,
      74                 :            :   tk::real t,
      75                 :            :   const tk::Fields& geoElem,
      76                 :            :   const std::unordered_map< std::size_t, std::set< std::size_t > >& elemblkid,
      77                 :            :   tk::Fields& R,
      78                 :            :   std::vector< int >& engSrcAdded ) const
      79                 :            : // *****************************************************************************
      80                 :            : //! Compute sources corresponding to a propagating front in user-defined box
      81                 :            : //! \param[in] nmat Number of materials
      82                 :            : //! \param[in] t Physical time
      83                 :            : //! \param[in] geoElem Element geometry array
      84                 :            : //! \param[in] elemblkid Element ids associated with mesh block ids where
      85                 :            : //!   user ICs are set
      86                 :            : //! \param[in,out] R Right-hand side vector
      87                 :            : //! \param[in,out] engSrcAdded Whether the energy source was added
      88                 :            : //! \details This function adds the energy source corresponding to a
      89                 :            : //!   spherically growing wave-front propagating with a user-specified
      90                 :            : //!   velocity, within a user-configured mesh block initial condition.
      91                 :            : //!   Example (SI) units of the quantities involved:
      92                 :            : //!    * internal energy content (energy per unit volume): J/m^3
      93                 :            : //!    * specific energy (internal energy per unit mass): J/kg
      94                 :            : // *****************************************************************************
      95                 :            : {
      96                 :            :   const auto& icmbk = g_inputdeck.get< tag::ic, tag::meshblock >();
      97         [ +  + ]:        800 :   for (const auto& mb : icmbk) { // for all blocks
      98                 :        400 :     auto blid = mb.get< tag::blockid >();
      99         [ +  - ]:        400 :     if (elemblkid.find(blid) != elemblkid.end()) { // if elements exist in blk
     100                 :            :       const auto& initiate = mb.template get< tag::initiate >();
     101         [ -  + ]:        400 :       if (initiate == ctr::InitiateType::LINEAR) { // if propagating src
     102                 :            : 
     103                 :            :         const auto& blkelems = tk::cref_find(elemblkid,blid);
     104                 :            : 
     105                 :          0 :         auto enc = mb.template get< tag::energy_content >();
     106                 :            :         Assert( enc > 0.0, "Box energy content must be nonzero" );
     107                 :            :         const auto& x0_front = mb.template get< tag::point >();
     108                 :            :         Assert(x0_front.size()==3, "Incorrectly sized front initial location");
     109                 :          0 :         auto blkmatid = mb.template get< tag::materialid >();
     110                 :            : 
     111                 :            :         // determine times at which sourcing is initialized and terminated
     112                 :          0 :         auto v_front = mb.template get< tag::front_speed >();
     113                 :          0 :         auto w_front = mb.template get< tag::front_width >();
     114                 :          0 :         auto tInit = mb.template get< tag::init_time >();
     115                 :            : 
     116         [ -  - ]:          0 :         if (t >= tInit) {
     117                 :            :           // current radius of front
     118                 :          0 :           tk::real r_front = v_front * (t-tInit);
     119                 :            :           // arbitrary shape form
     120                 :          0 :           auto amplE = enc * v_front / w_front;
     121                 :            : 
     122 [ -  - ][ -  - ]:          0 :           for (auto e : blkelems) {
     123                 :          0 :             std::array< tk::real, 3 > node{{ geoElem(e,1), geoElem(e,2),
     124                 :          0 :               geoElem(e,3) }};
     125                 :            : 
     126                 :          0 :             auto r_e = std::sqrt(
     127                 :          0 :               (node[0]-x0_front[0])*(node[0]-x0_front[0]) +
     128                 :          0 :               (node[1]-x0_front[1])*(node[1]-x0_front[1]) +
     129                 :          0 :               (node[2]-x0_front[2])*(node[2]-x0_front[2]) );
     130                 :            : 
     131                 :            :             // if element centroid lies within spherical shell add sources
     132 [ -  - ][ -  - ]:          0 :             if (r_e >= r_front && r_e <= r_front+w_front) {
     133                 :            :               // Add the source term to the rhs
     134                 :          0 :               R(e, energyDofIdx(nmat,blkmatid-1,1,0)) += geoElem(e,0) * amplE;
     135                 :          0 :               engSrcAdded[e] = 1;
     136                 :            :             }
     137                 :            :           }
     138                 :            :         }
     139                 :            : 
     140                 :            :       }
     141                 :            :     }
     142                 :            :   }
     143                 :        400 : }