Quinoa all test code coverage report
Current view: top level - PDE/Transport - DGTransport.hpp (source / functions) Hit Total Coverage
Commit: -128-NOTFOUND Lines: 70 83 84.3 %
Date: 2025-02-13 15:35:21 Functions: 32 75 42.7 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 68 170 40.0 %

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/PDE/Transport/DGTransport.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     Scalar transport using disccontinous Galerkin discretization
       9                 :            :   \details   This file implements the physics operators governing transported
      10                 :            :      scalars using disccontinuous Galerkin discretization.
      11                 :            : */
      12                 :            : // *****************************************************************************
      13                 :            : #ifndef DGTransport_h
      14                 :            : #define DGTransport_h
      15                 :            : 
      16                 :            : #include <vector>
      17                 :            : #include <array>
      18                 :            : #include <limits>
      19                 :            : #include <cmath>
      20                 :            : #include <unordered_set>
      21                 :            : #include <map>
      22                 :            : 
      23                 :            : #include "Macro.hpp"
      24                 :            : #include "Exception.hpp"
      25                 :            : #include "Vector.hpp"
      26                 :            : #include "UnsMesh.hpp"
      27                 :            : #include "Integrate/Basis.hpp"
      28                 :            : #include "Integrate/Quadrature.hpp"
      29                 :            : #include "Integrate/Initialize.hpp"
      30                 :            : #include "Integrate/Mass.hpp"
      31                 :            : #include "Integrate/Surface.hpp"
      32                 :            : #include "Integrate/Boundary.hpp"
      33                 :            : #include "Integrate/Volume.hpp"
      34                 :            : #include "Riemann/Upwind.hpp"
      35                 :            : #include "Reconstruction.hpp"
      36                 :            : #include "Limiter.hpp"
      37                 :            : #include "PrefIndicator.hpp"
      38                 :            : #include "EoS/EOS.hpp"
      39                 :            : #include "FunctionPrototypes.hpp"
      40                 :            : #include "ConfigureTransport.hpp"
      41                 :            : 
      42                 :            : namespace inciter {
      43                 :            : 
      44                 :            : extern ctr::InputDeck g_inputdeck;
      45                 :            : 
      46                 :            : namespace dg {
      47                 :            : 
      48                 :            : //! \brief Transport equation used polymorphically with tk::DGPDE
      49                 :            : //! \details The template argument(s) specify policies and are used to configure
      50                 :            : //!   the behavior of the class. The policies are:
      51                 :            : //!   - Physics - physics configuration, see PDE/Transport/Physics.h
      52                 :            : //!   - Problem - problem configuration, see PDE/Transport/Problem.h
      53                 :            : //! \note The default physics is DGAdvection, set in
      54                 :            : //!    inciter::deck::check_transport()
      55                 :            : template< class Physics, class Problem >
      56 [ -  - ][ +  - ]:        296 : class Transport {
         [ +  - ][ +  - ]
                 [ -  - ]
      57                 :            : 
      58                 :            :   private:
      59                 :            :     using eq = tag::transport;
      60                 :            : 
      61                 :            :   public:
      62                 :            :     //! Constructor
      63                 :        148 :     explicit Transport() :
      64                 :            :       m_physics( Physics() ),
      65                 :            :       m_problem( Problem() ),
      66         [ +  - ]:        148 :       m_ncomp( g_inputdeck.get< tag::ncomp >() )
      67                 :            :     {
      68                 :            :       // associate boundary condition configurations with state functions, the
      69                 :            :       // order in which the state functions listed matters, see ctr::bc::Keys
      70 [ +  - ][ +  - ]:       1924 :       brigand::for_each< ctr::bclist::Keys >( ConfigBC( m_bc,
         [ +  - ][ +  + ]
         [ +  - ][ +  + ]
         [ +  - ][ +  - ]
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
      71                 :            :         // BC State functions
      72                 :            :         { dirichlet
      73                 :            :         , invalidBC  // Symmetry BC not implemented
      74                 :            :         , outlet
      75                 :            :         , invalidBC  // Characteristic BC not implemented
      76                 :            :         , extrapolate
      77                 :            :         , invalidBC },      // No slip wall BC not implemented
      78                 :            :         // BC Gradient functions
      79                 :            :         { noOpGrad
      80                 :            :         , noOpGrad
      81                 :            :         , noOpGrad
      82                 :            :         , noOpGrad
      83                 :            :         , noOpGrad
      84                 :            :         , noOpGrad }
      85                 :            :         ) );
      86                 :            : 
      87                 :            :       // Inlet BC has a different structure than above BCs, so it must be 
      88                 :            :       // handled differently than with ConfigBC
      89 [ +  - ][ +  - ]:        296 :       ConfigInletBC(m_bc, inlet, noOpGrad);
                 [ -  - ]
      90                 :            : 
      91         [ -  - ]:          0 :       m_problem.errchk( m_ncomp );
      92                 :        148 :     }
      93                 :            : 
      94                 :            :     //! Find the number of primitive quantities required for this PDE system
      95                 :            :     //! \return The number of primitive quantities required to be stored for
      96                 :            :     //!   this PDE system
      97                 :            :     std::size_t nprim() const
      98                 :            :     {
      99                 :            :       // transport does not need/store any primitive quantities currently
     100                 :            :       return 0;
     101                 :            :     }
     102                 :            : 
     103                 :            :     //! Find the number of materials set up for this PDE system
     104                 :            :     //! \return The number of materials set up for this PDE system
     105                 :            :     std::size_t nmat() const
     106                 :            :     {
     107                 :            :       return m_ncomp;
     108                 :            :     }
     109                 :            : 
     110                 :            :     //! Assign number of DOFs per equation in the PDE system
     111                 :            :     //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
     112                 :            :     //!   equation
     113                 :            :     void numEquationDofs(std::vector< std::size_t >& numEqDof) const
     114                 :            :     {
     115                 :            :       // all equation-dofs initialized to ndofs
     116 [ -  - ][ +  + ]:       1060 :       for (std::size_t i=0; i<m_ncomp; ++i) {
         [ +  + ][ +  + ]
                 [ -  - ]
     117                 :        530 :         numEqDof.push_back(g_inputdeck.get< tag::ndof >());
     118                 :            :       }
     119                 :            :     }
     120                 :            : 
     121                 :            :     //! Find how 'stiff equations', which we currently
     122                 :            :     //! have none for Transport
     123                 :            :     //! \return number of stiff equations
     124                 :            :     std::size_t nstiffeq() const
     125                 :            :     { return 0; }
     126                 :            : 
     127                 :            :     //! Find how 'nonstiff equations', which we currently
     128                 :            :     //! don't use for Transport
     129                 :            :     //! \return number of non-stiff equations
     130                 :            :     std::size_t nnonstiffeq() const
     131                 :            :     { return 0; }
     132                 :            : 
     133                 :            :     //! Locate the stiff equations. Unused for transport.
     134                 :            :     //! \param[out] stiffEqIdx list
     135                 :            :     void setStiffEqIdx( std::vector< std::size_t >& stiffEqIdx ) const
     136                 :            :     {
     137                 :          0 :       stiffEqIdx.resize(0);
     138                 :            :     }
     139                 :            : 
     140                 :            :     //! Locate the nonstiff equations. Unused for transport.
     141                 :            :     //! \param[out] nonStiffEqIdx list
     142                 :            :     void setNonStiffEqIdx( std::vector< std::size_t >& nonStiffEqIdx ) const
     143                 :            :     {
     144                 :          0 :       nonStiffEqIdx.resize(0);
     145                 :            :     }
     146                 :            : 
     147                 :            :     //! Determine elements that lie inside the user-defined IC box
     148                 :            :     void IcBoxElems( const tk::Fields&,
     149                 :            :       std::size_t,
     150                 :            :       std::vector< std::unordered_set< std::size_t > >& ) const
     151                 :            :     {}
     152                 :            : 
     153                 :            :     //! Initalize the transport equations for DG
     154                 :            :     //! \param[in] L Element mass matrix
     155                 :            :     //! \param[in] inpoel Element-node connectivity
     156                 :            :     //! \param[in] coord Array of nodal coordinates
     157                 :            :     //! \param[in,out] unk Array of unknowns
     158                 :            :     //! \param[in] t Physical time
     159                 :            :     //! \param[in] nielem Number of internal elements
     160                 :            :     void
     161         [ +  - ]:        654 :     initialize(
     162                 :            :       const tk::Fields& L,
     163                 :            :       const std::vector< std::size_t >& inpoel,
     164                 :            :       const tk::UnsMesh::Coords& coord,
     165                 :            :       const std::vector< std::unordered_set< std::size_t > >& /*inbox*/,
     166                 :            :       const std::unordered_map< std::size_t, std::set< std::size_t > >&,
     167                 :            :       tk::Fields& unk,
     168                 :            :       tk::real t,
     169                 :            :       const std::size_t nielem ) const
     170                 :            :     {
     171         [ +  - ]:        654 :       tk::initialize( m_ncomp, m_mat_blk, L, inpoel, coord,
     172                 :            :                       Problem::initialize, unk, t, nielem );
     173                 :        654 :     }
     174                 :            : 
     175                 :            :     //! Compute density constraint for a given material
     176                 :            :     // //! \param[in] nelem Number of elements
     177                 :            :     // //! \param[in] unk Array of unknowns
     178                 :            :     //! \param[out] densityConstr Density Constraint: rho/(rho0*det(g))
     179                 :            :     void computeDensityConstr( std::size_t /*nelem*/,
     180                 :            :                                tk::Fields& /*unk*/,
     181                 :            :                                std::vector< tk::real >& densityConstr) const
     182                 :            :     {
     183                 :        993 :       densityConstr.resize(0);
     184                 :            :     }
     185                 :            : 
     186                 :            :     //! Compute the left hand side mass matrix
     187                 :            :     //! \param[in] geoElem Element geometry array
     188                 :            :     //! \param[in,out] l Block diagonal mass matrix
     189                 :            :     void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
     190                 :        654 :       const auto ndof = g_inputdeck.get< tag::ndof >();
     191                 :        654 :       tk::mass( m_ncomp, ndof, geoElem, l );
     192                 :            :     }
     193                 :            : 
     194                 :            :     //! Update the interface cells to first order dofs
     195                 :            :     //! \details This function resets the high-order terms in interface cells,
     196                 :            :     //!   and is currently not used in transport.
     197                 :            :     void updateInterfaceCells( tk::Fields&,
     198                 :            :       std::size_t,
     199                 :            :       std::vector< std::size_t >&,
     200                 :            :       std::vector< std::size_t >& ) const {}
     201                 :            : 
     202                 :            :     //! Update the primitives for this PDE system
     203                 :            :     //! \details This function computes and stores the dofs for primitive
     204                 :            :     //!   quantities, which are currently unused for transport.
     205                 :            :     void updatePrimitives( const tk::Fields&,
     206                 :            :                            const tk::Fields&,
     207                 :            :                            const tk::Fields&,
     208                 :            :                            tk::Fields&,
     209                 :            :                            std::size_t,
     210                 :            :                            std::vector< std::size_t >& ) const {}
     211                 :            : 
     212                 :            :     //! Clean up the state of trace materials for this PDE system
     213                 :            :     //! \details This function cleans up the state of materials present in trace
     214                 :            :     //!   quantities in each cell. This is currently unused for transport.
     215                 :            :     void cleanTraceMaterial( tk::real,
     216                 :            :                              const tk::Fields&,
     217                 :            :                              tk::Fields&,
     218                 :            :                              tk::Fields&,
     219                 :            :                              std::size_t ) const {}
     220                 :            : 
     221                 :            :     //! Reconstruct second-order solution from first-order
     222                 :            : //    //! \param[in] t Physical time
     223                 :            : //    //! \param[in] geoFace Face geometry array
     224                 :            : //    //! \param[in] geoElem Element geometry array
     225                 :            : //    //! \param[in] fd Face connectivity and boundary conditions object
     226                 :            : //    //! \param[in] esup Elements-surrounding-nodes connectivity
     227                 :            : //    //! \param[in] inpoel Element-node connectivity
     228                 :            : //    //! \param[in] coord Array of nodal coordinates
     229                 :            : //    //! \param[in,out] U Solution vector at recent time step
     230                 :            : //    //! \param[in,out] P Primitive vector at recent time step
     231                 :      19125 :     void reconstruct( tk::real,
     232                 :            :                       const tk::Fields&,
     233                 :            :                       const tk::Fields&,
     234                 :            :                       const inciter::FaceData&,
     235                 :            :                       const std::map< std::size_t, std::vector< std::size_t > >&,
     236                 :            :                       const std::vector< std::size_t >&,
     237                 :            :                       const tk::UnsMesh::Coords&,
     238                 :            :                       tk::Fields&,
     239                 :            :                       tk::Fields&,
     240                 :            :                       const bool,
     241                 :            :                       const std::vector< std::size_t >& ) const
     242                 :            :     {
     243                 :            :       // do reconstruction only if P0P1
     244         [ +  + ]:      19125 :       if (g_inputdeck.get< tag::rdof >() == 4 &&
     245         [ -  + ]:      18750 :         g_inputdeck.get< tag::ndof >() == 1)
     246 [ -  - ][ -  - ]:          0 :         Throw("P0P1 not supported for Transport.");
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
     247                 :      19125 :     }
     248                 :            : 
     249                 :            :     //! Limit second-order solution
     250                 :            :     //! \param[in] t Physical time
     251                 :            :     //! \param[in] geoFace Face geometry array
     252                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     253                 :            :     //! \param[in] esup Elements surrounding points
     254                 :            :     //! \param[in] inpoel Element-node connectivity
     255                 :            :     //! \param[in] coord Array of nodal coordinates
     256                 :            :     //! \param[in] ndofel Vector of local number of degrees of freedome
     257                 :            : //    //! \param[in] gid Local->global node id map
     258                 :            : //    //! \param[in] bid Local chare-boundary node ids (value) associated to
     259                 :            : //    //!   global node ids (key)
     260                 :            : //    //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
     261                 :            : //    //!   variables
     262                 :            :     //! \param[in,out] U Solution vector at recent time step
     263                 :      19125 :     void limit( [[maybe_unused]] tk::real t,
     264                 :            :                 [[maybe_unused]] const bool pref,
     265                 :            :                 [[maybe_unused]] const tk::Fields& geoFace,
     266                 :            :                 const tk::Fields&,
     267                 :            :                 const inciter::FaceData& fd,
     268                 :            :                 const std::map< std::size_t, std::vector< std::size_t > >& esup,
     269                 :            :                 const std::vector< std::size_t >& inpoel,
     270                 :            :                 const tk::UnsMesh::Coords& coord,
     271                 :            :                 const std::vector< std::size_t >& ndofel,
     272                 :            :                 const std::vector< std::size_t >&,
     273                 :            :                 const std::unordered_map< std::size_t, std::size_t >&,
     274                 :            :                 const std::vector< std::vector<tk::real> >&,
     275                 :            :                 const std::vector< std::vector<tk::real> >&,
     276                 :            :                 const std::vector< std::vector<tk::real> >&,
     277                 :            :                 tk::Fields& U,
     278                 :            :                 tk::Fields&,
     279                 :            :                 std::vector< std::size_t >& ) const
     280                 :            :     {
     281                 :      19125 :       const auto limiter = g_inputdeck.get< tag::limiter >();
     282                 :            : 
     283         [ +  + ]:      19125 :       if (limiter == ctr::LimiterType::WENOP1)
     284                 :       6000 :         WENO_P1( fd.Esuel(), U );
     285         [ +  + ]:      13125 :       else if (limiter == ctr::LimiterType::SUPERBEEP1)
     286                 :      12750 :         Superbee_P1( fd.Esuel(), inpoel, ndofel, coord, U );
     287         [ -  + ]:        375 :       else if (limiter == ctr::LimiterType::VERTEXBASEDP1)
     288                 :          0 :         VertexBasedTransport_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
     289                 :            :           coord, U );
     290                 :      19125 :     }
     291                 :            : 
     292                 :            :     //! Update the conservative variable solution for this PDE system
     293                 :            :     //! \details This function computes the updated dofs for conservative
     294                 :            :     //!   quantities based on the limited solution and is currently not used in
     295                 :            :     //!   transport.
     296                 :            :     void CPL( const tk::Fields&,
     297                 :            :               const tk::Fields&,
     298                 :            :               const std::vector< std::size_t >&,
     299                 :            :               const tk::UnsMesh::Coords&,
     300                 :            :               tk::Fields&,
     301                 :            :               std::size_t ) const {}
     302                 :            : 
     303                 :            :     //! Return cell-average deformation gradient tensor (no-op for transport)
     304                 :            :     //! \details This function is a no-op in transport.
     305                 :            :     std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
     306                 :            :       const tk::Fields&,
     307                 :            :       std::size_t ) const
     308                 :            :     {
     309                 :            :       return {};
     310                 :            :     }
     311                 :            : 
     312                 :            :     //! Reset the high order solution for p-adaptive scheme
     313                 :            :     //! \details This function reset the high order coefficient for p-adaptive
     314                 :            :     //!   solution polynomials and is currently not used in transport.
     315                 :            :     void resetAdapSol( const inciter::FaceData&,
     316                 :            :                        tk::Fields&,
     317                 :            :                        tk::Fields&,
     318                 :            :                        const std::vector< std::size_t >& ) const {}
     319                 :            : 
     320                 :            :     //! Compute right hand side
     321                 :            :     //! \param[in] t Physical time
     322                 :            :     //! \param[in] pref Indicator for p-adaptive algorithm
     323                 :            :     //! \param[in] geoFace Face geometry array
     324                 :            :     //! \param[in] geoElem Element geometry array
     325                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     326                 :            :     //! \param[in] inpoel Element-node connectivity
     327                 :            :     //! \param[in] coord Array of nodal coordinates
     328                 :            :     //! \param[in] U Solution vector at recent time step
     329                 :            :     //! \param[in] P Primitive vector at recent time step
     330                 :            :     //! \param[in] ndofel Vector of local number of degrees of freedom
     331                 :            :     //! \param[in] dt Delta time
     332                 :            :     //! \param[in,out] R Right-hand side vector computed
     333                 :      44280 :     void rhs( tk::real t,
     334                 :            :               const bool pref,
     335                 :            :               const tk::Fields& geoFace,
     336                 :            :               const tk::Fields& geoElem,
     337                 :            :               const inciter::FaceData& fd,
     338                 :            :               const std::vector< std::size_t >& inpoel,
     339                 :            :               const std::vector< std::unordered_set< std::size_t > >&,
     340                 :            :               const tk::UnsMesh::Coords& coord,
     341                 :            :               const tk::Fields& U,
     342                 :            :               const tk::Fields& P,
     343                 :            :               const std::vector< std::size_t >& ndofel,
     344                 :            :               const tk::real dt,
     345                 :            :               tk::Fields& R ) const
     346                 :            :     {
     347                 :      44280 :       const auto ndof = g_inputdeck.get< tag::ndof >();
     348                 :      44280 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     349                 :      44280 :       const auto intsharp = g_inputdeck.get< tag::transport,
     350                 :            :         tag::intsharp >();
     351                 :            : 
     352                 :            :       Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
     353                 :            :               "vector and primitive vector at recent time step incorrect" );
     354                 :            :       Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
     355                 :            :               "vector and right-hand side at recent time step incorrect" );
     356                 :            :       Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
     357                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     358                 :            :       Assert( P.nprop() == 0, "Number of components in primitive "
     359                 :            :               "vector must equal "+ std::to_string(0) );
     360                 :            :       Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
     361                 :            :               "side vector must equal "+ std::to_string(ndof*m_ncomp) );
     362                 :            :       Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
     363                 :            :               "Mismatch in inpofa size" );
     364                 :            : 
     365                 :            :       // set rhs to zero
     366                 :            :       R.fill(0.0);
     367                 :            : 
     368                 :            :       // empty vector for non-conservative terms. This vector is unused for
     369                 :            :       // linear transport since, there are no non-conservative terms in the
     370                 :            :       // system of PDEs.
     371                 :      44280 :       std::vector< std::vector < tk::real > > riemannDeriv;
     372                 :            : 
     373                 :            :       // compute internal surface flux integrals
     374 [ +  - ][ +  - ]:      44280 :       std::vector< std::size_t > solidx(1, 0);
     375 [ +  - ][ +  - ]:      88560 :       tk::surfInt( pref, m_ncomp, m_mat_blk, t, ndof, rdof,
                 [ -  - ]
     376                 :            :                    inpoel, solidx, coord, fd, geoFace, geoElem, Upwind::flux,
     377                 :            :                    Problem::prescribedVelocity, U, P, ndofel, dt, R,
     378                 :            :                    riemannDeriv, intsharp );
     379                 :            : 
     380         [ +  + ]:      44280 :       if(ndof > 1)
     381                 :            :         // compute volume integrals
     382 [ +  - ][ +  - ]:      57375 :         tk::volInt( m_ncomp, t, m_mat_blk, ndof, rdof,
                 [ -  - ]
     383         [ +  - ]:      19125 :                     fd.Esuel().size()/4, inpoel, coord, geoElem, flux,
     384                 :            :                     Problem::prescribedVelocity, U, P, ndofel, R, intsharp );
     385                 :            : 
     386                 :            :       // compute boundary surface flux integrals
     387         [ +  + ]:     324840 :       for (const auto& b : m_bc)
     388 [ +  - ][ +  - ]:     841680 :         tk::bndSurfInt( pref, m_ncomp, m_mat_blk, ndof, rdof,
         [ -  - ][ -  - ]
     389                 :            :           std::get<0>(b), fd, geoFace, geoElem, inpoel, coord, t, Upwind::flux,
     390                 :            :           Problem::prescribedVelocity, std::get<1>(b), U, P, ndofel, R,
     391                 :            :           riemannDeriv, intsharp );
     392                 :      44280 :     }
     393                 :            : 
     394                 :            :     //! Evaluate the adaptive indicator and mark the ndof for each element
     395                 :            :     //! \param[in] nunk Number of unknowns
     396                 :            :     //! \param[in] coord Array of nodal coordinates
     397                 :            :     //! \param[in] inpoel Element-node connectivity
     398                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     399                 :            :     //! \param[in] unk Array of unknowns
     400                 :            :     //! \param[in] prim Array of primitive quantities
     401                 :            :     //! \param[in] indicator p-refinement indicator type
     402                 :            :     //! \param[in] ndof Number of degrees of freedom in the solution
     403                 :            :     //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
     404                 :            :     //! \param[in] tolref Tolerance for p-refinement
     405                 :            :     //! \param[in,out] ndofel Vector of local number of degrees of freedome
     406         [ -  - ]:          0 :     void eval_ndof( std::size_t nunk,
     407                 :            :                     [[maybe_unused]] const tk::UnsMesh::Coords& coord,
     408                 :            :                     [[maybe_unused]] const std::vector< std::size_t >& inpoel,
     409                 :            :                     const inciter::FaceData& fd,
     410                 :            :                     const tk::Fields& unk,
     411                 :            :                     [[maybe_unused]] const tk::Fields& prim,
     412                 :            :                     inciter::ctr::PrefIndicatorType indicator,
     413                 :            :                     std::size_t ndof,
     414                 :            :                     std::size_t ndofmax,
     415                 :            :                     tk::real tolref,
     416                 :            :                     std::vector< std::size_t >& ndofel ) const
     417                 :            :     {
     418                 :            :       const auto& esuel = fd.Esuel();
     419                 :            : 
     420         [ -  - ]:          0 :       if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
     421                 :          0 :         spectral_decay( 1, nunk, esuel, unk, ndof, ndofmax, tolref, ndofel );
     422                 :            :       else
     423 [ -  - ][ -  - ]:          0 :         Throw( "No such adaptive indicator type" );
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
     424                 :          0 :     }
     425                 :            : 
     426                 :            :     //! Compute the minimum time step size
     427                 :            : //     //! \param[in] U Solution vector at recent time step
     428                 :            : //     //! \param[in] coord Mesh node coordinates
     429                 :            : //     //! \param[in] inpoel Mesh element connectivity
     430                 :            :     //! \return Minimum time step size
     431                 :            :     tk::real dt( const std::array< std::vector< tk::real >, 3 >& /*coord*/,
     432                 :            :                  const std::vector< std::size_t >& /*inpoel*/,
     433                 :            :                  const inciter::FaceData& /*fd*/,
     434                 :            :                  const tk::Fields& /*geoFace*/,
     435                 :            :                  const tk::Fields& /*geoElem*/,
     436                 :            :                  const std::vector< std::size_t >& /*ndofel*/,
     437                 :            :                  const tk::Fields& /*U*/,
     438                 :            :                  const tk::Fields&,
     439                 :            :                  const std::size_t /*nielem*/ ) const
     440                 :            :     {
     441                 :            :       tk::real mindt = std::numeric_limits< tk::real >::max();
     442                 :            :       return mindt;
     443                 :            :     }
     444                 :            : 
     445                 :            :     //! Compute stiff terms for a single element, not implemented here
     446                 :            :     // //! \param[in] e Element number
     447                 :            :     // //! \param[in] geoElem Element geometry array
     448                 :            :     // //! \param[in] inpoel Element-node connectivity
     449                 :            :     // //! \param[in] coord Array of nodal coordinates
     450                 :            :     // //! \param[in] U Solution vector at recent time step
     451                 :            :     // //! \param[in] P Primitive vector at recent time step
     452                 :            :     // //! \param[in] ndofel Vector of local number of degrees of freedom
     453                 :            :     // //! \param[in,out] R Right-hand side vector computed
     454                 :            :     void stiff_rhs( std::size_t /*e*/,
     455                 :            :                     const tk::Fields& /*geoElem*/,
     456                 :            :                     const std::vector< std::size_t >& /*inpoel*/,
     457                 :            :                     const tk::UnsMesh::Coords& /*coord*/,
     458                 :            :                     const tk::Fields& /*U*/,
     459                 :            :                     const tk::Fields& /*P*/,
     460                 :            :                     const std::vector< std::size_t >& /*ndofel*/,
     461                 :            :                     tk::Fields& /*R*/ ) const
     462                 :            :     {}
     463                 :            : 
     464                 :            :     //! Return a map that associates user-specified strings to functions
     465                 :            :     //! \return Map that associates user-specified strings to functions that
     466                 :            :     //!  compute relevant quantities to be output to file
     467         [ +  - ]:       1986 :     std::map< std::string, tk::GetVarFn > OutVarFn() const {
     468                 :            :       std::map< std::string, tk::GetVarFn > OutFnMap;
     469 [ +  - ][ +  - ]:       1986 :       OutFnMap["material_indicator"] = transport::matIndicatorOutVar;
     470                 :            : 
     471                 :       1986 :       return OutFnMap;
     472                 :            :     }
     473                 :            : 
     474                 :            :     //! Return analytic field names to be output to file
     475                 :            :     //! \return Vector of strings labelling analytic fields output in file
     476                 :        917 :     std::vector< std::string > analyticFieldNames() const {
     477                 :            :       std::vector< std::string > n;
     478                 :        917 :       auto depvar = g_inputdeck.get< tag::depvar >()[0];
     479         [ +  + ]:       1834 :       for (ncomp_t c=0; c<m_ncomp; ++c)
     480 [ +  - ][ +  - ]:       1834 :         n.push_back( depvar + std::to_string(c) + "_analytic" );
         [ -  + ][ -  + ]
         [ -  - ][ -  - ]
     481                 :        917 :       return n;
     482                 :            :     }
     483                 :            : 
     484                 :            :     //! Return surface field output going to file
     485                 :            :     std::vector< std::vector< tk::real > >
     486                 :            :     surfOutput( const std::map< int, std::vector< std::size_t > >&,
     487                 :            :                 tk::Fields& ) const
     488                 :            :     {
     489                 :            :       std::vector< std::vector< tk::real > > s; // punt for now
     490                 :            :       return s;
     491                 :            :     }
     492                 :            : 
     493                 :            :     //! Return time history field names to be output to file
     494                 :            :     //! \return Vector of strings labelling time history fields output in file
     495                 :            :     std::vector< std::string > histNames() const {
     496                 :            :       std::vector< std::string > s; // punt for now
     497                 :            :       return s;
     498                 :            :     }
     499                 :            : 
     500                 :            :     //! Return names of integral variables to be output to diagnostics file
     501                 :            :     //! \return Vector of strings labelling integral variables output
     502         [ +  - ]:         38 :     std::vector< std::string > names() const {
     503                 :            :       std::vector< std::string > n;
     504                 :            :       const auto& depvar =
     505                 :            :       g_inputdeck.get< tag::depvar >().at(0);
     506                 :            :       // construct the name of the numerical solution for all components
     507         [ +  + ]:         76 :       for (ncomp_t c=0; c<m_ncomp; ++c)
     508 [ +  - ][ -  + ]:         76 :         n.push_back( depvar + std::to_string(c) );
                 [ -  - ]
     509                 :         38 :       return n;
     510                 :            :     }
     511                 :            : 
     512                 :            :     //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
     513                 :            :     //! \param[in] xi X-coordinate at which to evaluate the analytic solution
     514                 :            :     //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
     515                 :            :     //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
     516                 :            :     //! \param[in] t Physical time at which to evaluate the analytic solution
     517                 :            :     //! \return Vector of analytic solution at given spatial location and time
     518                 :            :     std::vector< tk::real >
     519                 :            :     analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
     520                 :     839522 :     { return Problem::analyticSolution( m_ncomp, m_mat_blk, xi, yi,
     521                 :            :                                         zi, t ); }
     522                 :            : 
     523                 :            :     //! Return analytic solution for conserved variables
     524                 :            :     //! \param[in] xi X-coordinate at which to evaluate the analytic solution
     525                 :            :     //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
     526                 :            :     //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
     527                 :            :     //! \param[in] t Physical time at which to evaluate the analytic solution
     528                 :            :     //! \return Vector of analytic solution at given location and time
     529                 :            :     std::vector< tk::real >
     530                 :            :     solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
     531                 :    2211605 :     { return Problem::initialize( m_ncomp, m_mat_blk, xi, yi, zi, t ); }
     532                 :            : 
     533                 :            :     //! Return time history field output evaluated at time history points
     534                 :            :     //! \param[in] h History point data
     535                 :            :     std::vector< std::vector< tk::real > >
     536                 :            :     histOutput( const std::vector< HistData >& h,
     537                 :            :                 const std::vector< std::size_t >&,
     538                 :            :                 const tk::UnsMesh::Coords&,
     539                 :            :                 const tk::Fields&,
     540                 :            :                 const tk::Fields& ) const
     541                 :            :     {
     542                 :          0 :       std::vector< std::vector< tk::real > > Up(h.size()); //punt for now
     543                 :            :       return Up;
     544                 :            :     }
     545                 :            : 
     546                 :            :     //! Return cell-averaged total component mass per unit volume for an element
     547                 :            :     //! \param[in] e Element id for which total energy is required
     548                 :            :     //! \param[in] unk Vector of conserved quantities
     549                 :            :     //! \return Cell-averaged total component mass per unit volume for given
     550                 :            :     //!   element. Since transport does not have an associated total energy,
     551                 :            :     //!   return total mass.
     552                 :            :     tk::real sp_totalenergy(std::size_t e, const tk::Fields& unk) const
     553                 :            :     {
     554                 :    1315427 :       const auto rdof = g_inputdeck.get< tag::rdof >();
     555                 :            : 
     556                 :            :       tk::real sp_m(0.0);
     557 [ -  - ][ +  + ]:    2630854 :       for (std::size_t c=0; c<m_ncomp; ++c) {
         [ +  + ][ +  + ]
                 [ -  - ]
     558                 :    1315427 :         sp_m += unk(e,c*rdof);
     559                 :            :       }
     560                 :            :       return sp_m;
     561                 :            :     }
     562                 :            : 
     563                 :            :   private:
     564                 :            :     const Physics m_physics;            //!< Physics policy
     565                 :            :     const Problem m_problem;            //!< Problem policy
     566                 :            :     const ncomp_t m_ncomp;              //!< Number of components in this PDE
     567                 :            :     //! BC configuration
     568                 :            :     BCStateFn m_bc;
     569                 :            :     //! \brief EOS material block - This PDE does not require an EOS block,
     570                 :            :     //! thus this variable has not been intialized.
     571                 :            :     std::vector< EOS > m_mat_blk;
     572                 :            : 
     573                 :            :     //! Evaluate physical flux function for this PDE system
     574                 :            :     //! \param[in] ncomp Number of scalar components in this PDE system
     575                 :            :     //! \param[in] ugp Numerical solution at the Gauss point at which to
     576                 :            :     //!   evaluate the flux
     577                 :            :     //! \param[in] v Prescribed velocity evaluated at the Gauss point at which
     578                 :            :     //!   to evaluate the flux
     579                 :            :     //! \return Flux vectors for all components in this PDE system
     580                 :            :     //! \note The function signature must follow tk::FluxFn
     581                 :            :     static tk::FluxFn::result_type
     582                 :   19672200 :     flux( ncomp_t ncomp,
     583                 :            :           const std::vector< EOS >&,
     584                 :            :           const std::vector< tk::real >& ugp,
     585                 :            :           const std::vector< std::array< tk::real, 3 > >& v )
     586                 :            : 
     587                 :            :     {
     588                 :            :       Assert( ugp.size() == ncomp, "Size mismatch" );
     589                 :            :       Assert( v.size() == ncomp, "Size mismatch" );
     590                 :            : 
     591                 :   19672200 :       std::vector< std::array< tk::real, 3 > > fl( ugp.size() );
     592                 :            : 
     593         [ +  + ]:   39344400 :       for (ncomp_t c=0; c<ncomp; ++c)
     594                 :   19672200 :         fl[c] = {{ v[c][0] * ugp[c], v[c][1] * ugp[c], v[c][2] * ugp[c] }};
     595                 :            : 
     596                 :   19672200 :       return fl;
     597                 :            :     }
     598                 :            : 
     599                 :            :     //! \brief Boundary state function providing the left and right state of a
     600                 :            :     //!   face at extrapolation boundaries
     601                 :            :     //! \param[in] ul Left (domain-internal) state
     602                 :            :     //! \return Left and right states for all scalar components in this PDE
     603                 :            :     //!   system
     604                 :            :     //! \note The function signature must follow tk::StateFn
     605                 :            :     static tk::StateFn::result_type
     606                 :    8834220 :     extrapolate( ncomp_t, const std::vector< EOS >&,
     607                 :            :                  const std::vector< tk::real >& ul, tk::real, tk::real,
     608                 :            :                  tk::real, tk::real, const std::array< tk::real, 3 >& )
     609                 :            :     {
     610                 :    8834220 :       return {{ ul, ul }};
     611                 :            :     }
     612                 :            : 
     613                 :            :     //! \brief Boundary state function providing the left and right state of a
     614                 :            :     //!   face at extrapolation boundaries
     615                 :            :     //! \param[in] ul Left (domain-internal) state
     616                 :            :     //! \return Left and right states for all scalar components in this PDE
     617                 :            :     //!   system
     618                 :            :     //! \note The function signature must follow tk::StateFn
     619                 :            :     static tk::StateFn::result_type
     620                 :     163200 :     inlet( ncomp_t, const std::vector< EOS >&,
     621                 :            :            const std::vector< tk::real >& ul, tk::real, tk::real, tk::real,
     622                 :            :            tk::real, const std::array< tk::real, 3 >& )
     623                 :            :     {
     624                 :     163200 :       auto ur = ul;
     625                 :            :       std::fill( begin(ur), end(ur), 0.0 );
     626         [ +  - ]:     163200 :       return {{ ul, std::move(ur) }};
     627                 :            :     }
     628                 :            : 
     629                 :            :     //! \brief Boundary state function providing the left and right state of a
     630                 :            :     //!   face at outlet boundaries
     631                 :            :     //! \param[in] ul Left (domain-internal) state
     632                 :            :     //! \return Left and right states for all scalar components in this PDE
     633                 :            :     //!   system
     634                 :            :     //! \note The function signature must follow tk::StateFn
     635                 :            :     static tk::StateFn::result_type
     636                 :     715200 :     outlet( ncomp_t, const std::vector< EOS >&,
     637                 :            :             const std::vector< tk::real >& ul, tk::real, tk::real, tk::real,
     638                 :            :             tk::real, const std::array< tk::real, 3 >& )
     639                 :            :     {
     640                 :     715200 :       return {{ ul, ul }};
     641                 :            :     }
     642                 :            : 
     643                 :            :     //! \brief Boundary state function providing the left and right state of a
     644                 :            :     //!   face at Dirichlet boundaries
     645                 :            :     //! \param[in] ncomp Number of scalar components in this PDE system
     646                 :            :     //! \param[in] ul Left (domain-internal) state
     647                 :            :     //! \param[in] x X-coordinate at which to compute the states
     648                 :            :     //! \param[in] y Y-coordinate at which to compute the states
     649                 :            :     //! \param[in] z Z-coordinate at which to compute the states
     650                 :            :     //! \param[in] t Physical time
     651                 :            :     //! \return Left and right states for all scalar components in this PDE
     652                 :            :     //!   system
     653                 :            :     //! \note The function signature must follow tk::StateFn
     654                 :            :     static tk::StateFn::result_type
     655                 :     972540 :     dirichlet( ncomp_t ncomp, 
     656                 :            :                const std::vector< EOS >& mat_blk,
     657                 :            :                const std::vector< tk::real >& ul, tk::real x, tk::real y,
     658                 :            :                tk::real z, tk::real t, const std::array< tk::real, 3 >& )
     659                 :            :     {
     660                 :     972540 :       return {{ ul, Problem::initialize( ncomp, mat_blk, x, y, z, t ) }};
     661                 :            :     }
     662                 :            : 
     663                 :            :   //----------------------------------------------------------------------------
     664                 :            :   // Boundary Gradient functions
     665                 :            :   //----------------------------------------------------------------------------
     666                 :            : 
     667                 :            :   //! \brief Boundary gradient function copying the left gradient to the right
     668                 :            :   //!   gradient at a face
     669                 :            :   //! \param[in] dul Left (domain-internal) state
     670                 :            :   //! \return Left and right states for all scalar components in this PDE
     671                 :            :   //!   system
     672                 :            :   //! \note The function signature must follow tk::StateFn.
     673                 :            :   static tk::StateFn::result_type
     674                 :          0 :   noOpGrad( ncomp_t,
     675                 :            :             const std::vector< EOS >&,
     676                 :            :             const std::vector< tk::real >& dul,
     677                 :            :             tk::real, tk::real, tk::real, tk::real,
     678                 :            :             const std::array< tk::real, 3 >& )
     679                 :            :   {
     680                 :          0 :     return {{ dul, dul }};
     681                 :            :   }
     682                 :            : };
     683                 :            : 
     684                 :            : } // dg::
     685                 :            : } // inciter::
     686                 :            : 
     687                 :            : #endif // DGTransport_h

Generated by: LCOV version 1.14