Quinoa all test code coverage report
Current view: top level - PDE/Transport - DGTransport.hpp (source / functions) Hit Total Coverage
Commit: Quinoa_v0.3-957-gb4f0efae0 Lines: 85 95 89.5 %
Date: 2021-11-09 15:14:18 Functions: 29 65 44.6 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 71 168 42.3 %

           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 "Inciter/Options/BC.hpp"
      27                 :            : #include "UnsMesh.hpp"
      28                 :            : #include "Integrate/Basis.hpp"
      29                 :            : #include "Integrate/Quadrature.hpp"
      30                 :            : #include "Integrate/Initialize.hpp"
      31                 :            : #include "Integrate/Mass.hpp"
      32                 :            : #include "Integrate/Surface.hpp"
      33                 :            : #include "Integrate/Boundary.hpp"
      34                 :            : #include "Integrate/Volume.hpp"
      35                 :            : #include "Riemann/Upwind.hpp"
      36                 :            : #include "Reconstruction.hpp"
      37                 :            : #include "Limiter.hpp"
      38                 :            : #include "PrefIndicator.hpp"
      39                 :            : 
      40                 :            : namespace inciter {
      41                 :            : 
      42                 :            : extern ctr::InputDeck g_inputdeck;
      43                 :            : 
      44                 :            : namespace dg {
      45                 :            : 
      46                 :            : //! \brief Transport equation used polymorphically with tk::DGPDE
      47                 :            : //! \details The template argument(s) specify policies and are used to configure
      48                 :            : //!   the behavior of the class. The policies are:
      49                 :            : //!   - Physics - physics configuration, see PDE/Transport/Physics.h
      50                 :            : //!   - Problem - problem configuration, see PDE/Transport/Problem.h
      51                 :            : //! \note The default physics is DGAdvection, set in
      52                 :            : //!    inciter::deck::check_transport()
      53                 :            : template< class Physics, class Problem >
      54 [ -  - ][ -  - ]:        326 : class Transport {
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
         [ +  - ][ +  - ]
         [ +  - ][ -  - ]
      55                 :            : 
      56                 :            :   private:
      57                 :            :     using eq = tag::transport;
      58                 :            : 
      59                 :            :   public:
      60                 :            :     //! Constructor
      61                 :            :     //! \param[in] c Equation system index (among multiple systems configured)
      62                 :        163 :     explicit Transport( ncomp_t c ) :
      63                 :            :       m_physics( Physics() ),
      64                 :            :       m_problem( Problem() ),
      65                 :            :       m_system( c ),
      66                 :            :       m_ncomp(
      67                 :            :         g_inputdeck.get< tag::component >().get< eq >().at(c) ),
      68                 :            :       m_offset(
      69 [ -  + ][ +  - ]:        326 :         g_inputdeck.get< tag::component >().offset< eq >(c) )
      70                 :            :     {
      71                 :            :       // associate boundary condition configurations with state functions, the
      72                 :            :       // order in which the state functions listed matters, see ctr::bc::Keys
      73 [ +  - ][ +  - ]:       1141 :       brigand::for_each< ctr::bc::Keys >( ConfigBC< eq >( m_system, m_bc,
         [ +  + ][ +  - ]
         [ -  - ][ -  - ]
      74                 :            :         { dirichlet
      75                 :            :         , invalidBC  // Symmetry BC not implemented
      76                 :            :         , inlet
      77                 :            :         , outlet
      78                 :            :         , invalidBC  // Characteristic BC not implemented
      79                 :            :         , extrapolate } ) );
      80         [ -  - ]:          0 :       m_problem.errchk( m_system, m_ncomp );
      81                 :        163 :     }
      82                 :            : 
      83                 :            :     //! Find the number of primitive quantities required for this PDE system
      84                 :            :     //! \return The number of primitive quantities required to be stored for
      85                 :            :     //!   this PDE system
      86                 :            :     std::size_t nprim() const
      87                 :            :     {
      88                 :            :       // transport does not need/store any primitive quantities currently
      89                 :            :       return 0;
      90                 :            :     }
      91                 :            : 
      92                 :            :     //! Find the number of materials set up for this PDE system
      93                 :            :     //! \return The number of materials set up for this PDE system
      94                 :            :     std::size_t nmat() const
      95                 :            :     {
      96                 :            :       return m_ncomp;
      97                 :            :     }
      98                 :            : 
      99                 :            :     //! Assign number of DOFs per equation in the PDE system
     100                 :            :     //! \param[in,out] numEqDof Array storing number of Dofs for each PDE
     101                 :            :     //!   equation
     102                 :            :     void numEquationDofs(std::vector< std::size_t >& numEqDof) const
     103                 :            :     {
     104                 :            :       // all equation-dofs initialized to ndofs
     105 [ -  - ][ +  + ]:       1282 :       for (std::size_t i=0; i<m_ncomp; ++i) {
         [ +  + ][ +  + ]
                 [ -  - ]
     106                 :        641 :         numEqDof.push_back(g_inputdeck.get< tag::discr, tag::ndof >());
     107                 :            :       }
     108                 :            :     }
     109                 :            : 
     110                 :            :     //! Determine elements that lie inside the user-defined IC box
     111                 :            :     void IcBoxElems( const tk::Fields&,
     112                 :            :       std::size_t,
     113                 :            :       std::vector< std::unordered_set< std::size_t > >& ) const
     114                 :            :     {}
     115                 :            : 
     116                 :            :     //! Initalize the transport equations for DG
     117                 :            :     //! \param[in] L Element mass matrix
     118                 :            :     //! \param[in] inpoel Element-node connectivity
     119                 :            :     //! \param[in] coord Array of nodal coordinates
     120                 :            :     //! \param[in,out] unk Array of unknowns
     121                 :            :     //! \param[in] t Physical time
     122                 :            :     //! \param[in] nielem Number of internal elements
     123                 :            :     void
     124         [ +  - ]:        765 :     initialize(
     125                 :            :       const tk::Fields& L,
     126                 :            :       const std::vector< std::size_t >& inpoel,
     127                 :            :       const tk::UnsMesh::Coords& coord,
     128                 :            :       const std::vector< std::unordered_set< std::size_t > >& /*inbox*/,
     129                 :            :       tk::Fields& unk,
     130                 :            :       tk::real t,
     131                 :            :       const std::size_t nielem ) const
     132                 :            :     {
     133         [ +  - ]:        765 :       tk::initialize( m_system, m_ncomp, m_offset, L, inpoel, coord,
     134                 :            :                       Problem::initialize, unk, t, nielem );
     135                 :        765 :     }
     136                 :            : 
     137                 :            :     //! Compute the left hand side mass matrix
     138                 :            :     //! \param[in] geoElem Element geometry array
     139                 :            :     //! \param[in,out] l Block diagonal mass matrix
     140                 :            :     void lhs( const tk::Fields& geoElem, tk::Fields& l ) const {
     141                 :        873 :       const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
     142                 :        873 :       tk::mass( m_ncomp, m_offset, ndof, geoElem, l );
     143                 :            :     }
     144                 :            : 
     145                 :            :     //! Update the interface cells to first order dofs
     146                 :            :     //! \details This function resets the high-order terms in interface cells,
     147                 :            :     //!   and is currently not used in transport.
     148                 :            :     void updateInterfaceCells( tk::Fields&,
     149                 :            :       std::size_t,
     150                 :            :       std::vector< std::size_t >& ) const {}
     151                 :            : 
     152                 :            :     //! Update the primitives for this PDE system
     153                 :            :     //! \details This function computes and stores the dofs for primitive
     154                 :            :     //!   quantities, which are currently unused for transport.
     155                 :            :     void updatePrimitives( const tk::Fields&,
     156                 :            :                            const tk::Fields&,
     157                 :            :                            const tk::Fields&,
     158                 :            :                            tk::Fields&,
     159                 :            :                            std::size_t ) const {}
     160                 :            : 
     161                 :            :     //! Clean up the state of trace materials for this PDE system
     162                 :            :     //! \details This function cleans up the state of materials present in trace
     163                 :            :     //!   quantities in each cell. This is currently unused for transport.
     164                 :            :     void cleanTraceMaterial( const tk::Fields&,
     165                 :            :                              tk::Fields&,
     166                 :            :                              tk::Fields&,
     167                 :            :                              std::size_t ) const {}
     168                 :            : 
     169                 :            :     //! Reconstruct second-order solution from first-order
     170                 :            :     //! \param[in] t Physical time
     171                 :            :     //! \param[in] geoFace Face geometry array
     172                 :            :     //! \param[in] geoElem Element geometry array
     173                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     174                 :            :     //! \param[in] esup Elements-surrounding-nodes connectivity
     175                 :            :     //! \param[in] inpoel Element-node connectivity
     176                 :            :     //! \param[in] coord Array of nodal coordinates
     177                 :            :     //! \param[in,out] U Solution vector at recent time step
     178                 :            :     //! \param[in,out] P Primitive vector at recent time step
     179                 :      25875 :     void reconstruct( tk::real t,
     180                 :            :                       const tk::Fields& geoFace,
     181                 :            :                       const tk::Fields& geoElem,
     182                 :            :                       const inciter::FaceData& fd,
     183                 :            :                       const std::map< std::size_t, std::vector< std::size_t > >&
     184                 :            :                         esup,
     185                 :            :                       const std::vector< std::size_t >& inpoel,
     186                 :            :                       const tk::UnsMesh::Coords& coord,
     187                 :            :                       tk::Fields& U,
     188                 :            :                       tk::Fields& P ) const
     189                 :            :     {
     190                 :      25875 :       const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
     191                 :            : 
     192                 :            :       // do reconstruction only if P0P1
     193 [ +  + ][ +  + ]:      25875 :       if (rdof == 4 && g_inputdeck.get< tag::discr, tag::ndof >() == 1) {
     194                 :       6750 :         const auto nelem = fd.Esuel().size()/4;
     195                 :       6750 :         const auto intsharp = g_inputdeck.get< tag::param, tag::transport,
     196                 :       6750 :           tag::intsharp >()[m_system];
     197                 :            : 
     198                 :            :         Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
     199                 :            :                 "vector must equal "+ std::to_string(rdof*m_ncomp) );
     200                 :            :         Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
     201                 :            :                 "Mismatch in inpofa size" );
     202                 :            : 
     203                 :            :         // allocate and initialize matrix and vector for reconstruction
     204                 :            :         std::vector< std::array< std::array< tk::real, 3 >, 3 > >
     205                 :       6750 :           lhs_ls( nelem, {{ {{0.0, 0.0, 0.0}},
     206                 :            :                             {{0.0, 0.0, 0.0}},
     207                 :            :                             {{0.0, 0.0, 0.0}} }} );
     208                 :            :         // specify how many variables need to be reconstructed
     209                 :       6750 :         std::array< std::size_t, 2 > varRange {{0, m_ncomp-1}};
     210                 :            : 
     211                 :            :         std::vector< std::vector< std::array< tk::real, 3 > > >
     212 [ +  - ][ +  - ]:      13500 :           rhs_ls( nelem, std::vector< std::array< tk::real, 3 > >
                 [ -  - ]
     213                 :            :             ( m_ncomp,
     214                 :            :               {{ 0.0, 0.0, 0.0 }} ) );
     215                 :            : 
     216                 :            :         // reconstruct x,y,z-derivatives of unknowns
     217                 :            :         // 0. get lhs matrix, which is only geometry dependent
     218         [ +  - ]:       6750 :         tk::lhsLeastSq_P0P1(fd, geoElem, geoFace, lhs_ls);
     219                 :            : 
     220                 :            :         // 1. internal face contributions
     221         [ +  - ]:       6750 :         tk::intLeastSq_P0P1( m_offset, rdof, fd, geoElem, U, rhs_ls, varRange );
     222                 :            : 
     223                 :            :         // 2. boundary face contributions
     224         [ +  + ]:      47250 :         for (const auto& b : m_bc)
     225                 :      40500 :           tk::bndLeastSqConservedVar_P0P1( m_system, m_ncomp, m_offset, rdof,
     226         [ +  - ]:      40500 :             b.first, fd, geoFace, geoElem, t, b.second, P, U, rhs_ls, varRange );
     227                 :            : 
     228                 :            :         // 3. solve 3x3 least-squares system
     229         [ +  - ]:       6750 :         tk::solveLeastSq_P0P1( m_offset, rdof, lhs_ls, rhs_ls, U, varRange );
     230                 :            : 
     231         [ +  + ]:    1646100 :         for (std::size_t e=0; e<nelem; ++e)
     232                 :            :         {
     233         [ +  - ]:    1639350 :           std::vector< std::size_t > matInt(m_ncomp, 0);
     234 [ +  - ][ -  - ]:    1639350 :           std::vector< tk::real > alAvg(m_ncomp, 0.0);
     235         [ +  + ]:    3278700 :           for (std::size_t k=0; k<m_ncomp; ++k)
     236                 :    1639350 :             alAvg[k] = U(e, k*rdof, m_offset);
     237         [ +  - ]:    1639350 :           auto intInd = interfaceIndicator(m_ncomp, alAvg, matInt);
     238         [ -  + ]:    1639350 :           if ((intsharp > 0) && intInd)
     239                 :            :           {
     240                 :            :             // Reconstruct second-order dofs of volume-fractions in Taylor space
     241                 :            :             // using nodal-stencils, for a good interface-normal estimate
     242         [ -  - ]:          0 :             tk::recoLeastSqExtStencil( rdof, m_offset, e, esup, inpoel, geoElem,
     243                 :            :               U, varRange );
     244                 :            :           }
     245                 :            :         }
     246                 :            : 
     247                 :            :         // 4. transform reconstructed derivatives to Dubiner dofs
     248         [ +  - ]:       6750 :         tk::transform_P0P1( m_offset, rdof, nelem, inpoel, coord, U, varRange );
     249                 :            :       }
     250                 :      25875 :     }
     251                 :            : 
     252                 :            :     //! Limit second-order solution
     253                 :            :     //! \param[in] t Physical time
     254                 :            :     //! \param[in] geoFace Face geometry array
     255                 :            :     //! \param[in] geoElem Element geometry array
     256                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     257                 :            :     //! \param[in] esup Elements surrounding points
     258                 :            :     //! \param[in] inpoel Element-node connectivity
     259                 :            :     //! \param[in] coord Array of nodal coordinates
     260                 :            :     //! \param[in] ndofel Vector of local number of degrees of freedome
     261                 :            : //    //! \param[in] gid Local->global node id map
     262                 :            : //    //! \param[in] bid Local chare-boundary node ids (value) associated to
     263                 :            : //    //!   global node ids (key)
     264                 :            : //    //! \param[in] uNodalExtrm Chare-boundary nodal extrema for conservative
     265                 :            : //    //!   variables
     266                 :            :     //! \param[in,out] U Solution vector at recent time step
     267                 :      25875 :     void limit( [[maybe_unused]] tk::real t,
     268                 :            :                 [[maybe_unused]] const tk::Fields& geoFace,
     269                 :            :                 const tk::Fields& geoElem,
     270                 :            :                 const inciter::FaceData& fd,
     271                 :            :                 const std::map< std::size_t, std::vector< std::size_t > >& esup,
     272                 :            :                 const std::vector< std::size_t >& inpoel,
     273                 :            :                 const tk::UnsMesh::Coords& coord,
     274                 :            :                 const std::vector< std::size_t >& ndofel,
     275                 :            :                 const std::vector< std::size_t >&,
     276                 :            :                 const std::unordered_map< std::size_t, std::size_t >&,
     277                 :            :                 const std::vector< std::vector<tk::real> >&,
     278                 :            :                 const std::vector< std::vector<tk::real> >&,
     279                 :            :                 tk::Fields& U,
     280                 :            :                 tk::Fields&,
     281                 :            :                 std::vector< std::size_t >& ) const
     282                 :            :     {
     283                 :      25875 :       const auto limiter = g_inputdeck.get< tag::discr, tag::limiter >();
     284                 :            : 
     285         [ +  + ]:      25875 :       if (limiter == ctr::LimiterType::WENOP1)
     286                 :       6000 :         WENO_P1( fd.Esuel(), m_offset, U );
     287         [ +  + ]:      19875 :       else if (limiter == ctr::LimiterType::SUPERBEEP1)
     288                 :      19500 :         Superbee_P1( fd.Esuel(), inpoel, ndofel, m_offset, coord, U );
     289         [ -  + ]:        375 :       else if (limiter == ctr::LimiterType::VERTEXBASEDP1)
     290                 :          0 :         VertexBasedTransport_P1( esup, inpoel, ndofel, fd.Esuel().size()/4,
     291                 :          0 :           m_system, m_offset, geoElem, coord, U );
     292                 :      25875 :     }
     293                 :            : 
     294                 :            :     //! Compute right hand side
     295                 :            :     //! \param[in] t Physical time
     296                 :            :     //! \param[in] geoFace Face geometry array
     297                 :            :     //! \param[in] geoElem Element geometry array
     298                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     299                 :            :     //! \param[in] inpoel Element-node connectivity
     300                 :            :     //! \param[in] coord Array of nodal coordinates
     301                 :            :     //! \param[in] U Solution vector at recent time step
     302                 :            :     //! \param[in] P Primitive vector at recent time step
     303                 :            :     //! \param[in] ndofel Vector of local number of degrees of freedom
     304                 :            :     //! \param[in,out] R Right-hand side vector computed
     305                 :      53010 :     void rhs( tk::real t,
     306                 :            :               const tk::Fields& geoFace,
     307                 :            :               const tk::Fields& geoElem,
     308                 :            :               const inciter::FaceData& fd,
     309                 :            :               const std::vector< std::size_t >& inpoel,
     310                 :            :               const std::vector< std::unordered_set< std::size_t > >&,
     311                 :            :               const tk::UnsMesh::Coords& coord,
     312                 :            :               const tk::Fields& U,
     313                 :            :               const tk::Fields& P,
     314                 :            :               const std::vector< std::size_t >& ndofel,
     315                 :            :               tk::Fields& R ) const
     316                 :            :     {
     317                 :      53010 :       const auto ndof = g_inputdeck.get< tag::discr, tag::ndof >();
     318                 :      53010 :       const auto rdof = g_inputdeck.get< tag::discr, tag::rdof >();
     319                 :      53010 :       const auto intsharp = g_inputdeck.get< tag::param, tag::transport,
     320                 :      53010 :         tag::intsharp >()[m_system];
     321                 :            : 
     322                 :            :       Assert( U.nunk() == P.nunk(), "Number of unknowns in solution "
     323                 :            :               "vector and primitive vector at recent time step incorrect" );
     324                 :            :       Assert( U.nunk() == R.nunk(), "Number of unknowns in solution "
     325                 :            :               "vector and right-hand side at recent time step incorrect" );
     326                 :            :       Assert( U.nprop() == rdof*m_ncomp, "Number of components in solution "
     327                 :            :               "vector must equal "+ std::to_string(rdof*m_ncomp) );
     328                 :            :       Assert( P.nprop() == 0, "Number of components in primitive "
     329                 :            :               "vector must equal "+ std::to_string(0) );
     330                 :            :       Assert( R.nprop() == ndof*m_ncomp, "Number of components in right-hand "
     331                 :            :               "side vector must equal "+ std::to_string(ndof*m_ncomp) );
     332                 :            :       Assert( fd.Inpofa().size()/3 == fd.Esuf().size()/2,
     333                 :            :               "Mismatch in inpofa size" );
     334                 :            : 
     335                 :            :       // set rhs to zero
     336                 :            :       R.fill(0.0);
     337                 :            : 
     338                 :            :       // empty vector for non-conservative terms. This vector is unused for
     339                 :            :       // linear transport since, there are no non-conservative terms in the
     340                 :            :       // system of PDEs.
     341                 :      53010 :       std::vector< std::vector < tk::real > > riemannDeriv;
     342                 :            : 
     343                 :      53010 :       std::vector< std::vector< tk::real > > vriem;
     344                 :      53010 :       std::vector< std::vector< tk::real > > riemannLoc;
     345                 :            : 
     346                 :            :       // compute internal surface flux integrals
     347 [ +  - ][ +  - ]:     106020 :       tk::surfInt( m_system, m_ncomp, m_offset, t, ndof, rdof, inpoel, coord,
                 [ -  - ]
     348                 :            :                    fd, geoFace, geoElem, Upwind::flux,
     349                 :            :                    Problem::prescribedVelocity, U, P, ndofel, R, vriem,
     350                 :            :                    riemannLoc, riemannDeriv, intsharp );
     351                 :            : 
     352         [ +  + ]:      53010 :       if(ndof > 1)
     353                 :            :         // compute volume integrals
     354 [ +  - ][ +  - ]:      57375 :         tk::volInt( m_system, m_ncomp, m_offset, t, ndof, rdof,
                 [ -  - ]
     355         [ +  - ]:      19125 :                     fd.Esuel().size()/4, inpoel, coord, geoElem, flux,
     356                 :            :                     Problem::prescribedVelocity, U, P, ndofel, R, intsharp );
     357                 :            : 
     358                 :            :       // compute boundary surface flux integrals
     359         [ +  + ]:     371070 :       for (const auto& b : m_bc)
     360 [ +  - ][ +  - ]:     954180 :         tk::bndSurfInt( m_system, m_ncomp, m_offset, ndof, rdof, b.first, fd,
                 [ -  - ]
     361                 :            :           geoFace, geoElem, inpoel, coord, t, Upwind::flux,
     362         [ +  - ]:     318060 :           Problem::prescribedVelocity, b.second, U, P, ndofel, R, vriem,
     363                 :            :           riemannLoc, riemannDeriv, intsharp );
     364                 :      53010 :     }
     365                 :            : 
     366                 :            :     //! Evaluate the adaptive indicator and mark the ndof for each element
     367                 :            :     //! \param[in] nunk Number of unknowns
     368                 :            :     //! \param[in] coord Array of nodal coordinates
     369                 :            :     //! \param[in] inpoel Element-node connectivity
     370                 :            :     //! \param[in] fd Face connectivity and boundary conditions object
     371                 :            :     //! \param[in] unk Array of unknowns
     372                 :            :     //! \param[in] indicator p-refinement indicator type
     373                 :            :     //! \param[in] ndof Number of degrees of freedom in the solution
     374                 :            :     //! \param[in] ndofmax Max number of degrees of freedom for p-refinement
     375                 :            :     //! \param[in] tolref Tolerance for p-refinement
     376                 :            :     //! \param[in,out] ndofel Vector of local number of degrees of freedome
     377         [ -  - ]:          0 :     void eval_ndof( std::size_t nunk,
     378                 :            :                     [[maybe_unused]] const tk::UnsMesh::Coords& coord,
     379                 :            :                     [[maybe_unused]] const std::vector< std::size_t >& inpoel,
     380                 :            :                     const inciter::FaceData& fd,
     381                 :            :                     const tk::Fields& unk,
     382                 :            :                     inciter::ctr::PrefIndicatorType indicator,
     383                 :            :                     std::size_t ndof,
     384                 :            :                     std::size_t ndofmax,
     385                 :            :                     tk::real tolref,
     386                 :            :                     std::vector< std::size_t >& ndofel ) const
     387                 :            :     {
     388                 :            :       const auto& esuel = fd.Esuel();
     389                 :            : 
     390         [ -  - ]:          0 :       if(indicator == inciter::ctr::PrefIndicatorType::SPECTRAL_DECAY)
     391                 :          0 :         spectral_decay( 1, nunk, esuel, unk, ndof, ndofmax, tolref, ndofel );
     392                 :            :       else
     393 [ -  - ][ -  - ]:          0 :         Throw( "No such adaptive indicator type" );
         [ -  - ][ -  - ]
         [ -  - ][ -  - ]
     394                 :          0 :     }
     395                 :            : 
     396                 :            :     //! Compute the minimum time step size
     397                 :            : //     //! \param[in] U Solution vector at recent time step
     398                 :            : //     //! \param[in] coord Mesh node coordinates
     399                 :            : //     //! \param[in] inpoel Mesh element connectivity
     400                 :            :     //! \return Minimum time step size
     401                 :            :     tk::real dt( const std::array< std::vector< tk::real >, 3 >& /*coord*/,
     402                 :            :                  const std::vector< std::size_t >& /*inpoel*/,
     403                 :            :                  const inciter::FaceData& /*fd*/,
     404                 :            :                  const tk::Fields& /*geoFace*/,
     405                 :            :                  const tk::Fields& /*geoElem*/,
     406                 :            :                  const std::vector< std::size_t >& /*ndofel*/,
     407                 :            :                  const tk::Fields& /*U*/,
     408                 :            :                  const tk::Fields&,
     409                 :            :                  const std::size_t /*nielem*/ ) const
     410                 :            :     {
     411                 :            :       tk::real mindt = std::numeric_limits< tk::real >::max();
     412                 :            :       return mindt;
     413                 :            :     }
     414                 :            : 
     415                 :            :     //! Return analytic field names to be output to file
     416                 :            :     //! \return Vector of strings labelling analytic fields output in file
     417                 :       1646 :     std::vector< std::string > analyticFieldNames() const {
     418                 :            :       std::vector< std::string > n;
     419                 :       1646 :       auto depvar = g_inputdeck.get< tag::param, eq, tag::depvar >()[m_system];
     420         [ +  + ]:       3292 :       for (ncomp_t c=0; c<m_ncomp; ++c)
     421 [ +  - ][ +  - ]:       3292 :         n.push_back( depvar + std::to_string(c) + "_analytic" );
         [ -  + ][ -  + ]
         [ -  - ][ -  - ]
     422                 :       1646 :       return n;
     423                 :            :     }
     424                 :            : 
     425                 :            :     //! Return surface field output going to file
     426                 :            :     std::vector< std::vector< tk::real > >
     427                 :            :     surfOutput( const std::map< int, std::vector< std::size_t > >&,
     428                 :            :                 tk::Fields& ) const
     429                 :            :     {
     430                 :            :       std::vector< std::vector< tk::real > > s; // punt for now
     431                 :            :       return s;
     432                 :            :     }
     433                 :            : 
     434                 :            :     //! Return time history field names to be output to file
     435                 :            :     //! \return Vector of strings labelling time history fields output in file
     436                 :            :     std::vector< std::string > histNames() const {
     437                 :            :       std::vector< std::string > s; // punt for now
     438                 :            :       return s;
     439                 :            :     }
     440                 :            : 
     441                 :            :     //! Return names of integral variables to be output to diagnostics file
     442                 :            :     //! \return Vector of strings labelling integral variables output
     443         [ -  + ]:         49 :     std::vector< std::string > names() const {
     444                 :            :       std::vector< std::string > n;
     445                 :            :       const auto& depvar =
     446         [ -  + ]:         49 :       g_inputdeck.get< tag::param, eq, tag::depvar >().at(m_system);
     447                 :            :       // construct the name of the numerical solution for all components
     448         [ +  + ]:         98 :       for (ncomp_t c=0; c<m_ncomp; ++c)
     449 [ +  - ][ -  + ]:         98 :         n.push_back( depvar + std::to_string(c) );
                 [ -  - ]
     450                 :         49 :       return n;
     451                 :            :     }
     452                 :            : 
     453                 :            :     //! Return analytic solution (if defined by Problem) at xi, yi, zi, t
     454                 :            :     //! \param[in] xi X-coordinate at which to evaluate the analytic solution
     455                 :            :     //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
     456                 :            :     //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
     457                 :            :     //! \param[in] t Physical time at which to evaluate the analytic solution
     458                 :            :     //! \return Vector of analytic solution at given spatial location and time
     459                 :            :     std::vector< tk::real >
     460                 :            :     analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
     461                 :            :     { return Problem::analyticSolution( m_system, m_ncomp, xi, yi, zi, t ); }
     462                 :            : 
     463                 :            :     //! Return analytic solution for conserved variables
     464                 :            :     //! \param[in] xi X-coordinate at which to evaluate the analytic solution
     465                 :            :     //! \param[in] yi Y-coordinate at which to evaluate the analytic solution
     466                 :            :     //! \param[in] zi Z-coordinate at which to evaluate the analytic solution
     467                 :            :     //! \param[in] t Physical time at which to evaluate the analytic solution
     468                 :            :     //! \return Vector of analytic solution at given location and time
     469                 :            :     std::vector< tk::real >
     470                 :            :     solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
     471                 :    7772173 :     { return Problem::initialize( m_system, m_ncomp, xi, yi, zi, t ); }
     472                 :            : 
     473                 :            :     //! Return time history field output evaluated at time history points
     474                 :            :     //! \param[in] h History point data
     475                 :            :     std::vector< std::vector< tk::real > >
     476                 :            :     histOutput( const std::vector< HistData >& h,
     477                 :            :                 const std::vector< std::size_t >&,
     478                 :            :                 const tk::UnsMesh::Coords&,
     479                 :            :                 const tk::Fields&,
     480                 :            :                 const tk::Fields& ) const
     481                 :            :     {
     482                 :          0 :       std::vector< std::vector< tk::real > > Up(h.size()); //punt for now
     483                 :            :       return Up;
     484                 :            :     }
     485                 :            : 
     486                 :            :   private:
     487                 :            :     const Physics m_physics;            //!< Physics policy
     488                 :            :     const Problem m_problem;            //!< Problem policy
     489                 :            :     const ncomp_t m_system;             //!< Equation system index
     490                 :            :     const ncomp_t m_ncomp;              //!< Number of components in this PDE
     491                 :            :     const ncomp_t m_offset;             //!< Offset this PDE operates from
     492                 :            :     //! BC configuration
     493                 :            :     BCStateFn m_bc;
     494                 :            : 
     495                 :            :     //! Evaluate physical flux function for this PDE system
     496                 :            :     //! \param[in] ncomp Number of scalar components in this PDE system
     497                 :            :     //! \param[in] ugp Numerical solution at the Gauss point at which to
     498                 :            :     //!   evaluate the flux
     499                 :            :     //! \param[in] v Prescribed velocity evaluated at the Gauss point at which
     500                 :            :     //!   to evaluate the flux
     501                 :            :     //! \return Flux vectors for all components in this PDE system
     502                 :            :     //! \note The function signature must follow tk::FluxFn
     503                 :            :     static tk::FluxFn::result_type
     504                 :   19672200 :     flux( ncomp_t,
     505                 :            :           ncomp_t ncomp,
     506                 :            :           const std::vector< tk::real >& ugp,
     507                 :            :           const std::vector< std::array< tk::real, 3 > >& v )
     508                 :            :     {
     509                 :            :       Assert( ugp.size() == ncomp, "Size mismatch" );
     510                 :            :       Assert( v.size() == ncomp, "Size mismatch" );
     511                 :            : 
     512                 :   19672200 :       std::vector< std::array< tk::real, 3 > > fl( ugp.size() );
     513                 :            : 
     514         [ +  + ]:   39344400 :       for (ncomp_t c=0; c<ncomp; ++c)
     515                 :   19672200 :         fl[c] = {{ v[c][0] * ugp[c], v[c][1] * ugp[c], v[c][2] * ugp[c] }};
     516                 :            : 
     517                 :   19672200 :       return fl;
     518                 :            :     }
     519                 :            : 
     520                 :            :     //! \brief Boundary state function providing the left and right state of a
     521                 :            :     //!   face at extrapolation boundaries
     522                 :            :     //! \param[in] ul Left (domain-internal) state
     523                 :            :     //! \return Left and right states for all scalar components in this PDE
     524                 :            :     //!   system
     525                 :            :     //! \note The function signature must follow tk::StateFn
     526                 :            :     static tk::StateFn::result_type
     527                 :   10720020 :     extrapolate( ncomp_t, ncomp_t, const std::vector< tk::real >& ul,
     528                 :            :                  tk::real, tk::real, tk::real, tk::real,
     529                 :            :                  const std::array< tk::real, 3 >& )
     530                 :            :     {
     531                 :   10720020 :       return {{ ul, ul }};
     532                 :            :     }
     533                 :            : 
     534                 :            :     //! \brief Boundary state function providing the left and right state of a
     535                 :            :     //!   face at extrapolation boundaries
     536                 :            :     //! \param[in] ul Left (domain-internal) state
     537                 :            :     //! \return Left and right states for all scalar components in this PDE
     538                 :            :     //!   system
     539                 :            :     //! \note The function signature must follow tk::StateFn
     540                 :            :     static tk::StateFn::result_type
     541                 :     319200 :     inlet( ncomp_t, ncomp_t, const std::vector< tk::real >& ul,
     542                 :            :            tk::real, tk::real, tk::real, tk::real,
     543                 :            :            const std::array< tk::real, 3 >& )
     544                 :            :     {
     545                 :     319200 :       auto ur = ul;
     546                 :            :       std::fill( begin(ur), end(ur), 0.0 );
     547         [ +  - ]:     319200 :       return {{ ul, std::move(ur) }};
     548                 :            :     }
     549                 :            : 
     550                 :            :     //! \brief Boundary state function providing the left and right state of a
     551                 :            :     //!   face at outlet boundaries
     552                 :            :     //! \param[in] ul Left (domain-internal) state
     553                 :            :     //! \return Left and right states for all scalar components in this PDE
     554                 :            :     //!   system
     555                 :            :     //! \note The function signature must follow tk::StateFn
     556                 :            :     static tk::StateFn::result_type
     557                 :     873600 :     outlet( ncomp_t, ncomp_t, const std::vector< tk::real >& ul,
     558                 :            :             tk::real, tk::real, tk::real, tk::real,
     559                 :            :             const std::array< tk::real, 3 >& )
     560                 :            :     {
     561                 :     873600 :       return {{ ul, ul }};
     562                 :            :     }
     563                 :            : 
     564                 :            :     //! \brief Boundary state function providing the left and right state of a
     565                 :            :     //!   face at Dirichlet boundaries
     566                 :            :     //! \param[in] system Equation system index
     567                 :            :     //! \param[in] ncomp Number of scalar components in this PDE system
     568                 :            :     //! \param[in] ul Left (domain-internal) state
     569                 :            :     //! \param[in] x X-coordinate at which to compute the states
     570                 :            :     //! \param[in] y Y-coordinate at which to compute the states
     571                 :            :     //! \param[in] z Z-coordinate at which to compute the states
     572                 :            :     //! \param[in] t Physical time
     573                 :            :     //! \return Left and right states for all scalar components in this PDE
     574                 :            :     //!   system
     575                 :            :     //! \note The function signature must follow tk::StateFn
     576                 :            :     static tk::StateFn::result_type
     577                 :    4135860 :     dirichlet( ncomp_t system, ncomp_t ncomp, const std::vector< tk::real >& ul,
     578                 :            :                tk::real x, tk::real y, tk::real z, tk::real t,
     579                 :            :                const std::array< tk::real, 3 >& )
     580                 :            :     {
     581                 :    4135860 :       return {{ ul, Problem::initialize( system, ncomp, x, y, z, t ) }};
     582                 :            :     }
     583                 :            : };
     584                 :            : 
     585                 :            : } // dg::
     586                 :            : } // inciter::
     587                 :            : 
     588                 :            : #endif // DGTransport_h

Generated by: LCOV version 1.14