Quinoa all test code coverage report
Current view: top level - PDE/MultiMat - DGMultiMat.hpp (source / functions) Hit Total Coverage
Commit: -128-NOTFOUND Lines: 281 558 50.4 %
Date: 2025-02-13 15:54:38 Functions: 64 720 8.9 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 212 826 25.7 %

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/MultiMat/DGMultiMat.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 discontinuous Galerkin
       9                 :            :     finite elements
      10                 :            :   \details   This file implements calls to the physics operators governing
      11                 :            :     compressible multi-material flow (with velocity equilibrium) using
      12                 :            :     discontinuous Galerkin discretizations.
      13                 :            : */
      14                 :            : // *****************************************************************************
      15                 :            : #ifndef DGMultiMat_h
      16                 :            : #define DGMultiMat_h
      17                 :            : 
      18                 :            : #include <cmath>
      19                 :            : #include <algorithm>
      20                 :            : #include <unordered_set>
      21                 :            : #include <map>
      22                 :            : #include <array>
      23                 :            : 
      24                 :            : #include "Macro.hpp"
      25                 :            : #include "Exception.hpp"
      26                 :            : #include "Vector.hpp"
      27                 :            : #include "ContainerUtil.hpp"
      28                 :            : #include "UnsMesh.hpp"
      29                 :            : #include "Inciter/InputDeck/InputDeck.hpp"
      30                 :            : #include "Integrate/Basis.hpp"
      31                 :            : #include "Integrate/Quadrature.hpp"
      32                 :            : #include "Integrate/Initialize.hpp"
      33                 :            : #include "Integrate/Mass.hpp"
      34                 :            : #include "Integrate/Surface.hpp"
      35                 :            : #include "Integrate/Boundary.hpp"
      36                 :            : #include "Integrate/Volume.hpp"
      37                 :            : #include "Integrate/MultiMatTerms.hpp"
      38                 :            : #include "Integrate/Source.hpp"
      39                 :            : #include "Integrate/SolidTerms.hpp"
      40                 :            : #include "RiemannChoice.hpp"
      41                 :            : #include "MultiMat/MultiMatIndexing.hpp"
      42                 :            : #include "Reconstruction.hpp"
      43                 :            : #include "Limiter.hpp"
      44                 :            : #include "Problem/FieldOutput.hpp"
      45                 :            : #include "Problem/BoxInitialization.hpp"
      46                 :            : #include "PrefIndicator.hpp"
      47                 :            : #include "MultiMat/BCFunctions.hpp"
      48                 :            : #include "MultiMat/MiscMultiMatFns.hpp"
      49                 :            : #include "EoS/GetMatProp.hpp"
      50                 :            : 
      51                 :            : namespace inciter {
      52                 :            : 
      53                 :            : extern ctr::InputDeck g_inputdeck;
      54                 :            : 
      55                 :            : namespace dg {
      56                 :            : 
      57                 :            : //! \brief MultiMat used polymorphically with tk::DGPDE
      58                 :            : //! \details The template arguments specify policies and are used to configure
      59                 :            : //!   the behavior of the class. The policies are:
      60                 :            : //!   - Physics - physics configuration, see PDE/MultiMat/Physics.h
      61                 :            : //!   - Problem - problem configuration, see PDE/MultiMat/Problem.h
      62                 :            : //! \note The default physics is Euler, set in inciter::deck::check_multimat()
      63                 :            : template< class Physics, class Problem >
      64                 :            : class MultiMat {
      65                 :            : 
      66                 :            :   private:
      67                 :            :     using eq = tag::multimat;
      68                 :            : 
      69                 :            :   public:
      70                 :            :     //! Constructor
      71                 :        118 :     explicit MultiMat() :
      72                 :        118 :       m_ncomp( g_inputdeck.get< tag::ncomp >() ),
      73                 :            :       m_nprim(nprim()),
      74                 :            :       m_riemann( multimatRiemannSolver(
      75                 :        118 :         g_inputdeck.get< tag::flux >() ) )
      76                 :            :     {
      77                 :            :       // associate boundary condition configurations with state functions
      78 [ +  - ][ +  - ]:       1534 :       brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  + ]
         [ +  + ][ -  - ]
                 [ -  - ]
      79                 :            :       // BC State functions
      80                 :            :       { dirichlet
      81                 :            :       , symmetry
      82                 :            :       , invalidBC         // Outlet BC not implemented
      83                 :            :       , farfield
      84                 :            :       , extrapolate
      85                 :            :       , noslipwall },
      86                 :            :       // BC Gradient functions
      87                 :            :       { noOpGrad
      88                 :            :       , symmetryGrad
      89                 :            :       , noOpGrad
      90                 :            :       , noOpGrad
      91                 :            :       , noOpGrad
      92                 :            :       , noOpGrad }
      93                 :            :       ) );
      94                 :            : 
      95                 :            :       // Inlet BC has a different structure than above BCs, so it must be 
      96                 :            :       // handled differently than with ConfigBC
      97 [ +  - ][ +  - ]:        118 :       ConfigInletBC(m_bc, inlet, zeroGrad);
                 [ +  - ]
      98                 :            : 
      99                 :            :       // EoS initialization
     100         [ +  - ]:        118 :       initializeMaterialEoS( m_mat_blk );
     101                 :        118 :     }
     102                 :            : 
     103                 :            :     //! Find the number of primitive quantities required for this PDE system
     104                 :            :     //! \return The number of primitive quantities required to be stored for
     105                 :            :     //!   this PDE system
     106                 :        183 :     std::size_t nprim() const
     107                 :            :     {
     108                 :        183 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     109                 :            :       const auto& solidx = inciter::g_inputdeck.get<
     110                 :        183 :         tag::matidxmap, tag::solidx >();
     111                 :            : 
     112                 :            :       // individual material pressures and three velocity components
     113                 :        183 :       std::size_t np(nmat+3);
     114                 :            : 
     115         [ +  + ]:        633 :       for (std::size_t k=0; k<nmat; ++k) {
     116         [ -  + ]:        450 :         if (solidx[k] > 0) {
     117                 :            :           // individual material Cauchy stress tensor components
     118                 :          0 :           np += 6;
     119                 :            :         }
     120                 :            :       }
     121                 :            : 
     122                 :        183 :       return np;
     123                 :            :     }
     124                 :            : 
     125                 :            :     //! Find the number of materials set up for this PDE system
     126                 :            :     //! \return The number of materials set up for this PDE system
     127                 :          0 :     std::size_t nmat() const
     128                 :            :     {
     129                 :          0 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     130                 :          0 :       return nmat;
     131                 :            :     }
     132                 :            : 
     133                 :            :     //! Assign number of DOFs per equation in the PDE system
     134                 :            :     //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
     135                 :            :     //!   equation
     136                 :         65 :     void numEquationDofs(std::vector< std::size_t >& numEqDof) const
     137                 :            :     {
     138                 :            :       // all equation-dofs initialized to ndofs first
     139         [ +  + ]:        665 :       for (std::size_t i=0; i<m_ncomp; ++i) {
     140                 :        600 :         numEqDof.push_back(g_inputdeck.get< tag::ndof >());
     141                 :            :       }
     142                 :            : 
     143                 :            :       // volume fractions are P0Pm (ndof = 1) for multi-material simulations
     144                 :         65 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     145         [ +  - ]:         65 :       if(nmat > 1)
     146         [ +  + ]:        200 :         for (std::size_t k=0; k<nmat; ++k)
     147                 :        135 :           numEqDof[volfracIdx(nmat, k)] = 1;
     148                 :         65 :     }
     149                 :            : 
     150                 :            :     //! Determine elements that lie inside the user-defined IC box
     151                 :            :     //! \param[in] geoElem Element geometry array
     152                 :            :     //! \param[in] nielem Number of internal elements
     153                 :            :     //! \param[in,out] inbox List of nodes at which box user ICs are set for
     154                 :            :     //!    each IC box
     155                 :         65 :     void IcBoxElems( const tk::Fields& geoElem,
     156                 :            :       std::size_t nielem,
     157                 :            :       std::vector< std::unordered_set< std::size_t > >& inbox ) const
     158                 :            :     {
     159                 :         65 :       tk::BoxElems< eq >(geoElem, nielem, inbox);
     160                 :         65 :     }
     161                 :            : 
     162                 :            :     //! Find how many 'stiff equations', which are the inverse
     163                 :            :     //! deformation equations because of plasticity
     164                 :            :     //! \return number of stiff equations
     165                 :        390 :     std::size_t nstiffeq() const
     166                 :            :     {
     167                 :        390 :       const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
     168                 :        390 :       std::size_t nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     169                 :        390 :       return 9*numSolids(nmat, solidx);
     170                 :            :     }
     171                 :            : 
     172                 :            :     //! Find how many 'non-stiff equations', which are the inverse
     173                 :            :     //! deformation equations because of plasticity
     174                 :            :     //! \return number of stiff equations
     175                 :        130 :     std::size_t nnonstiffeq() const
     176                 :            :     {
     177                 :        130 :       return m_ncomp-nstiffeq();
     178                 :            :     }
     179                 :            : 
     180                 :            :     //! Locate the stiff equations.
     181                 :            :     //! \param[out] stiffEqIdx list with pointers to stiff equations
     182                 :          0 :     void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const
     183                 :            :     {
     184 [ -  - ][ -  - ]:          0 :       stiffEqIdx.resize(nstiffeq(), 0);
     185                 :          0 :       const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
     186                 :          0 :       std::size_t nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     187                 :          0 :       std::size_t icnt = 0;
     188         [ -  - ]:          0 :       for (std::size_t k=0; k<nmat; ++k)
     189         [ -  - ]:          0 :         if (solidx[k] > 0)
     190         [ -  - ]:          0 :           for (std::size_t i=0; i<3; ++i)
     191         [ -  - ]:          0 :             for (std::size_t j=0; j<3; ++j)
     192                 :            :             {
     193                 :          0 :               stiffEqIdx[icnt] =
     194                 :          0 :                 inciter::deformIdx(nmat, solidx[k], i, j);
     195                 :          0 :               icnt++;
     196                 :            :             }
     197                 :          0 :     }
     198                 :            : 
     199                 :            :     //! Locate the nonstiff equations.
     200                 :            :     //! \param[out] nonStiffEqIdx list with pointers to nonstiff equations
     201                 :          0 :     void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const
     202                 :            :     {
     203 [ -  - ][ -  - ]:          0 :       nonStiffEqIdx.resize(nnonstiffeq(), 0);
     204         [ -  - ]:          0 :       for (std::size_t icomp=0; icomp<nnonstiffeq(); icomp++)
     205                 :          0 :         nonStiffEqIdx[icomp] = icomp;
     206                 :          0 :     }
     207                 :            : 
     208                 :            :     //! Initialize the compressible flow equations, prepare for time integration
     209                 :            :     //! \param[in] L Block diagonal mass matrix
     210                 :            :     //! \param[in] inpoel Element-node connectivity
     211                 :            :     //! \param[in] coord Array of nodal coordinates
     212                 :            :     //! \param[in] inbox List of elements at which box user ICs are set for
     213                 :            :     //!   each IC box
     214                 :            :     //! \param[in] elemblkid Element ids associated with mesh block ids where
     215                 :            :     //!   user ICs are set
     216                 :            :     //! \param[in,out] unk Array of unknowns
     217                 :            :     //! \param[in] t Physical time
     218                 :            :     //! \param[in] nielem Number of internal elements
     219                 :         65 :     void initialize( const tk::Fields& L,
     220                 :            :       const std::vector< std::size_t >& inpoel,
     221                 :            :       const tk::UnsMesh::Coords& coord,
     222                 :            :       const std::vector< std::unordered_set< std::size_t > >& inbox,
     223                 :            :       const std::unordered_map< std::size_t, std::set< std::size_t > >&
     224                 :            :         elemblkid,
     225                 :            :       tk::Fields& unk,
     226                 :            :       tk::real t,
     227                 :            :       const std::size_t nielem ) const
     228                 :            :     {
     229 [ +  - ][ +  - ]:         65 :       tk::initialize( m_ncomp, m_mat_blk, L, inpoel, coord,
     230                 :            :                       Problem::initialize, unk, t, nielem );
     231                 :            : 
     232                 :         65 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     233                 :         65 :       const auto& ic = g_inputdeck.get< tag::ic >();
     234                 :         65 :       const auto& icbox = ic.get< tag::box >();
     235                 :         65 :       const auto& icmbk = ic.get< tag::meshblock >();
     236                 :            : 
     237                 :         65 :       const auto& bgpre = ic.get< tag::pressure >();
     238                 :         65 :       const auto& bgtemp = ic.get< tag::temperature >();
     239                 :            : 
     240                 :            :       // Set initial conditions inside user-defined IC boxes and mesh blocks
     241         [ +  - ]:        130 :       std::vector< tk::real > s(m_ncomp, 0.0);
     242         [ +  + ]:      22511 :       for (std::size_t e=0; e<nielem; ++e) {
     243                 :            :         // inside user-defined box
     244         [ -  + ]:      22446 :         if (!icbox.empty()) {
     245                 :          0 :           std::size_t bcnt = 0;
     246         [ -  - ]:          0 :           for (const auto& b : icbox) {   // for all boxes
     247 [ -  - ][ -  - ]:          0 :             if (inbox.size() > bcnt && inbox[bcnt].find(e) != inbox[bcnt].end())
         [ -  - ][ -  - ]
     248                 :            :             {
     249         [ -  - ]:          0 :               std::vector< tk::real > box
     250                 :          0 :                 { b.template get< tag::xmin >(), b.template get< tag::xmax >(),
     251                 :          0 :                   b.template get< tag::ymin >(), b.template get< tag::ymax >(),
     252                 :          0 :                   b.template get< tag::zmin >(), b.template get< tag::zmax >() };
     253                 :          0 :               auto V_ex = (box[1]-box[0]) * (box[3]-box[2]) * (box[5]-box[4]);
     254         [ -  - ]:          0 :               for (std::size_t c=0; c<m_ncomp; ++c) {
     255                 :          0 :                 auto mark = c*rdof;
     256         [ -  - ]:          0 :                 s[c] = unk(e,mark);
     257                 :            :                 // set high-order DOFs to zero
     258         [ -  - ]:          0 :                 for (std::size_t i=1; i<rdof; ++i)
     259         [ -  - ]:          0 :                   unk(e,mark+i) = 0.0;
     260                 :            :               }
     261         [ -  - ]:          0 :               initializeBox<ctr::boxList>( m_mat_blk, V_ex, t, b, bgpre,
     262                 :            :                 bgtemp, s );
     263                 :            :               // store box-initialization in solution vector
     264         [ -  - ]:          0 :               for (std::size_t c=0; c<m_ncomp; ++c) {
     265                 :          0 :                 auto mark = c*rdof;
     266         [ -  - ]:          0 :                 unk(e,mark) = s[c];
     267                 :            :               }
     268                 :            :             }
     269                 :          0 :             ++bcnt;
     270                 :            :           }
     271                 :            :         }
     272                 :            : 
     273                 :            :         // inside user-specified mesh blocks
     274         [ -  + ]:      22446 :         if (!icmbk.empty()) {
     275         [ -  - ]:          0 :           for (const auto& b : icmbk) { // for all blocks
     276                 :          0 :             auto blid = b.get< tag::blockid >();
     277                 :          0 :             auto V_ex = b.get< tag::volume >();
     278 [ -  - ][ -  - ]:          0 :             if (elemblkid.find(blid) != elemblkid.end()) {
     279         [ -  - ]:          0 :               const auto& elset = tk::cref_find(elemblkid, blid);
     280 [ -  - ][ -  - ]:          0 :               if (elset.find(e) != elset.end()) {
     281         [ -  - ]:          0 :                 initializeBox<ctr::meshblockList>( m_mat_blk, V_ex, t, b,
     282                 :            :                   bgpre, bgtemp, s );
     283                 :            :                 // store initialization in solution vector
     284         [ -  - ]:          0 :                 for (std::size_t c=0; c<m_ncomp; ++c) {
     285                 :          0 :                   auto mark = c*rdof;
     286         [ -  - ]:          0 :                   unk(e,mark) = s[c];
     287                 :            :                 }
     288                 :            :               }
     289                 :            :             }
     290                 :            :           }
     291                 :            :         }
     292                 :            :       }
     293                 :         65 :     }
     294                 :            : 
     295                 :            :     //! Compute density constraint for a given material
     296                 :            :     //! \param[in] nelem Number of elements
     297                 :            :     //! \param[in] unk Array of unknowns
     298                 :            :     //! \param[out] densityConstr Density Constraint: rho/(rho0*det(g))
     299                 :        155 :     void computeDensityConstr( std::size_t nelem,
     300                 :            :                                tk::Fields& unk,
     301                 :            :                                std::vector< tk::real >& densityConstr) const
     302                 :            :     {
     303                 :        155 :       const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
     304                 :        155 :       std::size_t rdof = g_inputdeck.get< tag::rdof >();
     305                 :        155 :       std::size_t nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     306         [ +  + ]:      64461 :       for (std::size_t e=0; e<nelem; ++e)
     307                 :      64306 :         densityConstr[e] = 0.0;
     308         [ +  + ]:        480 :       for (std::size_t imat=0; imat<nmat; ++imat)
     309         [ -  + ]:        325 :         if (solidx[imat] > 0)
     310                 :            :         {
     311         [ -  - ]:          0 :           for (std::size_t e=0; e<nelem; ++e)
     312                 :            :           {
     313                 :            :             // Retrieve unknowns
     314         [ -  - ]:          0 :             tk::real arho = unk(e, densityDofIdx(nmat, imat, rdof, 0));
     315                 :            :             std::array< std::array< tk::real, 3 >, 3 > g;
     316         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
     317         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
     318         [ -  - ]:          0 :                 g[i][j] = unk(e, deformDofIdx(nmat, solidx[imat], i, j, rdof, 0));
     319                 :            :             // Compute determinant of g
     320                 :          0 :             tk::real detg = tk::determinant(g);
     321                 :            :             // Compute constraint measure
     322         [ -  - ]:          0 :             densityConstr[e] += arho/(m_mat_blk[imat].compute< EOS::rho0 >()*detg);
     323                 :            :           }
     324                 :            :         }
     325                 :            :         else
     326                 :            :         {
     327         [ +  + ]:     150795 :           for (std::size_t e=0; e<nelem; ++e)
     328                 :            :           {
     329                 :            :             // Retrieve alpha and add it to the constraint measure
     330                 :     150470 :             tk::real alpha = unk(e, volfracDofIdx(nmat, imat, rdof, 0));
     331                 :     150470 :             densityConstr[e] += alpha;
     332                 :            :           }
     333                 :            :         }
     334                 :        155 :     }
     335                 :            : 
     336                 :            :     //! Compute the left hand side block-diagonal mass matrix
     337                 :            :     //! \param[in] geoElem Element geometry array
     338                 :            :     //! \param[in,out] l Block diagonal mass matrix
     339                 :         65 :     void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
     340                 :         65 :       const auto ndof = g_inputdeck.get< tag::ndof >();
     341                 :            :       // Unlike Compflow and Transport, there is a weak reconstruction about
     342                 :            :       // conservative variable after limiting function which will require the
     343                 :            :       // size of left hand side vector to be rdof
     344                 :         65 :       tk::mass( m_ncomp, ndof, geoElem, l );
     345                 :         65 :     }
     346                 :            : 
     347                 :            :     //! Update the interface cells to first order dofs
     348                 :            :     //! \param[in] unk Array of unknowns
     349                 :            :     //! \param[in] nielem Number of internal elements
     350                 :            :     //! \param[in,out] ndofel Array of dofs
     351                 :            :     //! \param[in,out] interface Vector of interface marker
     352                 :            :     //! \details This function resets the high-order terms in interface cells.
     353                 :       5205 :     void updateInterfaceCells( tk::Fields& unk,
     354                 :            :       std::size_t nielem,
     355                 :            :       std::vector< std::size_t >& ndofel,
     356                 :            :       std::vector< std::size_t >& interface ) const
     357                 :            :     {
     358                 :       5205 :       auto intsharp =
     359                 :       5205 :         g_inputdeck.get< tag::multimat, tag::intsharp >();
     360                 :            :       // If this cell is not material interface, return this function
     361         [ +  + ]:       5205 :       if(not intsharp)  return;
     362                 :            : 
     363                 :        375 :       auto rdof = g_inputdeck.get< tag::rdof >();
     364                 :        375 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     365                 :            :       const auto& solidx = g_inputdeck.get<
     366                 :        375 :         tag::matidxmap, tag::solidx >();
     367                 :            : 
     368         [ +  + ]:     227775 :       for (std::size_t e=0; e<nielem; ++e) {
     369         [ +  - ]:     454800 :         std::vector< std::size_t > matInt(nmat, 0);
     370         [ +  - ]:     454800 :         std::vector< tk::real > alAvg(nmat, 0.0);
     371         [ +  + ]:     682200 :         for (std::size_t k=0; k<nmat; ++k)
     372         [ +  - ]:     454800 :           alAvg[k] = unk(e, volfracDofIdx(nmat,k,rdof,0));
     373         [ +  - ]:     227400 :         auto intInd = interfaceIndicator(nmat, alAvg, matInt);
     374                 :            : 
     375                 :            :         // interface cells cannot be high-order
     376         [ +  + ]:     227400 :         if (intInd) {
     377                 :      10834 :           interface[e] = 1;
     378         [ +  + ]:      32502 :           for (std::size_t k=0; k<nmat; ++k) {
     379         [ +  - ]:      21668 :             if (matInt[k]) {
     380         [ +  + ]:      86672 :               for (std::size_t i=1; i<rdof; ++i) {
     381         [ +  - ]:      65004 :                 unk(e, densityDofIdx(nmat,k,rdof,i)) = 0.0;
     382         [ +  - ]:      65004 :                 unk(e, energyDofIdx(nmat,k,rdof,i)) = 0.0;
     383                 :            :               }
     384         [ -  + ]:      21668 :               if (solidx[k] > 0) {
     385         [ -  - ]:          0 :                 for (std::size_t i=0; i<3; ++i)
     386         [ -  - ]:          0 :                   for (std::size_t j=0; j<3; ++j)
     387         [ -  - ]:          0 :                     for (std::size_t idof=1; idof<rdof; ++idof) {
     388         [ -  - ]:          0 :                       unk(e, deformDofIdx(nmat,solidx[k],i,j,rdof,idof)) = 0.0;
     389                 :            :                     }
     390                 :            :               }
     391                 :            :             }
     392                 :            :           }
     393         [ +  + ]:      43336 :           for (std::size_t idir=0; idir<3; ++idir) {
     394         [ +  + ]:     130008 :             for (std::size_t i=1; i<rdof; ++i) {
     395         [ +  - ]:      97506 :               unk(e, momentumDofIdx(nmat,idir,rdof,i)) = 0.0;
     396                 :            :             }
     397                 :            :           }
     398                 :            :         } else {
     399                 :            :           // If the cell is marked as interface cell in the previous timestep
     400                 :            :           // and does not marked as interface for the current timestep, DGP2
     401                 :            :           // will be applied for the current timestep in p-adaptive process
     402                 :            :           // Please note this block is added since the spectral decay indicator
     403                 :            :           // does not applied to P0 cells.
     404         [ +  + ]:     216566 :           if (interface[e] == 1) {
     405 [ +  - ][ -  + ]:          6 :             if(ndofel[e] < 10 && rdof == 10) {
                 [ -  + ]
     406                 :          0 :               ndofel[e] = 10;
     407         [ -  - ]:          0 :               for (std::size_t k=0; k<nmat; ++k) {
     408         [ -  - ]:          0 :                 for (std::size_t i=1; i<rdof; ++i) {
     409         [ -  - ]:          0 :                   unk(e, densityDofIdx(nmat,k,rdof,i)) = 0.0;
     410         [ -  - ]:          0 :                   unk(e, energyDofIdx(nmat,k,rdof,i)) = 0.0;
     411                 :            :                 }
     412                 :            :               }
     413         [ -  - ]:          0 :               for (std::size_t idir=0; idir<3; ++idir) {
     414         [ -  - ]:          0 :                 for (std::size_t i=1; i<rdof; ++i) {
     415         [ -  - ]:          0 :                   unk(e, momentumDofIdx(nmat,idir,rdof,i)) = 0.0;
     416                 :            :                 }
     417                 :            :               }
     418                 :            :             }
     419                 :            :           }
     420                 :     216566 :           interface[e] = 0;
     421                 :            :         }
     422                 :            :       }
     423                 :            :     }
     424                 :            : 
     425                 :            :     //! Update the primitives for this PDE system
     426                 :            :     //! \param[in] unk Array of unknowns
     427                 :            :     //! \param[in] L The left hand side block-diagonal mass matrix
     428                 :            :     //! \param[in] geoElem Element geometry array
     429                 :            :     //! \param[in,out] prim Array of primitives
     430                 :            :     //! \param[in] nielem Number of internal elements
     431                 :            :     //! \param[in] ndofel Array of dofs
     432                 :            :     //! \details This function computes and stores the dofs for primitive
     433                 :            :     //!   quantities, which are required for obtaining reconstructed states used
     434                 :            :     //!   in the Riemann solver. See /PDE/Riemann/AUSM.hpp, where the
     435                 :            :     //!   normal velocity for advection is calculated from independently
     436                 :            :     //!   reconstructed velocities.
     437                 :       5270 :     void updatePrimitives( const tk::Fields& unk,
     438                 :            :                            const tk::Fields& L,
     439                 :            :                            const tk::Fields& geoElem,
     440                 :            :                            tk::Fields& prim,
     441                 :            :                            std::size_t nielem,
     442                 :            :                            std::vector< std::size_t >& ndofel ) const
     443                 :            :     {
     444                 :       5270 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     445                 :       5270 :       const auto ndof = g_inputdeck.get< tag::ndof >();
     446                 :       5270 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     447                 :       5270 :       const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
     448                 :            : 
     449 [ -  + ][ -  - ]:       5270 :       Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     450                 :            :               "vector and primitive vector at recent time step incorrect" );
     451 [ -  + ][ -  - ]:       5270 :       Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
     452                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     453 [ -  + ][ -  - ]:       5270 :       Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
         [ -  - ][ -  - ]
                 [ -  - ]
     454                 :            :               "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
     455                 :            : 
     456         [ +  + ]:    2189396 :       for (std::size_t e=0; e<nielem; ++e)
     457                 :            :       {
     458         [ +  - ]:    4368252 :         std::vector< tk::real > R(m_nprim*ndof, 0.0);
     459                 :            : 
     460         [ +  - ]:    2184126 :         auto ng = tk::NGvol(ndof);
     461                 :            : 
     462                 :            :         // arrays for quadrature points
     463                 :    4368252 :         std::array< std::vector< tk::real >, 3 > coordgp;
     464                 :    4368252 :         std::vector< tk::real > wgp;
     465                 :            : 
     466         [ +  - ]:    2184126 :         coordgp[0].resize( ng );
     467         [ +  - ]:    2184126 :         coordgp[1].resize( ng );
     468         [ +  - ]:    2184126 :         coordgp[2].resize( ng );
     469         [ +  - ]:    2184126 :         wgp.resize( ng );
     470                 :            : 
     471         [ +  - ]:    2184126 :         tk::GaussQuadratureTet( ng, coordgp, wgp );
     472                 :            : 
     473                 :            :         // Local degree of freedom
     474                 :    2184126 :         auto dof_el = ndofel[e];
     475                 :            : 
     476                 :            :         // Loop over quadrature points in element e
     477         [ +  + ]:    6399692 :         for (std::size_t igp=0; igp<ng; ++igp)
     478                 :            :         {
     479                 :            :           // Compute the basis function
     480         [ +  - ]:    8431132 :           auto B =
     481                 :    4215566 :             tk::eval_basis( dof_el, coordgp[0][igp], coordgp[1][igp], coordgp[2][igp] );
     482                 :            : 
     483         [ +  - ]:    4215566 :           auto w = wgp[igp] * geoElem(e, 0);
     484                 :            : 
     485         [ +  - ]:    8431132 :           auto state = tk::eval_state( m_ncomp, rdof, dof_el, e, unk, B );
     486                 :            : 
     487                 :            :           // bulk density at quadrature point
     488                 :    4215566 :           tk::real rhob(0.0);
     489         [ +  + ]:   13746884 :           for (std::size_t k=0; k<nmat; ++k)
     490                 :    9531318 :             rhob += state[densityIdx(nmat, k)];
     491                 :            : 
     492                 :            :           // velocity vector at quadrature point
     493                 :            :           std::array< tk::real, 3 >
     494                 :    4215566 :             vel{ state[momentumIdx(nmat, 0)]/rhob,
     495                 :    4215566 :                  state[momentumIdx(nmat, 1)]/rhob,
     496                 :    4215566 :                  state[momentumIdx(nmat, 2)]/rhob };
     497                 :            : 
     498         [ +  - ]:    8431132 :           std::vector< tk::real > pri(m_nprim, 0.0);
     499                 :            : 
     500                 :            :           // Evaluate material pressure at quadrature point
     501         [ +  + ]:   13746884 :           for(std::size_t imat = 0; imat < nmat; imat++)
     502                 :            :           {
     503                 :    9531318 :             auto alphamat = state[volfracIdx(nmat, imat)];
     504                 :    9531318 :             auto arhomat = state[densityIdx(nmat, imat)];
     505                 :    9531318 :             auto arhoemat = state[energyIdx(nmat, imat)];
     506         [ +  - ]:    9531318 :             auto gmat = getDeformGrad(nmat, imat, state);
     507                 :    9531318 :             pri[pressureIdx(nmat,imat)] = m_mat_blk[imat].compute<
     508         [ +  - ]:    9531318 :               EOS::pressure >( arhomat, vel[0], vel[1], vel[2], arhoemat,
     509                 :            :               alphamat, imat, gmat );
     510                 :            : 
     511         [ +  - ]:    9531318 :             pri[pressureIdx(nmat,imat)] = constrain_pressure( m_mat_blk,
     512                 :            :               pri[pressureIdx(nmat,imat)], arhomat, alphamat, imat);
     513                 :            : 
     514         [ -  + ]:    9531318 :             if (solidx[imat] > 0) {
     515         [ -  - ]:          0 :               auto asigmat = m_mat_blk[imat].computeTensor< EOS::CauchyStress >(
     516                 :            :               arhomat, vel[0], vel[1], vel[2], arhoemat,
     517                 :            :               alphamat, imat, gmat );
     518                 :            : 
     519                 :          0 :               pri[stressIdx(nmat,solidx[imat],0)] = asigmat[0][0];
     520                 :          0 :               pri[stressIdx(nmat,solidx[imat],1)] = asigmat[1][1];
     521                 :          0 :               pri[stressIdx(nmat,solidx[imat],2)] = asigmat[2][2];
     522                 :          0 :               pri[stressIdx(nmat,solidx[imat],3)] = asigmat[0][1];
     523                 :          0 :               pri[stressIdx(nmat,solidx[imat],4)] = asigmat[0][2];
     524                 :          0 :               pri[stressIdx(nmat,solidx[imat],5)] = asigmat[1][2];
     525                 :            :             }
     526                 :            :           }
     527                 :            : 
     528                 :            :           // Evaluate bulk velocity at quadrature point
     529         [ +  + ]:   16862264 :           for (std::size_t idir=0; idir<3; ++idir) {
     530                 :   12646698 :             pri[velocityIdx(nmat,idir)] = vel[idir];
     531                 :            :           }
     532                 :            : 
     533         [ +  + ]:   26393582 :           for(std::size_t k = 0; k < m_nprim; k++)
     534                 :            :           {
     535                 :   22178016 :             auto mark = k * ndof;
     536         [ +  + ]:   82445532 :             for(std::size_t idof = 0; idof < dof_el; idof++)
     537                 :   60267516 :               R[mark+idof] += w * pri[k] * B[idof];
     538                 :            :           }
     539                 :            :         }
     540                 :            : 
     541                 :            :         // Update the DG solution of primitive variables
     542         [ +  + ]:   14204942 :         for(std::size_t k = 0; k < m_nprim; k++)
     543                 :            :         {
     544                 :   12020816 :           auto mark = k * ndof;
     545                 :   12020816 :           auto rmark = k * rdof;
     546         [ +  + ]:   31659532 :           for(std::size_t idof = 0; idof < dof_el; idof++)
     547                 :            :           {
     548 [ +  - ][ +  - ]:   19638716 :             prim(e, rmark+idof) = R[mark+idof] / L(e, mark+idof);
     549 [ +  - ][ +  + ]:   19638716 :             if(fabs(prim(e, rmark+idof)) < 1e-16)
     550         [ +  - ]:    5168733 :               prim(e, rmark+idof) = 0;
     551                 :            :           }
     552                 :            :         }
     553                 :            :       }
     554                 :       5270 :     }
     555                 :            : 
     556                 :            :     //! Clean up the state of trace materials for this PDE system
     557                 :            :     //! \param[in] t Physical time
     558                 :            :     //! \param[in] geoElem Element geometry array
     559                 :            :     //! \param[in,out] unk Array of unknowns
     560                 :            :     //! \param[in,out] prim Array of primitives
     561                 :            :     //! \param[in] nielem Number of internal elements
     562                 :            :     //! \details This function cleans up the state of materials present in trace
     563                 :            :     //!   quantities in each cell. Specifically, the state of materials with
     564                 :            :     //!   very low volume-fractions in a cell is replaced by the state of the
     565                 :            :     //!   material which is present in the largest quantity in that cell. This
     566                 :            :     //!   becomes necessary when shocks pass through cells which contain a very
     567                 :            :     //!   small amount of material. The state of that tiny material might
     568                 :            :     //!   become unphysical and cause solution to diverge; thus requiring such
     569                 :            :     //!   a "reset".
     570                 :       5205 :     void cleanTraceMaterial( tk::real t,
     571                 :            :                              const tk::Fields& geoElem,
     572                 :            :                              tk::Fields& unk,
     573                 :            :                              tk::Fields& prim,
     574                 :            :                              std::size_t nielem ) const
     575                 :            :     {
     576                 :       5205 :       [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
     577                 :       5205 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     578                 :            : 
     579 [ -  + ][ -  - ]:       5205 :       Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     580                 :            :               "vector and primitive vector at recent time step incorrect" );
     581 [ -  + ][ -  - ]:       5205 :       Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
     582                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     583 [ -  + ][ -  - ]:       5205 :       Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
         [ -  - ][ -  - ]
                 [ -  - ]
     584                 :            :               "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
     585                 :            : 
     586                 :       5205 :       auto neg_density = cleanTraceMultiMat(t, nielem, m_mat_blk, geoElem, nmat,
     587                 :            :         unk, prim);
     588                 :            : 
     589 [ -  + ][ -  - ]:       5205 :       if (neg_density) Throw("Negative partial density.");
         [ -  - ][ -  - ]
     590                 :       5205 :     }
     591                 :            : 
     592                 :            :     //! Reconstruct second-order solution from first-order
     593                 :            :     //! \param[in] geoElem Element geometry array
     594                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     595                 :            :     //! \param[in] esup Elements-surrounding-nodes connectivity
     596                 :            :     //! \param[in] inpoel Element-node connectivity
     597                 :            :     //! \param[in] coord Array of nodal coordinates
     598                 :            :     //! \param[in,out] U Solution vector at recent time step
     599                 :            :     //! \param[in,out] P Vector of primitives at recent time step
     600                 :            :     //! \param[in] pref Indicator for p-adaptive algorithm
     601                 :            :     //! \param[in] ndofel Vector of local number of degrees of freedome
     602                 :       4080 :     void reconstruct( tk::real,
     603                 :            :                       const tk::Fields&,
     604                 :            :                       const tk::Fields& geoElem,
     605                 :            :                       const inciter::FaceData& fd,
     606                 :            :                       const std::map< std::size_t, std::vector< std::size_t > >&
     607                 :            :                         esup,
     608                 :            :                       const std::vector< std::size_t >& inpoel,
     609                 :            :                       const tk::UnsMesh::Coords& coord,
     610                 :            :                       tk::Fields& U,
     611                 :            :                       tk::Fields& P,
     612                 :            :                       const bool pref,
     613                 :            :                       const std::vector< std::size_t >& ndofel ) const
     614                 :            :     {
     615                 :       4080 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     616                 :       4080 :       const auto ndof = g_inputdeck.get< tag::ndof >();
     617                 :            : 
     618                 :       4080 :       bool is_p0p1(false);
     619 [ +  - ][ +  + ]:       4080 :       if (rdof == 4 && ndof == 1)
     620                 :       3300 :         is_p0p1 = true;
     621                 :            : 
     622                 :       4080 :       const auto nelem = fd.Esuel().size()/4;
     623                 :       4080 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     624                 :            : 
     625 [ -  + ][ -  - ]:       4080 :       Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
     626                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     627                 :            : 
     628                 :            :       //----- reconstruction of conserved quantities -----
     629                 :            :       //--------------------------------------------------
     630                 :            : 
     631         [ +  + ]:     845460 :       for (std::size_t e=0; e<nelem; ++e)
     632                 :            :       {
     633                 :            :         // 1. specify how many variables need to be reconstructed
     634                 :    1682760 :         std::vector< std::size_t > vars;
     635                 :            :         // for p-adaptive DG
     636         [ -  + ]:     841380 :         if (pref) {
     637                 :            :           // If DG is applied, reconstruct only volume fractions
     638         [ -  - ]:          0 :           if(ndofel[e] > 1) {
     639 [ -  - ][ -  - ]:          0 :             for (std::size_t k=0; k<nmat; ++k) vars.push_back(volfracIdx(nmat, k));
     640                 :            :           }
     641                 :            :           else  // If P0P1 is applied for this element
     642 [ -  - ][ -  - ]:          0 :             for (std::size_t c=0; c<m_ncomp; ++c) vars.push_back(c);
     643                 :            :         }
     644                 :            :         else {
     645                 :            :           // for P0P1, reconstruct all variables
     646         [ +  + ]:     841380 :           if (is_p0p1)
     647 [ +  + ][ +  - ]:    3411000 :             for (std::size_t c=0; c<m_ncomp; ++c) vars.push_back(c);
     648                 :            :           // for high-order DG, reconstruct only volume fractions
     649         [ +  - ]:     500280 :           else if (ndof > 1)
     650 [ +  + ][ +  - ]:    1500840 :             for (std::size_t k=0; k<nmat; ++k) vars.push_back(volfracIdx(nmat, k));
     651                 :            :         }
     652                 :            : 
     653                 :            :         // 2. solve 3x3 least-squares system
     654                 :            :         // Reconstruct second-order dofs of volume-fractions in Taylor space
     655                 :            :         // using nodal-stencils, for a good interface-normal estimate
     656         [ +  - ]:     841380 :         tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, U, vars );
     657                 :            : 
     658                 :            :         // 3. transform reconstructed derivatives to Dubiner dofs
     659         [ +  - ]:     841380 :         tk::transform_P0P1( rdof, e, inpoel, coord, U, vars );
     660                 :            :       }
     661                 :            : 
     662                 :            :       //----- reconstruction of primitive quantities -----
     663                 :            :       //--------------------------------------------------
     664                 :            :       // For multimat, conserved and primitive quantities are reconstructed
     665                 :            :       // separately.
     666                 :            : 
     667         [ +  + ]:     845460 :       for (std::size_t e=0; e<nelem; ++e)
     668                 :            :       {
     669                 :            :         // There are two conditions that requires the reconstruction of the
     670                 :            :         // primitive variables:
     671                 :            :         //   1. p-adaptive is triggered and P0P1 scheme is applied to specific
     672                 :            :         //      elements
     673                 :            :         //   2. p-adaptive is not triggered and P0P1 scheme is applied to the
     674                 :            :         //      whole computation domain
     675 [ -  + ][ -  - ]:     841380 :         if ((pref && ndofel[e] == 1) || (!pref && is_p0p1)) {
         [ +  - ][ +  + ]
                 [ +  + ]
     676                 :     682200 :           std::vector< std::size_t > vars;
     677 [ +  + ][ +  - ]:    2046600 :           for (std::size_t c=0; c<m_nprim; ++c) vars.push_back(c);
     678                 :            : 
     679                 :            :           // 1.
     680                 :            :           // Reconstruct second-order dofs of volume-fractions in Taylor space
     681                 :            :           // using nodal-stencils, for a good interface-normal estimate
     682         [ +  - ]:     341100 :           tk::recoLeastSqExtStencil( rdof, e, esup, inpoel, geoElem, P, vars );
     683                 :            : 
     684                 :            :           // 2.
     685         [ +  - ]:     341100 :           tk::transform_P0P1(rdof, e, inpoel, coord, P, vars);
     686                 :            :         }
     687                 :            :       }
     688                 :            : 
     689                 :       4080 :     }
     690                 :            : 
     691                 :            :     //! Limit second-order solution, and primitive quantities separately
     692                 :            :     //! \param[in] t Physical time
     693                 :            :     //! \param[in] pref Indicator for p-adaptive algorithm
     694                 :            :     //! \param[in] geoFace Face geometry array
     695                 :            :     //! \param[in] geoElem Element geometry array
     696                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     697                 :            :     //! \param[in] esup Elements-surrounding-nodes connectivity
     698                 :            :     //! \param[in] inpoel Element-node connectivity
     699                 :            :     //! \param[in] coord Array of nodal coordinates
     700                 :            :     //! \param[in] ndofel Vector of local number of degrees of freedome
     701                 :            :     //! \param[in] gid Local->global node id map
     702                 :            :     //! \param[in] bid Local chare-boundary node ids (value) associated to
     703                 :            :     //!   global node ids (key)
     704                 :            :     //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
     705                 :            :     //!   variables
     706                 :            :     //! \param[in] pNodalExtrm Chare-boundary nodal extrema for primitive
     707                 :            :     //!   variables
     708                 :            :     //! \param[in] mtInv Inverse of Taylor mass matrix
     709                 :            :     //! \param[in,out] U Solution vector at recent time step
     710                 :            :     //! \param[in,out] P Vector of primitives at recent time step
     711                 :            :     //! \param[in,out] shockmarker Vector of shock-marker values
     712                 :       4080 :     void limit( [[maybe_unused]] tk::real t,
     713                 :            :                 const bool pref,
     714                 :            :                 const tk::Fields& geoFace,
     715                 :            :                 const tk::Fields& geoElem,
     716                 :            :                 const inciter::FaceData& fd,
     717                 :            :                 const std::map< std::size_t, std::vector< std::size_t > >& esup,
     718                 :            :                 const std::vector< std::size_t >& inpoel,
     719                 :            :                 const tk::UnsMesh::Coords& coord,
     720                 :            :                 const std::vector< std::size_t >& ndofel,
     721                 :            :                 const std::vector< std::size_t >& gid,
     722                 :            :                 const std::unordered_map< std::size_t, std::size_t >& bid,
     723                 :            :                 const std::vector< std::vector<tk::real> >& uNodalExtrm,
     724                 :            :                 const std::vector< std::vector<tk::real> >& pNodalExtrm,
     725                 :            :                 const std::vector< std::vector<tk::real> >& mtInv,
     726                 :            :                 tk::Fields& U,
     727                 :            :                 tk::Fields& P,
     728                 :            :                 std::vector< std::size_t >& shockmarker ) const
     729                 :            :     {
     730 [ -  + ][ -  - ]:       4080 :       Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     731                 :            :               "vector and primitive vector at recent time step incorrect" );
     732                 :            : 
     733                 :       4080 :       const auto limiter = g_inputdeck.get< tag::limiter >();
     734                 :       4080 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     735                 :       4080 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     736                 :            :       const auto& solidx = g_inputdeck.get<
     737                 :       4080 :         tag::matidxmap, tag::solidx >();
     738                 :            : 
     739                 :            :       // limit vectors of conserved and primitive quantities
     740         [ -  + ]:       4080 :       if (limiter == ctr::LimiterType::SUPERBEEP1)
     741                 :            :       {
     742                 :          0 :         SuperbeeMultiMat_P1( fd.Esuel(), inpoel, ndofel,
     743                 :            :           coord, solidx, U, P, nmat );
     744                 :            :       }
     745 [ +  - ][ +  - ]:       4080 :       else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 4)
     746                 :            :       {
     747         [ +  - ]:       8160 :         VertexBasedMultiMat_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
     748                 :       4080 :           m_mat_blk, fd, geoFace, geoElem, coord, flux, solidx, U, P,
     749                 :            :           nmat, shockmarker );
     750                 :            :       }
     751 [ -  - ][ -  - ]:          0 :       else if (limiter == ctr::LimiterType::VERTEXBASEDP1 && rdof == 10)
     752                 :            :       {
     753         [ -  - ]:          0 :         VertexBasedMultiMat_P2( pref, esup, inpoel, ndofel, fd.Esuel().size()/4,
     754                 :          0 :           m_mat_blk, fd, geoFace, geoElem, coord, gid, bid,
     755                 :            :           uNodalExtrm, pNodalExtrm, mtInv, flux, solidx, U, P, nmat,
     756                 :            :           shockmarker );
     757                 :            :       }
     758         [ -  - ]:          0 :       else if (limiter != ctr::LimiterType::NOLIMITER)
     759                 :            :       {
     760 [ -  - ][ -  - ]:          0 :         Throw("Limiter type not configured for multimat.");
                 [ -  - ]
     761                 :            :       }
     762                 :       4080 :     }
     763                 :            : 
     764                 :            :     //! Apply CPL to the conservative variable solution for this PDE system
     765                 :            :     //! \param[in] prim Array of primitive variables
     766                 :            :     //! \param[in] geoElem Element geometry array
     767                 :            :     //! \param[in] inpoel Element-node connectivity
     768                 :            :     //! \param[in] coord Array of nodal coordinates
     769                 :            :     //! \param[in,out] unk Array of conservative variables
     770                 :            :     //! \param[in] nielem Number of internal elements
     771                 :            :     //! \details This function applies CPL to obtain consistent dofs for
     772                 :            :     //!   conservative quantities based on the limited primitive quantities.
     773                 :            :     //!   See Pandare et al. (2023). On the Design of Stable,
     774                 :            :     //!   Consistent, and Conservative High-Order Methods for Multi-Material
     775                 :            :     //!   Hydrodynamics. J Comp Phys, 112313.
     776                 :         30 :     void CPL( const tk::Fields& prim,
     777                 :            :       const tk::Fields& geoElem,
     778                 :            :       const std::vector< std::size_t >& inpoel,
     779                 :            :       const tk::UnsMesh::Coords& coord,
     780                 :            :       tk::Fields& unk,
     781                 :            :       std::size_t nielem ) const
     782                 :            :     {
     783                 :         30 :       [[maybe_unused]] const auto rdof = g_inputdeck.get< tag::rdof >();
     784                 :         30 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     785                 :            : 
     786 [ -  + ][ -  - ]:         30 :       Assert( unk.nunk() == prim.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     787                 :            :               "vector and primitive vector at recent time step incorrect" );
     788 [ -  + ][ -  - ]:         30 :       Assert( unk.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
     789                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     790 [ -  + ][ -  - ]:         30 :       Assert( prim.nprop() == rdof*m_nprim, "Number of components in vector of "
         [ -  - ][ -  - ]
                 [ -  - ]
     791                 :            :               "primitive quantities must equal "+ std::to_string(rdof*m_nprim) );
     792                 :            : 
     793                 :         30 :       correctLimConservMultiMat(nielem, m_mat_blk, nmat, inpoel,
     794                 :            :         coord, geoElem, prim, unk);
     795                 :         30 :     }
     796                 :            : 
     797                 :            :     //! Return cell-average deformation gradient tensor
     798                 :            :     //! \param[in] unk Solution vector at recent time step
     799                 :            :     //! \param[in] nielem Number of internal elements
     800                 :            :     //! \details This function returns the bulk cell-average inverse
     801                 :            :     //!   deformation gradient tensor
     802                 :          0 :     std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
     803                 :            :       const tk::Fields& unk,
     804                 :            :       std::size_t nielem ) const
     805                 :            :     {
     806                 :          0 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     807                 :          0 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     808                 :          0 :       const auto& solidx = g_inputdeck.get< tag::matidxmap, tag::solidx >();
     809                 :            : 
     810                 :          0 :       std::array< std::vector< tk::real >, 9 > gb;
     811 [ -  - ][ -  - ]:          0 :       if (inciter::haveSolid(nmat, solidx)) {
     812         [ -  - ]:          0 :         for (auto& gij : gb)
     813         [ -  - ]:          0 :           gij.resize(nielem, 0.0);
     814         [ -  - ]:          0 :         for (std::size_t e=0; e<nielem; ++e) {
     815         [ -  - ]:          0 :           for (std::size_t k=0; k<nmat; ++k) {
     816         [ -  - ]:          0 :             if (solidx[k] > 0) {
     817         [ -  - ]:          0 :               for (std::size_t i=0; i<3; ++i)
     818         [ -  - ]:          0 :                 for (std::size_t j=0; j<3; ++j)
     819         [ -  - ]:          0 :                   gb[3*i+j][e] += unk(e, volfracDofIdx(nmat,k,rdof,0)) *
     820         [ -  - ]:          0 :                     unk(e,deformDofIdx(nmat,solidx[k],i,j,rdof,0));
     821                 :            :             }
     822                 :            :           }
     823                 :            :         }
     824                 :            :       }
     825                 :            : 
     826                 :          0 :       return gb;
     827                 :            :     }
     828                 :            : 
     829                 :            : 
     830                 :            :     //! Reset the high order solution for p-adaptive scheme
     831                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     832                 :            :     //! \param[in,out] unk Solution vector at recent time step
     833                 :            :     //! \param[in,out] prim Primitive vector at recent time step
     834                 :            :     //! \param[in] ndofel Vector of local number of degrees of freedome
     835                 :            :     //! \details This function reset the high order coefficient for p-adaptive
     836                 :            :     //!   solution polynomials. Unlike compflow class, the high order of fv
     837                 :            :     //!   solution will not be reset since p0p1 is the base scheme for
     838                 :            :     //!   multi-material p-adaptive DG method.
     839                 :          0 :     void resetAdapSol( const inciter::FaceData& fd,
     840                 :            :                        tk::Fields& unk,
     841                 :            :                        tk::Fields& prim,
     842                 :            :                        const std::vector< std::size_t >& ndofel ) const
     843                 :            :     {
     844                 :          0 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     845                 :          0 :       const auto ncomp = unk.nprop() / rdof;
     846                 :          0 :       const auto nprim = prim.nprop() / rdof;
     847                 :            : 
     848         [ -  - ]:          0 :       for(std::size_t e = 0; e < fd.Esuel().size()/4; e++)
     849                 :            :       {
     850         [ -  - ]:          0 :         if(ndofel[e] < 10)
     851                 :            :         {
     852         [ -  - ]:          0 :           for (std::size_t c=0; c<ncomp; ++c)
     853                 :            :           {
     854                 :          0 :             auto mark = c*rdof;
     855                 :          0 :             unk(e, mark+4) = 0.0;
     856                 :          0 :             unk(e, mark+5) = 0.0;
     857                 :          0 :             unk(e, mark+6) = 0.0;
     858                 :          0 :             unk(e, mark+7) = 0.0;
     859                 :          0 :             unk(e, mark+8) = 0.0;
     860                 :          0 :             unk(e, mark+9) = 0.0;
     861                 :            :           }
     862         [ -  - ]:          0 :           for (std::size_t c=0; c<nprim; ++c)
     863                 :            :           {
     864                 :          0 :             auto mark = c*rdof;
     865                 :          0 :             prim(e, mark+4) = 0.0;
     866                 :          0 :             prim(e, mark+5) = 0.0;
     867                 :          0 :             prim(e, mark+6) = 0.0;
     868                 :          0 :             prim(e, mark+7) = 0.0;
     869                 :          0 :             prim(e, mark+8) = 0.0;
     870                 :          0 :             prim(e, mark+9) = 0.0;
     871                 :            :           }
     872                 :            :         }
     873                 :            :       }
     874                 :          0 :     }
     875                 :            : 
     876                 :            :     //! Compute right hand side
     877                 :            :     //! \param[in] t Physical time
     878                 :            :     //! \param[in] pref Indicator for p-adaptive algorithm
     879                 :            :     //! \param[in] geoFace Face geometry array
     880                 :            :     //! \param[in] geoElem Element geometry array
     881                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     882                 :            :     //! \param[in] inpoel Element-node connectivity
     883                 :            :     //! \param[in] coord Array of nodal coordinates
     884                 :            :     //! \param[in] U Solution vector at recent time step
     885                 :            :     //! \param[in] P Primitive vector at recent time step
     886                 :            :     //! \param[in] ndofel Vector of local number of degrees of freedom
     887                 :            :     //! \param[in] dt Delta time
     888                 :            :     //! \param[in,out] R Right-hand side vector computed
     889                 :       5205 :     void rhs( tk::real t,
     890                 :            :               const bool pref,
     891                 :            :               const tk::Fields& geoFace,
     892                 :            :               const tk::Fields& geoElem,
     893                 :            :               const inciter::FaceData& fd,
     894                 :            :               const std::vector< std::size_t >& inpoel,
     895                 :            :               const std::vector< std::unordered_set< std::size_t > >&,
     896                 :            :               const tk::UnsMesh::Coords& coord,
     897                 :            :               const tk::Fields& U,
     898                 :            :               const tk::Fields& P,
     899                 :            :               const std::vector< std::size_t >& ndofel,
     900                 :            :               const tk::real dt,
     901                 :            :               tk::Fields& R ) const
     902                 :            :     {
     903                 :       5205 :       const auto ndof = g_inputdeck.get< tag::ndof >();
     904                 :       5205 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     905                 :       5205 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
     906                 :       5205 :       const auto intsharp =
     907                 :       5205 :         g_inputdeck.get< tag::multimat, tag::intsharp >();
     908                 :            :       const auto& solidx = inciter::g_inputdeck.get<
     909                 :       5205 :         tag::matidxmap, tag::solidx >();
     910         [ +  - ]:       5205 :       auto nsld = numSolids(nmat, solidx);
     911                 :            : 
     912                 :       5205 :       const auto nelem = fd.Esuel().size()/4;
     913                 :            : 
     914 [ -  + ][ -  - ]:       5205 :       Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     915                 :            :               "vector and primitive vector at recent time step incorrect" );
     916 [ -  + ][ -  - ]:       5205 :       Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
     917                 :            :               "vector and right-hand side at recent time step incorrect" );
     918 [ -  + ][ -  - ]:       5205 :       Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
     919                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     920 [ -  + ][ -  - ]:       5205 :       Assert( P.nprop() == rdof*m_nprim, "Number of components in primitive "
         [ -  - ][ -  - ]
                 [ -  - ]
     921                 :            :               "vector must equal "+ std::to_string(rdof*m_nprim) );
     922 [ -  + ][ -  - ]:       5205 :       Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
         [ -  - ][ -  - ]
                 [ -  - ]
     923                 :            :               "side vector must equal "+ std::to_string(ndof*m_ncomp) );
     924 [ -  + ][ -  - ]:       5205 :       Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
         [ -  - ][ -  - ]
     925                 :            :               "Mismatch in inpofa size" );
     926                 :            : 
     927                 :            :       // set rhs to zero
     928         [ +  - ]:       5205 :       R.fill(0.0);
     929                 :            : 
     930                 :            :       // Allocate space for Riemann derivatives used in non-conservative terms.
     931                 :            :       // The following Riemann derivatives are stored, in order:
     932                 :            :       // 1) 3*nmat terms: derivatives of partial pressure of each material,
     933                 :            :       //    for the energy equations.
     934                 :            :       // 2) ndof terms: derivatives of Riemann velocity times the basis
     935                 :            :       //    function, for the volume fraction equations.
     936                 :            :       // 3) nmat*3*3*9 terms: 3 derivatives of u_l*g_ij for each material, for
     937                 :            :       //    the deformation gradient equations.
     938                 :            :       // 4) 3*nsld terms: 3 derivatives of \alpha \sigma_ij for each solid
     939                 :            :       //    material, for the energy equations.
     940                 :            :       std::vector< std::vector< tk::real > >
     941 [ +  - ][ +  - ]:      15615 :         riemannDeriv(3*nmat+ndof+3*nsld, std::vector<tk::real>(U.nunk(),0.0));
     942                 :            : 
     943                 :            :       // configure a no-op lambda for prescribed velocity
     944                 :    9919425 :       auto velfn = []( ncomp_t, tk::real, tk::real, tk::real, tk::real ){
     945                 :    9919425 :         return tk::VelFn::result_type(); };
     946                 :            : 
     947                 :            :       // compute internal surface flux integrals
     948         [ +  - ]:      10410 :       tk::surfInt( pref, nmat, m_mat_blk, t, ndof, rdof, inpoel, solidx,
     949         [ +  - ]:       5205 :                    coord, fd, geoFace, geoElem, m_riemann, velfn, U, P, ndofel,
     950                 :            :                    dt, R, riemannDeriv, intsharp );
     951                 :            : 
     952                 :            :       // compute optional source term
     953 [ +  - ][ +  - ]:       5205 :       tk::srcInt( m_mat_blk, t, ndof, fd.Esuel().size()/4, inpoel,
     954                 :            :                   coord, geoElem, Problem::src, ndofel, R, nmat );
     955                 :            : 
     956         [ +  + ]:       5205 :       if(ndof > 1)
     957                 :            :         // compute volume integrals
     958 [ +  - ][ +  - ]:        780 :         tk::volInt( nmat, t, m_mat_blk, ndof, rdof, nelem,
                 [ +  - ]
     959                 :            :                     inpoel, coord, geoElem, flux, velfn, U, P, ndofel, R,
     960                 :            :                     intsharp );
     961                 :            : 
     962                 :            :       // compute boundary surface flux integrals
     963         [ +  + ]:      36435 :       for (const auto& b : m_bc)
     964 [ +  - ][ +  - ]:      62460 :         tk::bndSurfInt( pref, nmat, m_mat_blk, ndof, rdof,
     965                 :      31230 :                         std::get<0>(b), fd, geoFace, geoElem, inpoel, coord, t,
     966                 :      62460 :                         m_riemann, velfn, std::get<1>(b), U, P, ndofel, R,
     967                 :            :                         riemannDeriv, intsharp );
     968                 :            : 
     969 [ -  + ][ -  - ]:       5205 :       Assert( riemannDeriv.size() == 3*nmat+ndof+3*nsld, "Size of "
         [ -  - ][ -  - ]
     970                 :            :               "Riemann derivative vector incorrect" );
     971                 :            : 
     972                 :            :       // get derivatives from riemannDeriv
     973         [ +  + ]:      46230 :       for (std::size_t k=0; k<riemannDeriv.size(); ++k)
     974                 :            :       {
     975 [ -  + ][ -  - ]:      41025 :         Assert( riemannDeriv[k].size() == U.nunk(), "Riemann derivative vector "
         [ -  - ][ -  - ]
     976                 :            :                 "for non-conservative terms has incorrect size" );
     977         [ +  + ]:   28165725 :         for (std::size_t e=0; e<U.nunk(); ++e)
     978         [ +  - ]:   28124700 :           riemannDeriv[k][e] /= geoElem(e, 0);
     979                 :            :       }
     980                 :            : 
     981                 :            :       // compute volume integrals of non-conservative terms
     982         [ +  - ]:       5205 :       tk::nonConservativeInt( pref, nmat, m_mat_blk, ndof, rdof, nelem,
     983                 :            :                               inpoel, coord, geoElem, U, P, riemannDeriv,
     984                 :            :                               ndofel, R, intsharp );
     985                 :            : 
     986                 :            :       // Compute integrals for inverse deformation correction in solid materials
     987 [ +  - ][ -  + ]:       5205 :       if (inciter::haveSolid(nmat, solidx) &&
         [ -  - ][ -  + ]
     988                 :          0 :         g_inputdeck.get< tag::multimat, tag::rho0constraint >())
     989         [ -  - ]:          0 :         tk::solidTermsVolInt( nmat, m_mat_blk, ndof, rdof, nelem,
     990                 :            :                               inpoel, coord, geoElem, U, P, ndofel,
     991                 :            :                               dt, R);
     992                 :            : 
     993                 :            :       // compute finite pressure relaxation terms
     994         [ +  + ]:       5205 :       if (g_inputdeck.get< tag::multimat, tag::prelax >())
     995                 :            :       {
     996                 :        375 :         const auto ct = g_inputdeck.get< tag::multimat,
     997                 :        375 :                                          tag::prelax_timescale >();
     998         [ +  - ]:        375 :         tk::pressureRelaxationInt( pref, nmat, m_mat_blk, ndof,
     999                 :            :                                    rdof, nelem, inpoel, coord, geoElem, U, P,
    1000                 :            :                                    ndofel, ct, R, intsharp );
    1001                 :            :       }
    1002                 :       5205 :     }
    1003                 :            : 
    1004                 :            :     //! Evaluate the adaptive indicator and mark the ndof for each element
    1005                 :            :     //! \param[in] nunk Number of unknowns
    1006                 :            :     //! \param[in] coord Array of nodal coordinates
    1007                 :            :     //! \param[in] inpoel Element-node connectivity
    1008                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
    1009                 :            :     //! \param[in] unk Array of unknowns
    1010                 :            :     //! \param[in] prim Array of primitive quantities
    1011                 :            :     //! \param[in] indicator p-refinement indicator type
    1012                 :            :     //! \param[in] ndof Number of degrees of freedom in the solution
    1013                 :            :     //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
    1014                 :            :     //! \param[in] tolref Tolerance for p-refinement
    1015                 :            :     //! \param[in,out] ndofel Vector of local number of degrees of freedome
    1016                 :          0 :     void eval_ndof( std::size_t nunk,
    1017                 :            :                     [[maybe_unused]] const tk::UnsMesh::Coords& coord,
    1018                 :            :                     [[maybe_unused]] const std::vector< std::size_t >& inpoel,
    1019                 :            :                     const inciter::FaceData& fd,
    1020                 :            :                     const tk::Fields& unk,
    1021                 :            :                     [[maybe_unused]] const tk::Fields& prim,
    1022                 :            :                     inciter::ctr::PrefIndicatorType indicator,
    1023                 :            :                     std::size_t ndof,
    1024                 :            :                     std::size_t ndofmax,
    1025                 :            :                     tk::real tolref,
    1026                 :            :                     std::vector< std::size_t >& ndofel ) const
    1027                 :            :     {
    1028                 :          0 :       const auto& esuel = fd.Esuel();
    1029                 :          0 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
    1030                 :            : 
    1031         [ -  - ]:          0 :       if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
    1032                 :          0 :         spectral_decay(nmat, nunk, esuel, unk, ndof, ndofmax, tolref, ndofel);
    1033                 :            :       else
    1034 [ -  - ][ -  - ]:          0 :         Throw( "No such adaptive indicator type" );
                 [ -  - ]
    1035                 :          0 :     }
    1036                 :            : 
    1037                 :            :     //! Compute the minimum time step size
    1038                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
    1039                 :            :     //! \param[in] geoFace Face geometry array
    1040                 :            :     //! \param[in] geoElem Element geometry array
    1041                 :            : //    //! \param[in] ndofel Vector of local number of degrees of freedom
    1042                 :            :     //! \param[in] U Solution vector at recent time step
    1043                 :            :     //! \param[in] P Vector of primitive quantities at recent time step
    1044                 :            :     //! \param[in] nielem Number of internal elements
    1045                 :            :     //! \return Minimum time step size
    1046                 :            :     //! \details The allowable dt is calculated by looking at the maximum
    1047                 :            :     //!   wave-speed in elements surrounding each face, times the area of that
    1048                 :            :     //!   face. Once the maximum of this quantity over the mesh is determined,
    1049                 :            :     //!   the volume of each cell is divided by this quantity. A minimum of this
    1050                 :            :     //!   ratio is found over the entire mesh, which gives the allowable dt.
    1051                 :        260 :     tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
    1052                 :            :                  const std::vector< std::size_t >&,
    1053                 :            :                  const inciter::FaceData& fd,
    1054                 :            :                  const tk::Fields& geoFace,
    1055                 :            :                  const tk::Fields& geoElem,
    1056                 :            :                  const std::vector< std::size_t >& /*ndofel*/,
    1057                 :            :                  const tk::Fields& U,
    1058                 :            :                  const tk::Fields& P,
    1059                 :            :                  const std::size_t nielem ) const
    1060                 :            :     {
    1061                 :        260 :       const auto ndof = g_inputdeck.get< tag::ndof >();
    1062                 :        260 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
    1063                 :            : 
    1064                 :        260 :       auto mindt = timeStepSizeMultiMat( m_mat_blk, fd.Esuf(), geoFace, geoElem,
    1065                 :            :         nielem, nmat, U, P);
    1066                 :            : 
    1067                 :        260 :       tk::real dgp = 0.0;
    1068         [ +  - ]:        260 :       if (ndof == 4)
    1069                 :            :       {
    1070                 :        260 :         dgp = 1.0;
    1071                 :            :       }
    1072         [ -  - ]:          0 :       else if (ndof == 10)
    1073                 :            :       {
    1074                 :          0 :         dgp = 2.0;
    1075                 :            :       }
    1076                 :            : 
    1077                 :            :       // Scale smallest dt with CFL coefficient and the CFL is scaled by (2*p+1)
    1078                 :            :       // where p is the order of the DG polynomial by linear stability theory.
    1079                 :        260 :       mindt /= (2.0*dgp + 1.0);
    1080                 :        260 :       return mindt;
    1081                 :            :     }
    1082                 :            : 
    1083                 :            :     //! Compute stiff terms for a single element
    1084                 :            :     //! \param[in] e Element number
    1085                 :            :     //! \param[in] geoElem Element geometry array
    1086                 :            :     //! \param[in] inpoel Element-node connectivity
    1087                 :            :     //! \param[in] coord Array of nodal coordinates
    1088                 :            :     //! \param[in] U Solution vector at recent time step
    1089                 :            :     //! \param[in] P Primitive vector at recent time step
    1090                 :            :     //! \param[in] ndofel Vector of local number of degrees of freedom
    1091                 :            :     //! \param[in,out] R Right-hand side vector computed
    1092                 :          0 :     void stiff_rhs( std::size_t e,
    1093                 :            :                     const tk::Fields& geoElem,
    1094                 :            :                     const std::vector< std::size_t >& inpoel,
    1095                 :            :                     const tk::UnsMesh::Coords& coord,
    1096                 :            :                     const tk::Fields& U,
    1097                 :            :                     const tk::Fields& P,
    1098                 :            :                     const std::vector< std::size_t >& ndofel,
    1099                 :            :                     tk::Fields& R ) const
    1100                 :            :     {
    1101                 :          0 :       const auto ndof = g_inputdeck.get< tag::ndof >();
    1102                 :          0 :       const auto rdof = g_inputdeck.get< tag::rdof >();
    1103                 :          0 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
    1104                 :          0 :       const auto intsharp =
    1105                 :          0 :         g_inputdeck.get< tag::multimat, tag::intsharp >();
    1106                 :            :       const auto& solidx = inciter::g_inputdeck.get<
    1107                 :          0 :         tag::matidxmap, tag::solidx >();
    1108                 :            : 
    1109 [ -  - ][ -  - ]:          0 :       Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
         [ -  - ][ -  - ]
    1110                 :            :               "vector and primitive vector at recent time step incorrect" );
    1111 [ -  - ][ -  - ]:          0 :       Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
         [ -  - ][ -  - ]
                 [ -  - ]
    1112                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
    1113 [ -  - ][ -  - ]:          0 :       Assert( P.nprop() == rdof*m_nprim, "Number of components in primitive "
         [ -  - ][ -  - ]
                 [ -  - ]
    1114                 :            :               "vector must equal "+ std::to_string(rdof*m_nprim) );
    1115 [ -  - ][ -  - ]:          0 :       Assert( R.nprop() == ndof*nstiffeq(), "Number of components in "
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
                 [ -  - ]
    1116                 :            :               "right-hand side must equal "+ std::to_string(ndof*nstiffeq()) );
    1117                 :            : 
    1118                 :            :       // set rhs to zero for element e
    1119 [ -  - ][ -  - ]:          0 :       for (std::size_t i=0; i<ndof*nstiffeq(); ++i)
    1120         [ -  - ]:          0 :         R(e, i) = 0.0;
    1121                 :            : 
    1122                 :          0 :       const auto& cx = coord[0];
    1123                 :          0 :       const auto& cy = coord[1];
    1124                 :          0 :       const auto& cz = coord[2];
    1125                 :            : 
    1126                 :          0 :       auto ncomp = U.nprop()/rdof;
    1127                 :          0 :       auto nprim = P.nprop()/rdof;
    1128                 :            : 
    1129         [ -  - ]:          0 :       auto ng = tk::NGvol(ndofel[e]);
    1130                 :            : 
    1131                 :            :       // arrays for quadrature points
    1132                 :          0 :       std::array< std::vector< tk::real >, 3 > coordgp;
    1133                 :          0 :       std::vector< tk::real > wgp;
    1134                 :            : 
    1135         [ -  - ]:          0 :       coordgp[0].resize( ng );
    1136         [ -  - ]:          0 :       coordgp[1].resize( ng );
    1137         [ -  - ]:          0 :       coordgp[2].resize( ng );
    1138         [ -  - ]:          0 :       wgp.resize( ng );
    1139                 :            : 
    1140         [ -  - ]:          0 :       tk::GaussQuadratureTet( ng, coordgp, wgp );
    1141                 :            : 
    1142                 :            :       // Extract the element coordinates
    1143                 :          0 :       std::array< std::array< tk::real, 3>, 4 > coordel {{
    1144                 :            :         {{ cx[ inpoel[4*e  ] ], cy[ inpoel[4*e  ] ], cz[ inpoel[4*e  ] ] }},
    1145                 :          0 :         {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
    1146                 :          0 :         {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
    1147                 :          0 :         {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
    1148                 :            :       }};
    1149                 :            : 
    1150                 :            :       // Gaussian quadrature
    1151         [ -  - ]:          0 :       for (std::size_t igp=0; igp<ng; ++igp)
    1152                 :            :       {
    1153                 :            :         // Compute the coordinates of quadrature point at physical domain
    1154         [ -  - ]:          0 :         auto gp = tk::eval_gp( igp, coordel, coordgp );
    1155                 :            : 
    1156                 :            :         // Compute the basis function
    1157         [ -  - ]:          0 :         auto B = tk::eval_basis( ndofel[e], coordgp[0][igp], coordgp[1][igp],
    1158                 :            :                              coordgp[2][igp] );
    1159                 :            : 
    1160         [ -  - ]:          0 :         auto state = tk::evalPolynomialSol(m_mat_blk, intsharp, ncomp, nprim,
    1161                 :            :           rdof, nmat, e, ndofel[e], inpoel, coord, geoElem, gp, B, U, P);
    1162                 :            : 
    1163                 :            :         // compute source
    1164                 :            :         // Loop through materials
    1165                 :          0 :         std::size_t ksld = 0;
    1166         [ -  - ]:          0 :         for (std::size_t k=0; k<nmat; ++k)
    1167                 :            :         {
    1168         [ -  - ]:          0 :           if (solidx[k] > 0)
    1169                 :            :           {
    1170                 :          0 :             tk::real alpha = state[inciter::volfracIdx(nmat, k)];
    1171                 :            :             std::array< std::array< tk::real, 3 >, 3 > g;
    1172                 :            :             // Compute the source terms
    1173         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1174         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
    1175                 :          0 :                 g[i][j] = state[inciter::deformIdx(nmat,solidx[k],i,j)];
    1176                 :            : 
    1177                 :            :             // Compute Lp
    1178                 :            :             // Reference: Ortega, A. L., Lombardini, M., Pullin, D. I., &
    1179                 :            :             // Meiron, D. I. (2014). Numerical simulation of elastic–plastic
    1180                 :            :             // solid mechanics using an Eulerian stretch tensor approach and
    1181                 :            :             // HLLD Riemann solver. Journal of Computational Physics, 257,
    1182                 :            :             // 414-441
    1183                 :            :             std::array< std::array< tk::real, 3 >, 3 > Lp;
    1184                 :            : 
    1185                 :            :             // 1. Compute dev(sigma)
    1186                 :          0 :             auto sigma_dev = m_mat_blk[k].computeTensor< EOS::CauchyStress >(
    1187         [ -  - ]:          0 :               0.0, 0.0, 0.0, 0.0, 0.0, alpha, k, g );
    1188                 :          0 :             tk::real apr = state[ncomp+inciter::pressureIdx(nmat, k)];
    1189         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i) sigma_dev[i][i] -= apr;
    1190         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1191         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
    1192                 :          0 :                 sigma_dev[i][j] /= alpha;
    1193                 :          0 :             tk::real sigma_trace =
    1194                 :          0 :               sigma_dev[0][0]+sigma_dev[1][1]+sigma_dev[2][2];
    1195         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1196                 :          0 :               sigma_dev[i][i] -= sigma_trace/3.0;
    1197                 :            : 
    1198                 :            :             // 2. Compute inv(g)
    1199                 :            :             double ginv[9];
    1200         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1201         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
    1202                 :          0 :                 ginv[3*i+j] = g[i][j];
    1203                 :            :             lapack_int ipiv[3];
    1204                 :            :             #ifndef NDEBUG
    1205                 :            :             lapack_int ierr =
    1206                 :            :             #endif
    1207         [ -  - ]:          0 :               LAPACKE_dgetrf(LAPACK_ROW_MAJOR, 3, 3, ginv, 3, ipiv);
    1208 [ -  - ][ -  - ]:          0 :             Assert(ierr==0, "Lapack error in LU factorization of g");
         [ -  - ][ -  - ]
    1209                 :            :             #ifndef NDEBUG
    1210                 :            :             lapack_int jerr =
    1211                 :            :             #endif
    1212         [ -  - ]:          0 :               LAPACKE_dgetri(LAPACK_ROW_MAJOR, 3, ginv, 3, ipiv);
    1213 [ -  - ][ -  - ]:          0 :             Assert(jerr==0, "Lapack error in inverting g");
         [ -  - ][ -  - ]
    1214                 :            : 
    1215                 :            :             // 3. Compute dev(sigma)*inv(g)
    1216                 :            :             std::array< std::array< tk::real, 3 >, 3 > aux_mat;
    1217         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1218         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
    1219                 :            :               {
    1220                 :          0 :                 tk::real sum = 0.0;
    1221         [ -  - ]:          0 :                 for (std::size_t l=0; l<3; ++l)
    1222                 :          0 :                   sum += sigma_dev[i][l]*ginv[3*l+j];
    1223                 :          0 :                 aux_mat[i][j] = sum;
    1224                 :            :               }
    1225                 :            : 
    1226                 :            :             // 4. Compute g*(dev(sigma)*inv(g))
    1227         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1228         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
    1229                 :            :               {
    1230                 :          0 :                 tk::real sum = 0.0;
    1231         [ -  - ]:          0 :                 for (std::size_t l=0; l<3; ++l)
    1232                 :          0 :                   sum += g[i][l]*aux_mat[l][j];
    1233                 :          0 :                 Lp[i][j] = sum;
    1234                 :            :               }
    1235                 :            : 
    1236                 :            :             // 5. Divide by 2*mu*tau
    1237                 :            :             // 'Perfect' plasticity
    1238         [ -  - ]:          0 :             tk::real yield_stress = getmatprop< tag::yield_stress >(k);
    1239                 :          0 :             tk::real equiv_stress = 0.0;
    1240         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1241         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
    1242                 :          0 :                 equiv_stress += sigma_dev[i][j]*sigma_dev[i][j];
    1243                 :          0 :             equiv_stress = std::sqrt(3.0*equiv_stress/2.0);
    1244                 :            :             // rel_factor = 1/tau <- Perfect plasticity for now.
    1245                 :          0 :             tk::real rel_factor = 0.0;
    1246         [ -  - ]:          0 :             if (equiv_stress >= yield_stress)
    1247                 :          0 :               rel_factor = 1.0e07;
    1248         [ -  - ]:          0 :             tk::real mu = getmatprop< tag::mu >(k);
    1249         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1250         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
    1251                 :          0 :                 Lp[i][j] *= rel_factor/(2.0*mu);
    1252                 :            : 
    1253                 :            :             // Compute the source terms
    1254         [ -  - ]:          0 :             std::vector< tk::real > s(9*ndof, 0.0);
    1255         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1256         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
    1257         [ -  - ]:          0 :                 for (std::size_t idof=0; idof<ndof; ++idof)
    1258                 :            :                 {
    1259                 :          0 :                   s[(i*3+j)*ndof+idof] = B[idof] * (Lp[i][0]*g[0][j]
    1260                 :          0 :                                                    +Lp[i][1]*g[1][j]
    1261                 :          0 :                                                    +Lp[i][2]*g[2][j]);
    1262                 :            :                 }
    1263                 :            : 
    1264         [ -  - ]:          0 :             auto wt = wgp[igp] * geoElem(e, 0);
    1265                 :            : 
    1266                 :            :             // Contribute to the right-hand-side
    1267         [ -  - ]:          0 :             for (std::size_t i=0; i<3; ++i)
    1268         [ -  - ]:          0 :               for (std::size_t j=0; j<3; ++j)
    1269         [ -  - ]:          0 :                 for (std::size_t idof=0; idof<ndof; ++idof)
    1270                 :            :                 {
    1271                 :          0 :                   std::size_t srcId = (i*3+j)*ndof+idof;
    1272                 :          0 :                   std::size_t dofId = solidTensorIdx(ksld,i,j)*ndof+idof;
    1273         [ -  - ]:          0 :                   R(e, dofId) += wt * s[srcId];
    1274                 :            :                 }
    1275                 :            : 
    1276                 :          0 :             ksld++;
    1277                 :            :           }
    1278                 :            :         }
    1279                 :            : 
    1280                 :            :         }
    1281                 :          0 :     }
    1282                 :            : 
    1283                 :            :     //! Extract the velocity field at cell nodes. Currently unused.
    1284                 :            :     //! \param[in] U Solution vector at recent time step
    1285                 :            :     //! \param[in] N Element node indices
    1286                 :            :     //! \return Array of the four values of the velocity field
    1287                 :            :     std::array< std::array< tk::real, 4 >, 3 >
    1288                 :            :     velocity( const tk::Fields& U,
    1289                 :            :               const std::array< std::vector< tk::real >, 3 >&,
    1290                 :            :               const std::array< std::size_t, 4 >& N ) const
    1291                 :            :     {
    1292                 :            :       const auto rdof = g_inputdeck.get< tag::rdof >();
    1293                 :            :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
    1294                 :            : 
    1295                 :            :       std::array< std::array< tk::real, 4 >, 3 > v;
    1296                 :            :       v[0] = U.extract( momentumDofIdx(nmat, 0, rdof, 0), N );
    1297                 :            :       v[1] = U.extract( momentumDofIdx(nmat, 1, rdof, 0), N );
    1298                 :            :       v[2] = U.extract( momentumDofIdx(nmat, 2, rdof, 0), N );
    1299                 :            : 
    1300                 :            :       std::vector< std::array< tk::real, 4 > > ar;
    1301                 :            :       ar.resize(nmat);
    1302                 :            :       for (std::size_t k=0; k<nmat; ++k)
    1303                 :            :         ar[k] = U.extract( densityDofIdx(nmat, k, rdof, 0), N );
    1304                 :            : 
    1305                 :            :       std::array< tk::real, 4 > r{{ 0.0, 0.0, 0.0, 0.0 }};
    1306                 :            :       for (std::size_t i=0; i<r.size(); ++i) {
    1307                 :            :         for (std::size_t k=0; k<nmat; ++k)
    1308                 :            :           r[i] += ar[k][i];
    1309                 :            :       }
    1310                 :            : 
    1311                 :            :       std::transform( r.begin(), r.end(), v[0].begin(), v[0].begin(),
    1312                 :            :                       []( tk::real s, tk::real& d ){ return d /= s; } );
    1313                 :            :       std::transform( r.begin(), r.end(), v[1].begin(), v[1].begin(),
    1314                 :            :                       []( tk::real s, tk::real& d ){ return d /= s; } );
    1315                 :            :       std::transform( r.begin(), r.end(), v[2].begin(), v[2].begin(),
    1316                 :            :                       []( tk::real s, tk::real& d ){ return d /= s; } );
    1317                 :            :       return v;
    1318                 :            :     }
    1319                 :            : 
    1320                 :            :     //! Return a map that associates user-specified strings to functions
    1321                 :            :     //! \return Map that associates user-specified strings to functions that
    1322                 :            :     //!   compute relevant quantities to be output to file
    1323                 :        310 :     std::map< std::string, tk::GetVarFn > OutVarFn() const
    1324                 :        310 :     { return MultiMatOutVarFn(); }
    1325                 :            : 
    1326                 :            :     //! Return analytic field names to be output to file
    1327                 :            :     //! \return Vector of strings labelling analytic fields output in file
    1328                 :          0 :     std::vector< std::string > analyticFieldNames() const {
    1329                 :          0 :       auto nmat = g_inputdeck.get< eq, tag::nmat >();
    1330                 :            : 
    1331                 :          0 :       return MultiMatFieldNames(nmat);
    1332                 :            :     }
    1333                 :            : 
    1334                 :            :     //! Return time history field names to be output to file
    1335                 :            :     //! \return Vector of strings labelling time history fields output in file
    1336                 :          0 :     std::vector< std::string > histNames() const {
    1337                 :          0 :       return MultiMatHistNames();
    1338                 :            :     }
    1339                 :            : 
    1340                 :            :     //! Return surface field output going to file
    1341                 :            :     std::vector< std::vector< tk::real > >
    1342                 :          0 :     surfOutput( const std::map< int, std::vector< std::size_t > >&,
    1343                 :            :                 tk::Fields& ) const
    1344                 :            :     {
    1345                 :          0 :       std::vector< std::vector< tk::real > > s; // punt for now
    1346                 :          0 :       return s;
    1347                 :            :     }
    1348                 :            : 
    1349                 :            :     //! Return time history field output evaluated at time history points
    1350                 :            :     //! \param[in] h History point data
    1351                 :            :     //! \param[in] inpoel Element-node connectivity
    1352                 :            :     //! \param[in] coord Array of nodal coordinates
    1353                 :            :     //! \param[in] U Array of unknowns
    1354                 :            :     //! \param[in] P Array of primitive quantities
    1355                 :            :     //! \return Vector of time history output of bulk flow quantities (density,
    1356                 :            :     //!   velocity, total energy, and pressure) evaluated at time history points
    1357                 :            :     std::vector< std::vector< tk::real > >
    1358                 :          0 :     histOutput( const std::vector< HistData >& h,
    1359                 :            :                 const std::vector< std::size_t >& inpoel,
    1360                 :            :                 const tk::UnsMesh::Coords& coord,
    1361                 :            :                 const tk::Fields& U,
    1362                 :            :                 const tk::Fields& P ) const
    1363                 :            :     {
    1364                 :          0 :       const auto rdof = g_inputdeck.get< tag::rdof >();
    1365                 :          0 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
    1366                 :            : 
    1367                 :          0 :       const auto& x = coord[0];
    1368                 :          0 :       const auto& y = coord[1];
    1369                 :          0 :       const auto& z = coord[2];
    1370                 :            : 
    1371         [ -  - ]:          0 :       std::vector< std::vector< tk::real > > Up(h.size());
    1372                 :            : 
    1373                 :          0 :       std::size_t j = 0;
    1374         [ -  - ]:          0 :       for (const auto& p : h) {
    1375                 :          0 :         auto e = p.get< tag::elem >();
    1376                 :          0 :         auto chp = p.get< tag::coord >();
    1377                 :            : 
    1378                 :            :         // Evaluate inverse Jacobian
    1379                 :          0 :         std::array< std::array< tk::real, 3>, 4 > cp{{
    1380                 :            :           {{ x[inpoel[4*e  ]], y[inpoel[4*e  ]], z[inpoel[4*e  ]] }},
    1381                 :          0 :           {{ x[inpoel[4*e+1]], y[inpoel[4*e+1]], z[inpoel[4*e+1]] }},
    1382                 :          0 :           {{ x[inpoel[4*e+2]], y[inpoel[4*e+2]], z[inpoel[4*e+2]] }},
    1383                 :          0 :           {{ x[inpoel[4*e+3]], y[inpoel[4*e+3]], z[inpoel[4*e+3]] }} }};
    1384                 :          0 :         auto J = tk::inverseJacobian( cp[0], cp[1], cp[2], cp[3] );
    1385                 :            : 
    1386                 :            :         // evaluate solution at history-point
    1387                 :          0 :         std::array< tk::real, 3 > dc{{chp[0]-cp[0][0], chp[1]-cp[0][1],
    1388                 :          0 :           chp[2]-cp[0][2]}};
    1389         [ -  - ]:          0 :         auto B = tk::eval_basis(rdof, tk::dot(J[0],dc), tk::dot(J[1],dc),
    1390                 :          0 :           tk::dot(J[2],dc));
    1391         [ -  - ]:          0 :         auto uhp = eval_state(m_ncomp, rdof, rdof, e, U, B);
    1392         [ -  - ]:          0 :         auto php = eval_state(m_nprim, rdof, rdof, e, P, B);
    1393                 :            : 
    1394                 :            :         // store solution in history output vector
    1395         [ -  - ]:          0 :         Up[j].resize(6, 0.0);
    1396         [ -  - ]:          0 :         for (std::size_t k=0; k<nmat; ++k) {
    1397                 :          0 :           Up[j][0] += uhp[densityIdx(nmat,k)];
    1398                 :          0 :           Up[j][4] += uhp[energyIdx(nmat,k)];
    1399                 :          0 :           Up[j][5] += php[pressureIdx(nmat,k)];
    1400                 :            :         }
    1401                 :          0 :         Up[j][1] = php[velocityIdx(nmat,0)];
    1402                 :          0 :         Up[j][2] = php[velocityIdx(nmat,1)];
    1403                 :          0 :         Up[j][3] = php[velocityIdx(nmat,2)];
    1404                 :          0 :         ++j;
    1405                 :            :       }
    1406                 :            : 
    1407                 :          0 :       return Up;
    1408                 :            :     }
    1409                 :            : 
    1410                 :            :     //! Return names of integral variables to be output to diagnostics file
    1411                 :            :     //! \return Vector of strings labelling integral variables output
    1412                 :         12 :     std::vector< std::string > names() const
    1413                 :            :     {
    1414                 :         12 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
    1415                 :         12 :       return MultiMatDiagNames(nmat);
    1416                 :            :     }
    1417                 :            : 
    1418                 :            :     //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
    1419                 :            :     //! \param[in] xi X-coordinate at which to evaluate the analytic solution
    1420                 :            :     //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
    1421                 :            :     //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
    1422                 :            :     //! \param[in] t Physical time at which to evaluate the analytic solution
    1423                 :            :     //! \return Vector of analytic solution at given location and time
    1424                 :            :     std::vector< tk::real >
    1425                 :          0 :     analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
    1426                 :          0 :     { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
    1427                 :            : 
    1428                 :            :     //! Return analytic solution for conserved variables
    1429                 :            :     //! \param[in] xi X-coordinate at which to evaluate the analytic solution
    1430                 :            :     //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
    1431                 :            :     //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
    1432                 :            :     //! \param[in] t Physical time at which to evaluate the analytic solution
    1433                 :            :     //! \return Vector of analytic solution at given location and time
    1434                 :            :     std::vector< tk::real >
    1435                 :    1322920 :     solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
    1436                 :    1322920 :     { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
    1437                 :            : 
    1438                 :            :     //! Return cell-averaged specific total energy for an element
    1439                 :            :     //! \param[in] e Element id for which total energy is required
    1440                 :            :     //! \param[in] unk Vector of conserved quantities
    1441                 :            :     //! \return Cell-averaged specific total energy for given element
    1442                 :     822640 :     tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
    1443                 :            :     {
    1444                 :     822640 :       const auto rdof = g_inputdeck.get< tag::rdof >();
    1445                 :     822640 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
    1446                 :            : 
    1447                 :     822640 :       tk::real sp_te(0.0);
    1448                 :            :       // sum each material total energy
    1449         [ +  + ]:    2858500 :       for (std::size_t k=0; k<nmat; ++k) {
    1450                 :    2035860 :         sp_te += unk(e, energyDofIdx(nmat,k,rdof,0));
    1451                 :            :       }
    1452                 :     822640 :       return sp_te;
    1453                 :            :     }
    1454                 :            : 
    1455                 :            :   private:
    1456                 :            :     //! Number of components in this PDE system
    1457                 :            :     const ncomp_t m_ncomp;
    1458                 :            :     //! Number of primitive quantities stored in this PDE system
    1459                 :            :     const ncomp_t m_nprim;
    1460                 :            :     //! Riemann solver
    1461                 :            :     tk::RiemannFluxFn m_riemann;
    1462                 :            :     //! BC configuration
    1463                 :            :     BCStateFn m_bc;
    1464                 :            :     //! EOS material block
    1465                 :            :     std::vector< EOS > m_mat_blk;
    1466                 :            : 
    1467                 :            :     //! Evaluate conservative part of physical flux function for this PDE system
    1468                 :            :     //! \param[in] ncomp Number of scalar components in this PDE system
    1469                 :            :     //! \param[in] ugp Numerical solution at the Gauss point at which to
    1470                 :            :     //!   evaluate the flux
    1471                 :            :     //! \return Flux vectors for all components in this PDE system
    1472                 :            :     //! \note The function signature must follow tk::FluxFn
    1473                 :            :     static tk::FluxFn::result_type
    1474                 :    7374000 :     flux( ncomp_t ncomp,
    1475                 :            :           const std::vector< EOS >& mat_blk,
    1476                 :            :           const std::vector< tk::real >& ugp,
    1477                 :            :           const std::vector< std::array< tk::real, 3 > >& )
    1478                 :            :     {
    1479                 :    7374000 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
    1480                 :            : 
    1481                 :    7374000 :       return tk::fluxTerms(ncomp, nmat, mat_blk, ugp);
    1482                 :            :     }
    1483                 :            : 
    1484                 :            :     //! \brief Boundary state function providing the left and right state of a
    1485                 :            :     //!   face at Dirichlet boundaries
    1486                 :            :     //! \param[in] ncomp Number of scalar components in this PDE system
    1487                 :            :     //! \param[in] mat_blk EOS material block
    1488                 :            :     //! \param[in] ul Left (domain-internal) state
    1489                 :            :     //! \param[in] x X-coordinate at which to compute the states
    1490                 :            :     //! \param[in] y Y-coordinate at which to compute the states
    1491                 :            :     //! \param[in] z Z-coordinate at which to compute the states
    1492                 :            :     //! \param[in] t Physical time
    1493                 :            :     //! \return Left and right states for all scalar components in this PDE
    1494                 :            :     //!   system
    1495                 :            :     //! \note The function signature must follow tk::StateFn. For multimat, the
    1496                 :            :     //!   left or right state is the vector of conserved quantities, followed by
    1497                 :            :     //!   the vector of primitive quantities appended to it.
    1498                 :            :     static tk::StateFn::result_type
    1499                 :      48000 :     dirichlet( ncomp_t ncomp,
    1500                 :            :                const std::vector< EOS >& mat_blk,
    1501                 :            :                const std::vector< tk::real >& ul, tk::real x, tk::real y,
    1502                 :            :                tk::real z, tk::real t, const std::array< tk::real, 3 >& )
    1503                 :            :     {
    1504                 :      48000 :       auto nmat = g_inputdeck.get< tag::multimat, tag::nmat >();
    1505                 :            :       const auto& solidx = g_inputdeck.get<
    1506                 :      48000 :         tag::matidxmap, tag::solidx >();
    1507                 :            : 
    1508         [ +  - ]:      48000 :       [[maybe_unused]] auto nsld = numSolids(nmat, solidx);
    1509                 :            : 
    1510         [ +  - ]:      96000 :       auto ur = Problem::initialize( ncomp, mat_blk, x, y, z, t );
    1511 [ -  + ][ -  - ]:      48000 :       Assert( ur.size() == ncomp, "Incorrect size for boundary state vector" );
         [ -  - ][ -  - ]
    1512                 :            : 
    1513         [ +  - ]:      48000 :       ur.resize(ul.size());
    1514                 :            : 
    1515                 :      48000 :       tk::real rho(0.0);
    1516         [ +  + ]:     192000 :       for (std::size_t k=0; k<nmat; ++k)
    1517                 :     144000 :         rho += ur[densityIdx(nmat, k)];
    1518                 :            : 
    1519                 :            :       // get primitives in boundary state
    1520                 :            : 
    1521                 :            :       // velocity
    1522                 :      48000 :       ur[ncomp+velocityIdx(nmat, 0)] = ur[momentumIdx(nmat, 0)] / rho;
    1523                 :      48000 :       ur[ncomp+velocityIdx(nmat, 1)] = ur[momentumIdx(nmat, 1)] / rho;
    1524                 :      48000 :       ur[ncomp+velocityIdx(nmat, 2)] = ur[momentumIdx(nmat, 2)] / rho;
    1525                 :            : 
    1526                 :            :       // material pressures
    1527         [ +  + ]:     192000 :       for (std::size_t k=0; k<nmat; ++k)
    1528                 :            :       {
    1529         [ +  - ]:     144000 :         auto gk = getDeformGrad(nmat, k, ur);
    1530         [ +  - ]:     288000 :         ur[ncomp+pressureIdx(nmat, k)] = mat_blk[k].compute< EOS::pressure >(
    1531                 :     144000 :           ur[densityIdx(nmat, k)], ur[ncomp+velocityIdx(nmat, 0)],
    1532                 :     144000 :           ur[ncomp+velocityIdx(nmat, 1)], ur[ncomp+velocityIdx(nmat, 2)],
    1533                 :     144000 :           ur[energyIdx(nmat, k)], ur[volfracIdx(nmat, k)], k, gk );
    1534                 :            :       }
    1535                 :            : 
    1536 [ -  + ][ -  - ]:      48000 :       Assert( ur.size() == ncomp+nmat+3+nsld*6, "Incorrect size for appended "
         [ -  - ][ -  - ]
    1537                 :            :               "boundary state vector" );
    1538                 :            : 
    1539         [ +  - ]:      96000 :       return {{ std::move(ul), std::move(ur) }};
    1540                 :            :     }
    1541                 :            : 
    1542                 :            :     // Other boundary condition types that do not depend on "Problem" should be
    1543                 :            :     // added in BCFunctions.hpp
    1544                 :            : };
    1545                 :            : 
    1546                 :            : } // dg::
    1547                 :            : 
    1548                 :            : } // inciter::
    1549                 :            : 
    1550                 :            : #endif // DGMultiMat_h

Generated by: LCOV version 1.14