Quinoa all test code coverage report
Current view: top level - PDE/MultiMat - FVMultiMat.hpp (source / functions) Hit Total Coverage
Commit: -128-NOTFOUND Lines: 173 281 61.6 %
Date: 2025-02-13 15:54:38 Functions: 46 500 9.2 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 126 436 28.9 %

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/MultiMat/FVMultiMat.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     Compressible multi-material flow using finite volumes
       9                 :            :   \details   This file implements calls to the physics operators governing
      10                 :            :     compressible multi-material flow (with velocity equilibrium) using finite
      11                 :            :     volume discretizations.
      12                 :            : */
      13                 :            : // *****************************************************************************
      14                 :            : #ifndef FVMultiMat_h
      15                 :            : #define FVMultiMat_h
      16                 :            : 
      17                 :            : #include <cmath>
      18                 :            : #include <algorithm>
      19                 :            : #include <unordered_set>
      20                 :            : #include <map>
      21                 :            : #include <array>
      22                 :            : 
      23                 :            : #include "Macro.hpp"
      24                 :            : #include "Exception.hpp"
      25                 :            : #include "Vector.hpp"
      26                 :            : #include "ContainerUtil.hpp"
      27                 :            : #include "UnsMesh.hpp"
      28                 :            : #include "Inciter/InputDeck/InputDeck.hpp"
      29                 :            : #include "Integrate/Basis.hpp"
      30                 :            : #include "Integrate/Quadrature.hpp"
      31                 :            : #include "Integrate/Initialize.hpp"
      32                 :            : #include "Integrate/Mass.hpp"
      33                 :            : #include "Integrate/Surface.hpp"
      34                 :            : #include "Integrate/Boundary.hpp"
      35                 :            : #include "Integrate/Volume.hpp"
      36                 :            : #include "Integrate/MultiMatTerms.hpp"
      37                 :            : #include "Integrate/Source.hpp"
      38                 :            : #include "RiemannChoice.hpp"
      39                 :            : #include "MultiMat/MultiMatIndexing.hpp"
      40                 :            : #include "Reconstruction.hpp"
      41                 :            : #include "Limiter.hpp"
      42                 :            : #include "Problem/FieldOutput.hpp"
      43                 :            : #include "Problem/BoxInitialization.hpp"
      44                 :            : #include "MultiMat/BCFunctions.hpp"
      45                 :            : #include "MultiMat/MiscMultiMatFns.hpp"
      46                 :            : 
      47                 :            : namespace inciter {
      48                 :            : 
      49                 :            : extern ctr::InputDeck g_inputdeck;
      50                 :            : 
      51                 :            : namespace fv {
      52                 :            : 
      53                 :            : //! \brief MultiMat used polymorphically with tk::FVPDE
      54                 :            : //! \details The template arguments specify policies and are used to configure
      55                 :            : //!   the behavior of the class. The policies are:
      56                 :            : //!   - Physics - physics configuration, see PDE/MultiMat/Physics.h
      57                 :            : //!   - Problem - problem configuration, see PDE/MultiMat/Problem.h
      58                 :            : //! \note The default physics is Euler, set in inciter::deck::check_multimat()
      59                 :            : template< class Physics, class Problem >
      60                 :            : class MultiMat {
      61                 :            : 
      62                 :            :   private:
      63                 :            :     using eq = tag::multimat;
      64                 :            : 
      65                 :            :   public:
      66                 :            :     //! Constructor
      67                 :         76 :     explicit MultiMat() :
      68                 :            :       m_physics(),
      69                 :         76 :       m_ncomp( g_inputdeck.get< tag::ncomp >() ),
      70                 :            :       m_riemann( multimatRiemannSolver(
      71                 :         76 :         g_inputdeck.get< tag::flux >() ) )
      72                 :            :     {
      73                 :            :       // associate boundary condition configurations with state functions
      74 [ +  - ][ +  - ]:        988 :       brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  + ]
         [ +  + ][ -  - ]
                 [ -  - ]
      75                 :            :         // BC State functions
      76                 :            :         { dirichlet
      77                 :            :         , symmetry
      78                 :            :         , invalidBC         // Outlet BC not implemented
      79                 :            :         , farfield
      80                 :            :         , extrapolate
      81                 :            :         , noslipwall },
      82                 :            :         // BC Gradient functions
      83                 :            :         { noOpGrad
      84                 :            :         , symmetryGrad
      85                 :            :         , noOpGrad
      86                 :            :         , noOpGrad
      87                 :            :         , noOpGrad
      88                 :            :         , noOpGrad }
      89                 :            :         ) );
      90                 :            : 
      91                 :            :       // Inlet BC has a different structure than above BCs, so it must be 
      92                 :            :       // handled differently than with ConfigBC
      93 [ +  - ][ +  - ]:         76 :       ConfigInletBC(m_bc, inlet, zeroGrad);
                 [ +  - ]
      94                 :            : 
      95                 :            :       // EoS initialization
      96         [ +  - ]:         76 :       initializeMaterialEoS( m_mat_blk );
      97                 :         76 :     }
      98                 :            : 
      99                 :            :     //! Find the number of primitive quantities required for this PDE system
     100                 :            :     //! \return The number of primitive quantities required to be stored for
     101                 :            :     //!   this PDE system
     102                 :      13168 :     std::size_t nprim() const
     103                 :            :     {
     104                 :      13168 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     105                 :            :       // multimat needs individual material pressures and velocities currently
     106                 :      13168 :       return (nmat+3);
     107                 :            :     }
     108                 :            : 
     109                 :            :     //! Find the number of materials set up for this PDE system
     110                 :            :     //! \return The number of materials set up for this PDE system
     111                 :          0 :     std::size_t nmat() const
     112                 :            :     {
     113                 :          0 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     114                 :          0 :       return nmat;
     115                 :            :     }
     116                 :            : 
     117                 :            :     //! Determine elements that lie inside the user-defined IC box
     118                 :            :     //! \param[in] geoElem Element geometry array
     119                 :            :     //! \param[in] nielem Number of internal elements
     120                 :            :     //! \param[in,out] inbox List of nodes at which box user ICs are set for
     121                 :            :     //!    each IC box
     122                 :        294 :     void IcBoxElems( const tk::Fields& geoElem,
     123                 :            :       std::size_t nielem,
     124                 :            :       std::vector< std::unordered_set< std::size_t > >& inbox ) const
     125                 :            :     {
     126                 :        294 :       tk::BoxElems< eq >(geoElem, nielem, inbox);
     127                 :        294 :     }
     128                 :            : 
     129                 :            :     //! Initalize the compressible flow equations, prepare for time integration
     130                 :            :     //! \param[in] L Block diagonal mass matrix
     131                 :            :     //! \param[in] inpoel Element-node connectivity
     132                 :            :     //! \param[in] coord Array of nodal coordinates
     133                 :            :     //! \param[in] inbox List of elements at which box user ICs are set for
     134                 :            :     //!   each IC box
     135                 :            :     //! \param[in] elemblkid Element ids associated with mesh block ids where
     136                 :            :     //!   user ICs are set
     137                 :            :     //! \param[in,out] unk Array of unknowns
     138                 :            :     //! \param[in] t Physical time
     139                 :            :     //! \param[in] nielem Number of internal elements
     140                 :        294 :     void initialize( const tk::Fields& L,
     141                 :            :       const std::vector< std::size_t >& inpoel,
     142                 :            :       const tk::UnsMesh::Coords& coord,
     143                 :            :       const std::vector< std::unordered_set< std::size_t > >& inbox,
     144                 :            :       const std::unordered_map< std::size_t, std::set< std::size_t > >&
     145                 :            :         elemblkid,
     146                 :            :       tk::Fields& unk,
     147                 :            :       tk::real t,
     148                 :            :       const std::size_t nielem ) const
     149                 :            :     {
     150 [ +  - ][ +  - ]:        294 :       tk::initialize( m_ncomp, m_mat_blk, L, inpoel, coord,
     151                 :            :                       Problem::initialize, unk, t, nielem );
     152                 :            : 
     153                 :        294 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     154                 :        294 :       const auto& ic = g_inputdeck.get< tag::ic >();
     155                 :        294 :       const auto& icbox = ic.get< tag::box >();
     156                 :        294 :       const auto& icmbk = ic.get< tag::meshblock >();
     157                 :            : 
     158                 :        294 :       const auto& bgpre = ic.get< tag::pressure >();
     159                 :        294 :       const auto& bgtemp = ic.get< tag::temperature >();
     160                 :            : 
     161                 :            :       // Set initial conditions inside user-defined IC boxes and mesh blocks
     162         [ +  - ]:        588 :       std::vector< tk::real > s(m_ncomp, 0.0);
     163         [ +  + ]:      16466 :       for (std::size_t e=0; e<nielem; ++e) {
     164                 :            :         // inside user-defined box
     165         [ -  + ]:      16172 :         if (!icbox.empty()) {
     166                 :          0 :           std::size_t bcnt = 0;
     167         [ -  - ]:          0 :           for (const auto& b : icbox) {   // for all boxes
     168 [ -  - ][ -  - ]:          0 :             if (inbox.size() > bcnt && inbox[bcnt].find(e) != inbox[bcnt].end())
         [ -  - ][ -  - ]
     169                 :            :             {
     170         [ -  - ]:          0 :               std::vector< tk::real > box
     171                 :          0 :                 { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
     172                 :          0 :                   b.template get< tag::ymin >(), b.template get< tag::ymax >(),
     173                 :          0 :                   b.template get< tag::zmin >(), b.template get< tag::zmax >() };
     174                 :          0 :               auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
     175         [ -  - ]:          0 :               for (std::size_t c=0; c<m_ncomp; ++c) {
     176                 :          0 :                 auto mark = c*rdof;
     177         [ -  - ]:          0 :                 s[c] = unk(e,mark);
     178                 :            :                 // set high-order DOFs to zero
     179         [ -  - ]:          0 :                 for (std::size_t i=1; i<rdof; ++i)
     180         [ -  - ]:          0 :                   unk(e,mark+i) = 0.0;
     181                 :            :               }
     182         [ -  - ]:          0 :               initializeBox<ctr::boxList>( m_mat_blk, V_ex, t, b, bgpre,
     183                 :            :                 bgtemp, s );
     184                 :            :               // store box-initialization in solution vector
     185         [ -  - ]:          0 :               for (std::size_t c=0; c<m_ncomp; ++c) {
     186                 :          0 :                 auto mark = c*rdof;
     187         [ -  - ]:          0 :                 unk(e,mark) = s[c];
     188                 :            :               }
     189                 :            :             }
     190                 :          0 :             ++bcnt;
     191                 :            :           }
     192                 :            :         }
     193                 :            : 
     194                 :            :         // inside user-specified mesh blocks
     195         [ -  + ]:      16172 :         for (const auto& b : icmbk) { // for all blocks
     196                 :          0 :           auto blid = b.get< tag::blockid >();
     197                 :          0 :           auto V_ex = b.get< tag::volume >();
     198 [ -  - ][ -  - ]:          0 :           if (elemblkid.find(blid) != elemblkid.end()) {
     199         [ -  - ]:          0 :             const auto& elset = tk::cref_find(elemblkid, blid);
     200 [ -  - ][ -  - ]:          0 :             if (elset.find(e) != elset.end()) {
     201         [ -  - ]:          0 :               initializeBox<ctr::meshblockList>( m_mat_blk, V_ex, t, b,
     202                 :            :                 bgpre, bgtemp, s );
     203                 :            :               // store initialization in solution vector
     204         [ -  - ]:          0 :               for (std::size_t c=0; c<m_ncomp; ++c) {
     205                 :          0 :                 auto mark = c*rdof;
     206         [ -  - ]:          0 :                 unk(e,mark) = s[c];
     207                 :            :               }
     208                 :            :             }
     209                 :            :           }
     210                 :            :         }
     211                 :            :       }
     212                 :        294 :     }
     213                 :            : 
     214                 :            :     //! Compute the left hand side block-diagonal mass matrix
     215                 :            :     //! \param[in] geoElem Element geometry array
     216                 :            :     //! \param[in,out] l Block diagonal mass matrix
     217                 :        294 :     void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
     218                 :        294 :       const auto nelem = geoElem.nunk();
     219         [ +  + ]:      52112 :       for (std::size_t e=0; e<nelem; ++e)
     220         [ +  + ]:     664538 :         for (ncomp_t c=0; c<m_ncomp; ++c)
     221                 :     612720 :           l(e, c) = geoElem(e,0);
     222                 :        294 :     }
     223                 :            : 
     224                 :            :     //! Update the primitives for this PDE system
     225                 :            :     //! \param[in] unk Array of unknowns
     226                 :            :     //! \param[in,out] prim Array of primitives
     227                 :            :     //! \param[in] nielem Number of internal elements
     228                 :            :     //! \details This function computes and stores the dofs for primitive
     229                 :            :     //!   quantities, which are required for obtaining reconstructed states used
     230                 :            :     //!   in the Riemann solver. See /PDE/Riemann/AUSM.hpp, where the
     231                 :            :     //!   normal velocity for advection is calculated from independently
     232                 :            :     //!   reconstructed velocities.
     233                 :       1562 :     void updatePrimitives( const tk::Fields& unk,
     234                 :            :                            tk::Fields& prim,
     235                 :            :                            std::size_t nielem ) const
     236                 :            :     {
     237                 :       1562 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     238                 :       1562 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     239                 :            : 
     240 [ -  + ][ -  - ]:       1562 :       Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     241                 :            :               "vector and primitive vector at recent time step incorrect" );
     242 [ -  + ][ -  - ]:       1562 :       Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
     243                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     244 [ -  + ][ -  - ]:       1562 :       Assert( prim.nprop() == rdof*nprim(), "Number of components in vector of "
         [ -  - ][ -  - ]
                 [ -  - ]
     245                 :            :               "primitive quantities must equal "+ std::to_string(rdof*nprim()) );
     246                 :            : 
     247         [ +  + ]:     221894 :       for (std::size_t e=0; e<nielem; ++e)
     248                 :            :       {
     249                 :            :         // cell-average bulk density
     250                 :     220332 :         tk::real rhob(0.0);
     251         [ +  + ]:     726696 :         for (std::size_t k=0; k<nmat; ++k)
     252                 :            :         {
     253         [ +  - ]:     506364 :           rhob += unk(e, densityDofIdx(nmat, k, rdof, 0));
     254                 :            :         }
     255                 :            : 
     256                 :            :         // cell-average velocity
     257                 :            :         std::array< tk::real, 3 >
     258         [ +  - ]:     220332 :           vel{{ unk(e, momentumDofIdx(nmat, 0, rdof, 0))/rhob,
     259         [ +  - ]:     220332 :                 unk(e, momentumDofIdx(nmat, 1, rdof, 0))/rhob,
     260         [ +  - ]:     220332 :                 unk(e, momentumDofIdx(nmat, 2, rdof, 0))/rhob }};
     261                 :            : 
     262         [ +  + ]:     881328 :         for (std::size_t idir=0; idir<3; ++idir)
     263                 :            :         {
     264         [ +  - ]:     660996 :           prim(e, velocityDofIdx(nmat, idir, rdof, 0)) = vel[idir];
     265         [ +  + ]:    2643984 :           for (std::size_t idof=1; idof<rdof; ++idof)
     266         [ +  - ]:    1982988 :             prim(e, velocityDofIdx(nmat, idir, rdof, idof)) = 0.0;
     267                 :            :         }
     268                 :            : 
     269                 :            :         // cell-average material pressure
     270         [ +  + ]:     726696 :         for (std::size_t k=0; k<nmat; ++k)
     271                 :            :         {
     272         [ +  - ]:     506364 :           tk::real arhomat = unk(e, densityDofIdx(nmat, k, rdof, 0));
     273         [ +  - ]:     506364 :           tk::real arhoemat = unk(e, energyDofIdx(nmat, k, rdof, 0));
     274         [ +  - ]:     506364 :           tk::real alphamat = unk(e, volfracDofIdx(nmat, k, rdof, 0));
     275 [ +  - ][ +  - ]:     506364 :           auto gmat = getDeformGrad(nmat, k, unk.extract(e));
     276         [ +  - ]:     506364 :           prim(e, pressureDofIdx(nmat, k, rdof, 0)) =
     277         [ +  - ]:     506364 :             m_mat_blk[k].compute< EOS::pressure >( arhomat, vel[0], vel[1],
     278                 :            :             vel[2], arhoemat, alphamat, k, gmat );
     279         [ +  - ]:     506364 :           prim(e, pressureDofIdx(nmat, k, rdof, 0)) =
     280                 :     506364 :             constrain_pressure( m_mat_blk,
     281 [ +  - ][ +  - ]:     506364 :             prim(e, pressureDofIdx(nmat, k, rdof, 0)), arhomat, alphamat, k);
     282         [ +  + ]:    2025456 :           for (std::size_t idof=1; idof<rdof; ++idof)
     283         [ +  - ]:    1519092 :             prim(e, pressureDofIdx(nmat, k, rdof, idof)) = 0.0;
     284                 :            :         }
     285                 :            :       }
     286                 :       1562 :     }
     287                 :            : 
     288                 :            :     //! Clean up the state of trace materials for this PDE system
     289                 :            :     //! \param[in] t Physical time
     290                 :            :     //! \param[in] geoElem Element geometry array
     291                 :            :     //! \param[in,out] unk Array of unknowns
     292                 :            :     //! \param[in,out] prim Array of primitives
     293                 :            :     //! \param[in] nielem Number of internal elements
     294                 :            :     //! \details This function cleans up the state of materials present in trace
     295                 :            :     //!   quantities in each cell. Specifically, the state of materials with
     296                 :            :     //!   very low volume-fractions in a cell is replaced by the state of the
     297                 :            :     //!   material which is present in the largest quantity in that cell. This
     298                 :            :     //!   becomes necessary when shocks pass through cells which contain a very
     299                 :            :     //!   small amount of material. The state of that tiny material might
     300                 :            :     //!   become unphysical and cause solution to diverge; thus requiring such
     301                 :            :     //!   a "reset".
     302                 :       1268 :     void cleanTraceMaterial( tk::real t,
     303                 :            :                              const tk::Fields& geoElem,
     304                 :            :                              tk::Fields& unk,
     305                 :            :                              tk::Fields& prim,
     306                 :            :                              std::size_t nielem ) const
     307                 :            :     {
     308                 :       1268 :       [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
     309                 :       1268 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     310                 :            : 
     311 [ -  + ][ -  - ]:       1268 :       Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     312                 :            :               "vector and primitive vector at recent time step incorrect" );
     313 [ -  + ][ -  - ]:       1268 :       Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
     314                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     315 [ -  + ][ -  - ]:       1268 :       Assert( prim.nprop() == rdof*nprim(), "Number of components in vector of "
         [ -  - ][ -  - ]
                 [ -  - ]
     316                 :            :               "primitive quantities must equal "+ std::to_string(rdof*nprim()) );
     317 [ -  + ][ -  - ]:       1268 :       Assert( (g_inputdeck.get< tag::ndof >()) <= 4, "High-order "
         [ -  - ][ -  - ]
     318                 :            :               "discretizations not set up for multimat cleanTraceMaterial()" );
     319                 :            : 
     320                 :       1268 :       auto neg_density = cleanTraceMultiMat(t, nielem, m_mat_blk, geoElem, nmat,
     321                 :            :         unk, prim);
     322                 :            : 
     323 [ -  + ][ -  - ]:       1268 :       if (neg_density) Throw("Negative partial density.");
         [ -  - ][ -  - ]
     324                 :       1268 :     }
     325                 :            : 
     326                 :            :     //! Reconstruct second-order solution from first-order
     327                 :            :     //! \param[in] geoElem Element geometry array
     328                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     329                 :            :     //! \param[in] esup Elements-surrounding-nodes connectivity
     330                 :            :     //! \param[in] inpoel Element-node connectivity
     331                 :            :     //! \param[in] coord Array of nodal coordinates
     332                 :            :     //! \param[in,out] U Solution vector at recent time step
     333                 :            :     //! \param[in,out] P Vector of primitives at recent time step
     334                 :       1268 :     void reconstruct( const tk::Fields& geoElem,
     335                 :            :                       const inciter::FaceData& fd,
     336                 :            :                       const std::map< std::size_t, std::vector< std::size_t > >&
     337                 :            :                         esup,
     338                 :            :                       const std::vector< std::size_t >& inpoel,
     339                 :            :                       const tk::UnsMesh::Coords& coord,
     340                 :            :                       tk::Fields& U,
     341                 :            :                       tk::Fields& P ) const
     342                 :            :     {
     343                 :       1268 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     344                 :       1268 :       const auto nelem = fd.Esuel().size()/4;
     345                 :       1268 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     346                 :            : 
     347 [ -  + ][ -  - ]:       1268 :       Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
     348                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     349                 :            : 
     350                 :            :       //----- reconstruction of conserved quantities -----
     351                 :            :       //--------------------------------------------------
     352                 :            :       // specify how many variables need to be reconstructed
     353                 :       2536 :       std::vector< std::size_t > vars;
     354         [ +  + ]:       4972 :       for (std::size_t k=0; k<nmat; ++k) {
     355         [ +  - ]:       3704 :         vars.push_back(volfracIdx(nmat,k));
     356         [ +  - ]:       3704 :         vars.push_back(densityIdx(nmat,k));
     357                 :            :       }
     358                 :            : 
     359         [ +  + ]:     205428 :       for (std::size_t e=0; e<nelem; ++e)
     360                 :            :       {
     361                 :            :         // 1. solve 3x3 least-squares system
     362                 :            :         // Reconstruct second-order dofs of volume-fractions in Taylor space
     363                 :            :         // using nodal-stencils, for a good interface-normal estimate
     364         [ +  - ]:     204160 :         tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, U, vars );
     365                 :            : 
     366                 :            :         // 2. transform reconstructed derivatives to Dubiner dofs
     367         [ +  - ]:     204160 :         tk::transform_P0P1(rdof, e, inpoel, coord, U, vars);
     368                 :            :       }
     369                 :            : 
     370                 :            :       //----- reconstruction of primitive quantities -----
     371                 :            :       //--------------------------------------------------
     372                 :            :       // For multimat, conserved and primitive quantities are reconstructed
     373                 :            :       // separately.
     374                 :       1268 :       vars.clear();
     375 [ +  + ][ +  - ]:       8776 :       for (std::size_t c=0; c<nprim(); ++c) vars.push_back(c);
     376         [ +  + ]:     205428 :       for (std::size_t e=0; e<nelem; ++e)
     377                 :            :       {
     378                 :            :         // 1.
     379                 :            :         // Reconstruct second-order dofs of volume-fractions in Taylor space
     380                 :            :         // using nodal-stencils, for a good interface-normal estimate
     381         [ +  - ]:     204160 :         tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, P, vars );
     382                 :            : 
     383                 :            :         // 2.
     384         [ +  - ]:     204160 :         tk::transform_P0P1(rdof, e, inpoel, coord, P, vars );
     385                 :            :       }
     386                 :       1268 :     }
     387                 :            : 
     388                 :            :     //! Limit second-order solution, and primitive quantities separately
     389                 :            :     //! \param[in] geoFace Face geometry array
     390                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     391                 :            :     //! \param[in] esup Elements-surrounding-nodes connectivity
     392                 :            :     //! \param[in] inpoel Element-node connectivity
     393                 :            :     //! \param[in] coord Array of nodal coordinates
     394                 :            :     //! \param[in] srcFlag Whether the energy source was added
     395                 :            :     //! \param[in,out] U Solution vector at recent time step
     396                 :            :     //! \param[in,out] P Vector of primitives at recent time step
     397                 :       1268 :     void limit( const tk::Fields& geoFace,
     398                 :            :                 const inciter::FaceData& fd,
     399                 :            :                 const std::map< std::size_t, std::vector< std::size_t > >& esup,
     400                 :            :                 const std::vector< std::size_t >& inpoel,
     401                 :            :                 const tk::UnsMesh::Coords& coord,
     402                 :            :                 const std::vector< int >& srcFlag,
     403                 :            :                 tk::Fields& U,
     404                 :            :                 tk::Fields& P ) const
     405                 :            :     {
     406 [ -  + ][ -  - ]:       1268 :       Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     407                 :            :               "vector and primitive vector at recent time step incorrect" );
     408                 :            : 
     409                 :       1268 :       const auto limiter = g_inputdeck.get< tag::limiter >();
     410                 :       1268 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     411                 :            :       const auto& solidx = g_inputdeck.get<
     412                 :       1268 :         tag::matidxmap, tag::solidx >();
     413                 :            : 
     414                 :            :       // limit vectors of conserved and primitive quantities
     415         [ +  - ]:       1268 :       if (limiter == ctr::LimiterType::VERTEXBASEDP1)
     416                 :            :       {
     417                 :       1268 :         VertexBasedMultiMat_FV( esup, inpoel, fd.Esuel().size()/4,
     418                 :            :           coord, srcFlag, solidx, U, P, nmat );
     419                 :       1268 :         PositivityPreservingMultiMat_FV( inpoel, fd.Esuel().size()/4, nmat,
     420                 :       1268 :           m_mat_blk, coord, geoFace, U, P );
     421                 :            :       }
     422         [ -  - ]:          0 :       else if (limiter != ctr::LimiterType::NOLIMITER)
     423                 :            :       {
     424 [ -  - ][ -  - ]:          0 :         Throw("Limiter type not configured for multimat.");
                 [ -  - ]
     425                 :            :       }
     426                 :       1268 :     }
     427                 :            : 
     428                 :            :     //! Apply CPL to the conservative variable solution for this PDE system
     429                 :            :     //! \param[in] prim Array of primitive variables
     430                 :            :     //! \param[in] geoElem Element geometry array
     431                 :            :     //! \param[in] inpoel Element-node connectivity
     432                 :            :     //! \param[in] coord Array of nodal coordinates
     433                 :            :     //! \param[in,out] unk Array of conservative variables
     434                 :            :     //! \param[in] nielem Number of internal elements
     435                 :            :     //! \details This function applies CPL to obtain consistent dofs for
     436                 :            :     //!   conservative quantities based on the limited primitive quantities.
     437                 :            :     //!   See Pandare et al. (2023). On the Design of Stable,
     438                 :            :     //!   Consistent, and Conservative High-Order Methods for Multi-Material
     439                 :            :     //!   Hydrodynamics. J Comp Phys, 112313.
     440                 :            :     void CPL( const tk::Fields& prim,
     441                 :            :       const tk::Fields& geoElem,
     442                 :            :       const std::vector< std::size_t >& inpoel,
     443                 :            :       const tk::UnsMesh::Coords& coord,
     444                 :            :       tk::Fields& unk,
     445                 :            :       std::size_t nielem ) const
     446                 :            :     {
     447                 :            :       [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
     448                 :            :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     449                 :            : 
     450                 :            :       Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
     451                 :            :               "vector and primitive vector at recent time step incorrect" );
     452                 :            :       Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
     453                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     454                 :            :       Assert( prim.nprop() == rdof*nprim(), "Number of components in vector of "
     455                 :            :               "primitive quantities must equal "+ std::to_string(rdof*nprim()) );
     456                 :            : 
     457                 :            :       correctLimConservMultiMat(nielem, m_mat_blk, nmat, inpoel,
     458                 :            :         coord, geoElem, prim, unk);
     459                 :            :     }
     460                 :            : 
     461                 :            :     //! Compute right hand side
     462                 :            :     //! \param[in] t Physical time
     463                 :            :     //! \param[in] geoFace Face geometry array
     464                 :            :     //! \param[in] geoElem Element geometry array
     465                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     466                 :            :     //! \param[in] inpoel Element-node connectivity
     467                 :            :     //! \param[in] coord Array of nodal coordinates
     468                 :            :     //! \param[in] elemblkid Element ids associated with mesh block ids where
     469                 :            :     //!   user ICs are set
     470                 :            :     //! \param[in] U Solution vector at recent time step
     471                 :            :     //! \param[in] P Primitive vector at recent time step
     472                 :            :     //! \param[in,out] R Right-hand side vector computed
     473                 :            :     //! \param[in,out] srcFlag Whether the energy source was added
     474                 :       1268 :     void rhs( tk::real t,
     475                 :            :               const tk::Fields& geoFace,
     476                 :            :               const tk::Fields& geoElem,
     477                 :            :               const inciter::FaceData& fd,
     478                 :            :               const std::vector< std::size_t >& inpoel,
     479                 :            :               const tk::UnsMesh::Coords& coord,
     480                 :            :               const std::unordered_map< std::size_t, std::set< std::size_t > >&
     481                 :            :                 elemblkid,
     482                 :            :               const tk::Fields& U,
     483                 :            :               const tk::Fields& P,
     484                 :            :               tk::Fields& R,
     485                 :            :               std::vector< int >& srcFlag ) const
     486                 :            :     {
     487                 :       1268 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     488                 :       1268 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     489                 :       1268 :       const auto intsharp =
     490                 :       1268 :         g_inputdeck.get< tag::multimat, tag::intsharp >();
     491                 :       1268 :       auto viscous = g_inputdeck.get< tag::multimat, tag::viscous >();
     492                 :            : 
     493                 :       1268 :       const auto nelem = fd.Esuel().size()/4;
     494                 :            : 
     495 [ -  + ][ -  - ]:       1268 :       Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     496                 :            :               "vector and primitive vector at recent time step incorrect" );
     497 [ -  + ][ -  - ]:       1268 :       Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     498                 :            :               "vector and right-hand side at recent time step incorrect" );
     499 [ -  + ][ -  - ]:       1268 :       Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
     500                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     501 [ -  + ][ -  - ]:       1268 :       Assert( P.nprop() == rdof*nprim(), "Number of components in primitive "
         [ -  - ][ -  - ]
                 [ -  - ]
     502                 :            :               "vector must equal "+ std::to_string(rdof*nprim()) );
     503 [ -  + ][ -  - ]:       1268 :       Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
         [ -  - ][ -  - ]
     504                 :            :               "Mismatch in inpofa size" );
     505                 :            : 
     506                 :            :       // set rhs to zero
     507         [ +  - ]:       1268 :       R.fill(0.0);
     508                 :            : 
     509                 :            :       // configure a no-op lambda for prescribed velocity
     510                 :     476216 :       auto velfn = []( ncomp_t, tk::real, tk::real, tk::real, tk::real ){
     511                 :     476216 :         return tk::VelFn::result_type(); };
     512                 :            : 
     513                 :            :       // compute internal surface flux (including non-conservative) integrals
     514         [ +  - ]:       2536 :       tk::surfIntFV( nmat, m_mat_blk, t, rdof, inpoel,
     515         [ +  - ]:       1268 :                      coord, fd, geoFace, geoElem, m_riemann, velfn, U, P,
     516                 :            :                      srcFlag, R, intsharp );
     517                 :            :       // compute internal surface viscous flux integrals
     518         [ +  - ]:       2536 :       tk::Fields T( U.nunk(), rdof*nmat );
     519         [ -  + ]:       1268 :       if (viscous) {
     520         [ -  - ]:          0 :         computeTemperaturesFV( m_mat_blk, nmat, inpoel, coord, geoElem,
     521                 :            :                                U, P, srcFlag, T );
     522                 :            : 
     523         [ -  - ]:          0 :         tk::surfIntViscousFV( nmat, m_mat_blk, rdof, inpoel,
     524                 :            :                               coord, fd, geoFace, geoElem, U, P, T,
     525                 :            :                               srcFlag, R, intsharp );
     526                 :            :       }
     527                 :            : 
     528                 :            :       // compute boundary surface flux (including non-conservative) integrals
     529         [ +  + ]:       8876 :       for (const auto& b : m_bc) {
     530 [ +  - ][ +  - ]:      15216 :         tk::bndSurfIntFV( nmat, m_mat_blk, rdof, std::get<0>(b),
     531                 :       7608 :                           fd, geoFace, geoElem, inpoel, coord, t, m_riemann,
     532                 :       7608 :                           velfn, std::get<1>(b), U, P, srcFlag, R, intsharp );
     533         [ -  + ]:       7608 :         if (viscous)
     534         [ -  - ]:          0 :           tk::bndSurfIntViscousFV( nmat, m_mat_blk, rdof, std::get<0>(b),
     535                 :            :                                    fd, geoFace, geoElem, inpoel, coord, t,
     536                 :          0 :                                    std::get<1>(b), std::get<2>(b), U, P, T,
     537                 :            :                                    srcFlag, R, intsharp );
     538                 :            :       }
     539                 :            : 
     540                 :            :       // compute optional source term
     541 [ +  - ][ +  - ]:       1268 :       tk::srcIntFV( m_mat_blk, t, fd.Esuel().size()/4,
     542                 :            :                     geoElem, Problem::src, R, nmat );
     543                 :            : 
     544                 :            :       // compute finite pressure relaxation terms
     545         [ +  + ]:       1268 :       if (g_inputdeck.get< tag::multimat, tag::prelax >())
     546                 :            :       {
     547                 :       1218 :         const auto ct = g_inputdeck.get< tag::multimat,
     548                 :       1218 :                                          tag::prelax_timescale >();
     549         [ +  - ]:       1218 :         tk::pressureRelaxationIntFV( nmat, m_mat_blk, rdof,
     550                 :            :                                      nelem, inpoel, coord, geoElem, U, P, ct,
     551                 :            :                                      R );
     552                 :            :       }
     553                 :            : 
     554                 :            :       // compute external (energy) sources
     555         [ -  - ]:       1268 :       m_physics.physSrc(nmat, t, geoElem, elemblkid, R, srcFlag);
     556                 :       1268 :     }
     557                 :            : 
     558                 :            :     //! Compute the minimum time step size
     559                 :            : //    //! \param[in] fd Face connectivity and boundary conditions object
     560                 :            : //    //! \param[in] geoFace Face geometry array
     561                 :            :     //! \param[in] geoElem Element geometry array
     562                 :            :     //! \param[in] U Solution vector at recent time step
     563                 :            :     //! \param[in] P Vector of primitive quantities at recent time step
     564                 :            :     //! \param[in] nielem Number of internal elements
     565                 :            :     //! \param[in] srcFlag Whether the energy source was added
     566                 :            :     //! \param[in,out] local_dte Time step size for each element (for local
     567                 :            :     //!   time stepping)
     568                 :            :     //! \return Minimum time step size
     569                 :            :     //! \details The allowable dt is calculated by looking at the maximum
     570                 :            :     //!   wave-speed in elements surrounding each face, times the area of that
     571                 :            :     //!   face. Once the maximum of this quantity over the mesh is determined,
     572                 :            :     //!   the volume of each cell is divided by this quantity. A minimum of this
     573                 :            :     //!   ratio is found over the entire mesh, which gives the allowable dt.
     574                 :         50 :     tk::real dt( const inciter::FaceData& /*fd*/,
     575                 :            :                  const tk::Fields& /*geoFace*/,
     576                 :            :                  const tk::Fields& geoElem,
     577                 :            :                  const tk::Fields& U,
     578                 :            :                  const tk::Fields& P,
     579                 :            :                  const std::size_t nielem,
     580                 :            :                  const std::vector< int >& srcFlag,
     581                 :            :                  std::vector< tk::real >& local_dte ) const
     582                 :            :     {
     583                 :         50 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     584                 :         50 :       auto viscous = g_inputdeck.get< tag::multimat, tag::viscous >();
     585                 :            : 
     586                 :            :       // obtain dt restrictions from all physics
     587         [ +  - ]:         50 :       auto dt_e = timeStepSizeMultiMatFV(m_mat_blk, geoElem, nielem, nmat, U,
     588                 :            :         P, local_dte);
     589         [ -  + ]:         50 :       if (viscous)
     590         [ -  - ]:          0 :         dt_e = std::min(dt_e, timeStepSizeViscousFV(geoElem, nielem, nmat, U));
     591         [ -  - ]:         50 :       auto dt_p = m_physics.dtRestriction(geoElem, nielem, srcFlag);
     592                 :            : 
     593                 :         50 :       return std::min(dt_e, dt_p);
     594                 :            :     }
     595                 :            : 
     596                 :            :     //! Extract the velocity field at cell nodes. Currently unused.
     597                 :            :     //! \param[in] U Solution vector at recent time step
     598                 :            :     //! \param[in] N Element node indices
     599                 :            :     //! \return Array of the four values of the velocity field
     600                 :            :     std::array< std::array< tk::real, 4 >, 3 >
     601                 :            :     velocity( const tk::Fields& U,
     602                 :            :               const std::array< std::vector< tk::real >, 3 >&,
     603                 :            :               const std::array< std::size_t, 4 >& N ) const
     604                 :            :     {
     605                 :            :       const auto rdof = g_inputdeck.get< tag::rdof >();
     606                 :            :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     607                 :            : 
     608                 :            :       std::array< std::array< tk::real, 4 >, 3 > v;
     609                 :            :       v[0] = U.extract( momentumDofIdx(nmat, 0, rdof, 0), N );
     610                 :            :       v[1] = U.extract( momentumDofIdx(nmat, 1, rdof, 0), N );
     611                 :            :       v[2] = U.extract( momentumDofIdx(nmat, 2, rdof, 0), N );
     612                 :            : 
     613                 :            :       std::vector< std::array< tk::real, 4 > > ar;
     614                 :            :       ar.resize(nmat);
     615                 :            :       for (std::size_t k=0; k<nmat; ++k)
     616                 :            :         ar[k] = U.extract( densityDofIdx(nmat, k, rdof, 0), N );
     617                 :            : 
     618                 :            :       std::array< tk::real, 4 > r{{ 0.0, 0.0, 0.0, 0.0 }};
     619                 :            :       for (std::size_t i=0; i<r.size(); ++i) {
     620                 :            :         for (std::size_t k=0; k<nmat; ++k)
     621                 :            :           r[i] += ar[k][i];
     622                 :            :       }
     623                 :            : 
     624                 :            :       std::transform( r.begin(), r.end(), v[0].begin(), v[0].begin(),
     625                 :            :                       []( tk::real s, tk::real& d ){ return d /= s; } );
     626                 :            :       std::transform( r.begin(), r.end(), v[1].begin(), v[1].begin(),
     627                 :            :                       []( tk::real s, tk::real& d ){ return d /= s; } );
     628                 :            :       std::transform( r.begin(), r.end(), v[2].begin(), v[2].begin(),
     629                 :            :                       []( tk::real s, tk::real& d ){ return d /= s; } );
     630                 :            :       return v;
     631                 :            :     }
     632                 :            : 
     633                 :            :     //! Return a map that associates user-specified strings to functions
     634                 :            :     //! \return Map that associates user-specified strings to functions that
     635                 :            :     //!   compute relevant quantities to be output to file
     636                 :          8 :     std::map< std::string, tk::GetVarFn > OutVarFn() const
     637                 :          8 :     { return MultiMatOutVarFn(); }
     638                 :            : 
     639                 :            :     //! Return analytic field names to be output to file
     640                 :            :     //! \return Vector of strings labelling analytic fields output in file
     641                 :          0 :     std::vector< std::string > analyticFieldNames() const {
     642                 :          0 :       auto nmat = g_inputdeck.get< eq, tag::nmat >();
     643                 :            : 
     644                 :          0 :       return MultiMatFieldNames(nmat);
     645                 :            :     }
     646                 :            : 
     647                 :            :     //! Return surface field names to be output to file
     648                 :            :     //! \return Vector of strings labelling surface fields output in file
     649                 :          4 :     std::vector< std::string > surfNames() const
     650                 :          4 :     { return MultiMatSurfNames(); }
     651                 :            : 
     652                 :            :     //! Return time history field names to be output to file
     653                 :            :     //! \return Vector of strings labelling time history fields output in file
     654                 :          0 :     std::vector< std::string > histNames() const {
     655                 :          0 :       return MultiMatHistNames();
     656                 :            :     }
     657                 :            : 
     658                 :            :     //! Return surface field output going to file
     659                 :            :     std::vector< std::vector< tk::real > >
     660                 :          4 :     surfOutput( const inciter::FaceData& fd,
     661                 :            :       const tk::Fields& U,
     662                 :            :       const tk::Fields& P ) const
     663                 :            :     {
     664                 :          4 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     665                 :          4 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     666                 :            : 
     667                 :          4 :       return MultiMatSurfOutput( nmat, rdof, fd, U, P );
     668                 :            :     }
     669                 :            : 
     670                 :            :     //! Return time history field output evaluated at time history points
     671                 :            :     //! \param[in] h History point data
     672                 :            :     //! \param[in] inpoel Element-node connectivity
     673                 :            :     //! \param[in] coord Array of nodal coordinates
     674                 :            :     //! \param[in] U Array of unknowns
     675                 :            :     //! \param[in] P Array of primitive quantities
     676                 :            :     //! \return Vector of time history output of bulk flow quantities (density,
     677                 :            :     //!   velocity, total energy, pressure, and volume fraction) evaluated at 
     678                 :            :     //!   time history points
     679                 :            :     std::vector< std::vector< tk::real > >
     680                 :          0 :     histOutput( const std::vector< HistData >& h,
     681                 :            :                 const std::vector< std::size_t >& inpoel,
     682                 :            :                 const tk::UnsMesh::Coords& coord,
     683                 :            :                 const tk::Fields& U,
     684                 :            :                 const tk::Fields& P ) const
     685                 :            :     {
     686                 :          0 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     687                 :          0 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     688                 :            : 
     689                 :          0 :       const auto& x = coord[0];
     690                 :          0 :       const auto& y = coord[1];
     691                 :          0 :       const auto& z = coord[2];
     692                 :            : 
     693         [ -  - ]:          0 :       std::vector< std::vector< tk::real > > Up(h.size());
     694                 :            : 
     695                 :          0 :       std::size_t j = 0;
     696         [ -  - ]:          0 :       for (const auto& p : h) {
     697                 :          0 :         auto e = p.get< tag::elem >();
     698                 :          0 :         auto chp = p.get< tag::coord >();
     699                 :            : 
     700                 :            :         // Evaluate inverse Jacobian
     701                 :          0 :         std::array< std::array< tk::real, 3>, 4 > cp{{
     702                 :            :           {{ x[inpoel[4*e  ]], y[inpoel[4*e  ]], z[inpoel[4*e  ]] }},
     703                 :          0 :           {{ x[inpoel[4*e+1]], y[inpoel[4*e+1]], z[inpoel[4*e+1]] }},
     704                 :          0 :           {{ x[inpoel[4*e+2]], y[inpoel[4*e+2]], z[inpoel[4*e+2]] }},
     705                 :          0 :           {{ x[inpoel[4*e+3]], y[inpoel[4*e+3]], z[inpoel[4*e+3]] }} }};
     706                 :          0 :         auto J = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
     707                 :            : 
     708                 :            :         // evaluate solution at history-point
     709                 :          0 :         std::array< tk::real, 3 > dc{{chp[0]-cp[0][0], chp[1]-cp[0][1],
     710                 :          0 :           chp[2]-cp[0][2]}};
     711         [ -  - ]:          0 :         auto B = tk::eval_basis(rdof, tk::dot(J[0],dc), tk::dot(J[1],dc),
     712                 :          0 :           tk::dot(J[2],dc));
     713         [ -  - ]:          0 :         auto uhp = eval_state(m_ncomp, rdof, rdof, e, U, B);
     714         [ -  - ]:          0 :         auto php = eval_state(nprim(), rdof, rdof, e, P, B);
     715                 :            : 
     716                 :            :         // store solution in history output vector
     717         [ -  - ]:          0 :         Up[j].resize(6+nmat, 0.0);
     718         [ -  - ]:          0 :         for (std::size_t k=0; k<nmat; ++k) {
     719                 :          0 :           Up[j][0] += uhp[densityIdx(nmat,k)];
     720                 :          0 :           Up[j][4] += uhp[energyIdx(nmat,k)];
     721                 :          0 :           Up[j][5] += php[pressureIdx(nmat,k)];
     722                 :          0 :           Up[j][6+k] = uhp[volfracIdx(nmat,k)];
     723                 :            :         }
     724                 :          0 :         Up[j][1] = php[velocityIdx(nmat,0)];
     725                 :          0 :         Up[j][2] = php[velocityIdx(nmat,1)];
     726                 :          0 :         Up[j][3] = php[velocityIdx(nmat,2)];
     727                 :          0 :         ++j;
     728                 :            :       }
     729                 :            : 
     730                 :          0 :       return Up;
     731                 :            :     }
     732                 :            : 
     733                 :            :     //! Return names of integral variables to be output to diagnostics file
     734                 :            :     //! \return Vector of strings labelling integral variables output
     735                 :         20 :     std::vector< std::string > names() const
     736                 :            :     {
     737                 :         20 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     738                 :         20 :       return MultiMatDiagNames(nmat);
     739                 :            :     }
     740                 :            : 
     741                 :            :     //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
     742                 :            :     //! \param[in] xi X-coordinate at which to evaluate the analytic solution
     743                 :            :     //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
     744                 :            :     //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
     745                 :            :     //! \param[in] t Physical time at which to evaluate the analytic solution
     746                 :            :     //! \return Vector of analytic solution at given location and time
     747                 :            :     std::vector< tk::real >
     748                 :          0 :     analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
     749                 :          0 :     { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
     750                 :            : 
     751                 :            :     //! Return analytic solution for conserved variables
     752                 :            :     //! \param[in] xi X-coordinate at which to evaluate the analytic solution
     753                 :            :     //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
     754                 :            :     //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
     755                 :            :     //! \param[in] t Physical time at which to evaluate the analytic solution
     756                 :            :     //! \return Vector of analytic solution at given location and time
     757                 :            :     std::vector< tk::real >
     758                 :          0 :     solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
     759                 :          0 :     { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
     760                 :            : 
     761                 :            :     //! Return cell-averaged specific total energy for an element
     762                 :            :     //! \param[in] e Element id for which total energy is required
     763                 :            :     //! \param[in] unk Vector of conserved quantities
     764                 :            :     //! \return Cell-averaged specific total energy for given element
     765                 :          0 :     tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
     766                 :            :     {
     767                 :          0 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     768                 :          0 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     769                 :            : 
     770                 :          0 :       tk::real sp_te(0.0);
     771                 :            :       // sum each material total energy
     772         [ -  - ]:          0 :       for (std::size_t k=0; k<nmat; ++k) {
     773                 :          0 :         sp_te += unk(e, energyDofIdx(nmat,k,rdof,0));
     774                 :            :       }
     775                 :          0 :       return sp_te;
     776                 :            :     }
     777                 :            : 
     778                 :            :     //! Compute relevant sound speed for output
     779                 :            :     //! \param[in] nielem Number of internal elements
     780                 :            :     //! \param[in] U Solution vector at recent time step
     781                 :            :     //! \param[in] P Primitive vector at recent time step
     782                 :            :     //! \param[in,out] ss Sound speed vector
     783                 :          4 :     void soundspeed(
     784                 :            :       std::size_t nielem,
     785                 :            :       const tk::Fields& U,
     786                 :            :       const tk::Fields& P,
     787                 :            :       std::vector< tk::real >& ss) const
     788                 :            :     {
     789 [ -  + ][ -  - ]:          4 :       Assert( ss.size() == nielem, "Size of sound speed vector incorrect " );
         [ -  - ][ -  - ]
     790                 :            : 
     791                 :          4 :       const auto ndof = g_inputdeck.get< tag::ndof >();
     792                 :          4 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     793                 :          4 :       const auto use_mass_avg =
     794                 :          4 :         g_inputdeck.get< tag::multimat, tag::dt_sos_massavg >();
     795                 :          4 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     796                 :          4 :       std::size_t ncomp = U.nprop()/rdof;
     797                 :          4 :       std::size_t nprim = P.nprop()/rdof;
     798                 :            : 
     799 [ +  - ][ +  - ]:          8 :       std::vector< tk::real > ugp(ncomp, 0.0), pgp(nprim, 0.0);
     800                 :            : 
     801         [ +  + ]:       6068 :       for (std::size_t e=0; e<nielem; ++e) {
     802                 :            :         // basis function at centroid
     803         [ +  - ]:      12128 :         std::vector< tk::real > B(rdof, 0.0);
     804                 :       6064 :         B[0] = 1.0;
     805                 :            : 
     806                 :            :         // get conserved quantities
     807         [ +  - ]:       6064 :         ugp = eval_state(ncomp, rdof, ndof, e, U, B);
     808                 :            :         // get primitive quantities
     809         [ +  - ]:       6064 :         pgp = eval_state(nprim, rdof, ndof, e, P, B);
     810                 :            : 
     811                 :            :         // acoustic speed (this should be consistent with time-step calculation)
     812                 :       6064 :         ss[e] = 0.0;
     813                 :       6064 :         tk::real mixtureDensity = 0.0;
     814         [ +  + ]:      18192 :         for (std::size_t k=0; k<nmat; ++k)
     815                 :            :         {
     816         [ -  + ]:      12128 :           if (use_mass_avg > 0)
     817                 :            :           {
     818                 :            :             // mass averaging SoS
     819                 :          0 :             ss[e] += ugp[densityIdx(nmat,k)]*
     820         [ -  - ]:          0 :               m_mat_blk[k].compute< EOS::soundspeed >(
     821                 :            :               ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
     822                 :            :               ugp[volfracIdx(nmat, k)], k );
     823                 :            : 
     824                 :          0 :             mixtureDensity += ugp[densityIdx(nmat,k)];
     825                 :            :           }
     826                 :            :           else
     827                 :            :           {
     828         [ +  + ]:      12128 :             if (ugp[volfracIdx(nmat, k)] > 1.0e-04)
     829                 :            :             {
     830         [ +  - ]:      12356 :               ss[e] = std::max( ss[e], m_mat_blk[k].compute< EOS::soundspeed >(
     831                 :            :                 ugp[densityIdx(nmat, k)], pgp[pressureIdx(nmat, k)],
     832                 :       6178 :                 ugp[volfracIdx(nmat, k)], k ) );
     833                 :            :             }
     834                 :            :           }
     835                 :            :         }
     836         [ -  + ]:       6064 :         if (use_mass_avg > 0) ss[e] /= mixtureDensity;
     837                 :            :       }
     838                 :          4 :     }
     839                 :            : 
     840                 :            :   private:
     841                 :            :     //! Physics policy
     842                 :            :     const Physics m_physics;
     843                 :            :     //! Number of components in this PDE system
     844                 :            :     const ncomp_t m_ncomp;
     845                 :            :     //! Riemann solver
     846                 :            :     tk::RiemannFluxFn m_riemann;
     847                 :            :     //! BC configuration
     848                 :            :     BCStateFn m_bc;
     849                 :            :     //! EOS material block
     850                 :            :     std::vector< EOS > m_mat_blk;
     851                 :            : 
     852                 :            :     //! Evaluate conservative part of physical flux function for this PDE system
     853                 :            :     //! \param[in] ncomp Number of scalar components in this PDE system
     854                 :            :     //! \param[in] ugp Numerical solution at the Gauss point at which to
     855                 :            :     //!   evaluate the flux
     856                 :            :     //! \return Flux vectors for all components in this PDE system
     857                 :            :     //! \note The function signature must follow tk::FluxFn
     858                 :            :     static tk::FluxFn::result_type
     859                 :            :     flux( ncomp_t ncomp,
     860                 :            :           const std::vector< EOS >& mat_blk,
     861                 :            :           const std::vector< tk::real >& ugp,
     862                 :            :           const std::vector< std::array< tk::real, 3 > >& )
     863                 :            :     {
     864                 :            :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     865                 :            : 
     866                 :            :       return tk::fluxTerms(ncomp, nmat, mat_blk, ugp);
     867                 :            :     }
     868                 :            : 
     869                 :            :     //! \brief Boundary state function providing the left and right state of a
     870                 :            :     //!   face at Dirichlet boundaries
     871                 :            :     //! \param[in] ncomp Number of scalar components in this PDE system
     872                 :            :     //! \param[in] ul Left (domain-internal) state
     873                 :            :     //! \param[in] x X-coordinate at which to compute the states
     874                 :            :     //! \param[in] y Y-coordinate at which to compute the states
     875                 :            :     //! \param[in] z Z-coordinate at which to compute the states
     876                 :            :     //! \param[in] t Physical time
     877                 :            :     //! \return Left and right states for all scalar components in this PDE
     878                 :            :     //!   system
     879                 :            :     //! \note The function signature must follow tk::StateFn. For multimat, the
     880                 :            :     //!   left or right state is the vector of conserved quantities, followed by
     881                 :            :     //!   the vector of primitive quantities appended to it.
     882                 :            :     static tk::StateFn::result_type
     883                 :          0 :     dirichlet( ncomp_t ncomp,
     884                 :            :                const std::vector< EOS >& mat_blk,
     885                 :            :                const std::vector< tk::real >& ul, tk::real x, tk::real y,
     886                 :            :                tk::real z, tk::real t, const std::array< tk::real, 3 >& )
     887                 :            :     {
     888                 :          0 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     889                 :            : 
     890         [ -  - ]:          0 :       auto ur = Problem::initialize( ncomp, mat_blk, x, y, z, t );
     891 [ -  - ][ -  - ]:          0 :       Assert( ur.size() == ncomp, "Incorrect size for boundary state vector" );
         [ -  - ][ -  - ]
     892                 :            : 
     893         [ -  - ]:          0 :       ur.resize(ul.size());
     894                 :            : 
     895                 :          0 :       tk::real rho(0.0);
     896         [ -  - ]:          0 :       for (std::size_t k=0; k<nmat; ++k)
     897                 :          0 :         rho += ur[densityIdx(nmat, k)];
     898                 :            : 
     899                 :            :       // get primitives in boundary state
     900                 :            : 
     901                 :            :       // velocity
     902                 :          0 :       ur[ncomp+velocityIdx(nmat, 0)] = ur[momentumIdx(nmat, 0)] / rho;
     903                 :          0 :       ur[ncomp+velocityIdx(nmat, 1)] = ur[momentumIdx(nmat, 1)] / rho;
     904                 :          0 :       ur[ncomp+velocityIdx(nmat, 2)] = ur[momentumIdx(nmat, 2)] / rho;
     905                 :            : 
     906                 :            :       // material pressures
     907         [ -  - ]:          0 :       for (std::size_t k=0; k<nmat; ++k)
     908                 :            :       {
     909         [ -  - ]:          0 :         auto gk = getDeformGrad(nmat, k, ur);
     910         [ -  - ]:          0 :         ur[ncomp+pressureIdx(nmat, k)] = mat_blk[k].compute< EOS::pressure >(
     911                 :          0 :           ur[densityIdx(nmat, k)], ur[ncomp+velocityIdx(nmat, 0)],
     912                 :          0 :           ur[ncomp+velocityIdx(nmat, 1)], ur[ncomp+velocityIdx(nmat, 2)],
     913                 :          0 :           ur[energyIdx(nmat, k)], ur[volfracIdx(nmat, k)], k, gk );
     914                 :            :       }
     915                 :            : 
     916 [ -  - ][ -  - ]:          0 :       Assert( ur.size() == ncomp+nmat+3, "Incorrect size for appended "
         [ -  - ][ -  - ]
     917                 :            :               "boundary state vector" );
     918                 :            : 
     919         [ -  - ]:          0 :       return {{ std::move(ul), std::move(ur) }};
     920                 :            :     }
     921                 :            : 
     922                 :            :     // Other boundary condition types that do not depend on "Problem" should be
     923                 :            :     // added in BCFunctions.hpp
     924                 :            : };
     925                 :            : 
     926                 :            : } // fv::
     927                 :            : 
     928                 :            : } // inciter::
     929                 :            : 
     930                 :            : #endif // FVMultiMat_h

Generated by: LCOV version 1.14