Quinoa all test code coverage report
Current view: top level - PDE/MultiMat - DGMultiMat.hpp (source / functions) Hit Total Coverage
Commit: -128-NOTFOUND Lines: 305 591 51.6 %
Date: 2025-10-21 17:21:13 Functions: 87 780 11.2 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 237 844 28.1 %

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

Generated by: LCOV version 1.14