Quinoa all test code coverage report
Current view: top level - Inciter - DG.cpp (source / functions) Hit Total Coverage
Commit: -128-NOTFOUND Lines: 822 1219 67.4 %
Date: 2025-10-21 17:21:13 Functions: 38 44 86.4 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 813 2048 39.7 %

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/Inciter/DG.cpp
       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     DG advances a system of PDEs with the discontinuous Galerkin scheme
       9                 :            :   \details   DG advances a system of partial differential equations (PDEs) using
      10                 :            :     discontinuous Galerkin (DG) finite element (FE) spatial discretization (on
      11                 :            :     tetrahedron elements) combined with Runge-Kutta (RK) time stepping.
      12                 :            :   \see The documentation in DG.h.
      13                 :            : */
      14                 :            : // *****************************************************************************
      15                 :            : 
      16                 :            : #include <algorithm>
      17                 :            : #include <numeric>
      18                 :            : #include <sstream>
      19                 :            : 
      20                 :            : #include "DG.hpp"
      21                 :            : #include "Discretization.hpp"
      22                 :            : #include "DGPDE.hpp"
      23                 :            : #include "DiagReducer.hpp"
      24                 :            : #include "DerivedData.hpp"
      25                 :            : #include "ElemDiagnostics.hpp"
      26                 :            : #include "Inciter/InputDeck/InputDeck.hpp"
      27                 :            : #include "Refiner.hpp"
      28                 :            : #include "Limiter.hpp"
      29                 :            : #include "Reorder.hpp"
      30                 :            : #include "Vector.hpp"
      31                 :            : #include "Around.hpp"
      32                 :            : #include "Integrate/Basis.hpp"
      33                 :            : #include "FieldOutput.hpp"
      34                 :            : #include "ChareStateCollector.hpp"
      35                 :            : #include "PDE/MultiMat/MultiMatIndexing.hpp"
      36                 :            : #include "Integrate/Quadrature.hpp"
      37                 :            : 
      38                 :            : #include <fstream>
      39                 :            : 
      40                 :            : // ignore old-style-casts required for lapack/blas calls
      41                 :            : #if defined(__clang__)
      42                 :            :   #pragma clang diagnostic ignored "-Wold-style-cast"
      43                 :            : #endif
      44                 :            : 
      45                 :            : // Lapacke forward declarations
      46                 :            : extern "C" {
      47                 :            : 
      48                 :            : using lapack_int = long;
      49                 :            : 
      50                 :            : #define LAPACK_ROW_MAJOR 101
      51                 :            : 
      52                 :            : extern lapack_int LAPACKE_dgesv( int, lapack_int, lapack_int, double*,
      53                 :            :   lapack_int, lapack_int*, double*, lapack_int );
      54                 :            : 
      55                 :            : }
      56                 :            : 
      57                 :            : namespace inciter {
      58                 :            : 
      59                 :            : extern ctr::InputDeck g_inputdeck;
      60                 :            : extern std::vector< DGPDE > g_dgpde;
      61                 :            : 
      62                 :            : //! Runge-Kutta coefficients
      63                 :            : static const std::array< std::array< tk::real, 3 >, 2 >
      64                 :            :   rkcoef{{ {{ 0.0, 3.0/4.0, 1.0/3.0 }}, {{ 1.0, 1.0/4.0, 2.0/3.0 }} }};
      65                 :            : 
      66                 :            : //! Implicit-Explicit Runge-Kutta Coefficients
      67                 :            : static const tk::real rk_gamma = (2.0-std::sqrt(2.0))/2.0;
      68                 :            : static const tk::real rk_delta = -2.0*std::sqrt(2.0)/3.0;
      69                 :            : static const tk::real c2 =
      70                 :            :   (27.0 + std::pow(2187.0-1458.0*std::sqrt(2.0),1.0/3.0)
      71                 :            :    + 9.0*std::pow(3.0+2.0*std::sqrt(2.0),1.0/3.0))/54.0;
      72                 :            : static const tk::real c3 = c2/(6.0*std::pow(c2,2.0)-3.0*c2+1.0);
      73                 :            : static const tk::real b2 = (3.0*c2-1.0)/(6.0*std::pow(c2,2.0));
      74                 :            : static const tk::real b3 =
      75                 :            :   (6.0*std::pow(c2,2.0)-3.0*c2+1.0)/(6.0*std::pow(c2,2.0));
      76                 :            : static const tk::real a22_impl = c2;
      77                 :            : static const tk::real a21_expl = c2;
      78                 :            : static const tk::real a32_expl = c3;
      79                 :            : static const tk::real a33_impl =
      80                 :            :   (1.0/6.0-b2*std::pow(c2,2.0)-b3*c2*c3)/(b3*(c3-c2));
      81                 :            : static const tk::real a32_impl = a33_impl-c3;
      82                 :            : static const std::array< std::array< tk::real, 3 >, 2 >
      83                 :            :   expl_rkcoef{{ {{ 0.0, 0.0, b2 }},
      84                 :            :                 {{ a21_expl, a32_expl, b3 }} }};
      85                 :            : static const std::array< std::array< tk::real, 3 >, 2>
      86                 :            :   impl_rkcoef{{ {{ 0.0, a32_impl, b2 }},
      87                 :            :                 {{ a22_impl, a33_impl, b3}} }};
      88                 :            : static const std::array< tk::real, 10 > mass_dubiner( tk::massMatrixDubiner() );
      89                 :            : 
      90                 :            : } // inciter::
      91                 :            : 
      92                 :            : extern tk::CProxy_ChareStateCollector stateProxy;
      93                 :            : 
      94                 :            : using inciter::DG;
      95                 :            : 
      96                 :        671 : DG::DG( const CProxy_Discretization& disc,
      97                 :            :         const CProxy_Ghosts& ghostsproxy,
      98                 :            :         const std::map< int, std::vector< std::size_t > >& bface,
      99                 :            :         const std::map< int, std::vector< std::size_t > >& /* bnode */,
     100                 :        671 :         const std::vector< std::size_t >& triinpoel ) :
     101                 :            :   m_disc( disc ),
     102                 :            :   m_ghosts( ghostsproxy ),
     103                 :            :   m_ndof_NodalExtrm( 3 ), // for the first order derivatives in 3 directions
     104                 :            :   m_nsol( 0 ),
     105                 :            :   m_ninitsol( 0 ),
     106                 :            :   m_nlim( 0 ),
     107                 :            :   m_nnod( 0 ),
     108                 :            :   m_nrefine( 0 ),
     109                 :            :   m_nsmooth( 0 ),
     110                 :            :   m_nreco( 0 ),
     111                 :            :   m_nnodalExtrema( 0 ),
     112         [ +  - ]:        671 :   m_nstiffeq( g_dgpde[Disc()->MeshId()].nstiffeq() ),
     113         [ +  - ]:        671 :   m_nnonstiffeq( g_dgpde[Disc()->MeshId()].nnonstiffeq() ),
     114                 :        671 :   m_u( Disc()->Inpoel().size()/4,
     115                 :       1342 :        g_inputdeck.get< tag::rdof >()*
     116                 :        671 :        g_inputdeck.get< tag::ncomp >() ),
     117                 :            :   m_un( m_u.nunk(), m_u.nprop() ),
     118                 :       1342 :   m_p( m_u.nunk(), g_inputdeck.get< tag::rdof >()*
     119 [ +  - ][ +  - ]:        671 :     g_dgpde[Disc()->MeshId()].nprim() ),
     120                 :            :   m_rhs( m_u.nunk(),
     121                 :       1342 :          g_inputdeck.get< tag::ndof >()*
     122                 :        671 :          g_inputdeck.get< tag::ncomp >() ),
     123                 :            :   m_rhsprev( m_u.nunk(), m_rhs.nprop() ),
     124                 :       1342 :   m_stiffrhs( m_u.nunk(), g_inputdeck.get< tag::ndof >()*
     125 [ +  - ][ +  - ]:        671 :               g_dgpde[Disc()->MeshId()].nstiffeq() ),
     126                 :       1342 :   m_stiffrhsprev( m_u.nunk(), g_inputdeck.get< tag::ndof >()*
     127 [ +  - ][ +  - ]:        671 :                   g_dgpde[Disc()->MeshId()].nstiffeq() ),
     128         [ +  - ]:        671 :   m_stiffEqIdx( g_dgpde[Disc()->MeshId()].nstiffeq() ),
     129         [ +  - ]:        671 :   m_nonStiffEqIdx( g_dgpde[Disc()->MeshId()].nnonstiffeq() ),
     130                 :            :   m_mtInv(
     131                 :        671 :     tk::invMassMatTaylorRefEl(g_inputdeck.get< tag::rdof >()) ),
     132                 :            :   m_uNodalExtrm(),
     133                 :            :   m_pNodalExtrm(),
     134                 :            :   m_uNodalExtrmc(),
     135                 :            :   m_pNodalExtrmc(),
     136         [ +  - ]:        671 :   m_npoin( Disc()->Coord()[0].size() ),
     137                 :            :   m_diag(),
     138                 :            :   m_nstage( 3 ),
     139                 :            :   m_stage( 0 ),
     140                 :            :   m_ndof(),
     141                 :            :   m_interface(),
     142                 :            :   m_numEqDof(),
     143                 :            :   m_uc(),
     144                 :            :   m_pc(),
     145                 :            :   m_ndofc(),
     146                 :            :   m_interfacec(),
     147                 :            :   m_initial( 1 ),
     148                 :            :   m_uElemfields( m_u.nunk(),
     149                 :        671 :                  g_inputdeck.get< tag::ncomp >() ),
     150                 :            :   m_pElemfields( m_u.nunk(),
     151                 :        671 :                  m_p.nprop() / g_inputdeck.get< tag::rdof >() ),
     152                 :            :   m_uNodefields( m_npoin,
     153                 :        671 :                  g_inputdeck.get< tag::ncomp >() ),
     154                 :            :   m_pNodefields( m_npoin,
     155                 :        671 :                  m_p.nprop() / g_inputdeck.get< tag::rdof >() ),
     156                 :            :   m_uNodefieldsc(),
     157                 :            :   m_pNodefieldsc(),
     158                 :            :   m_outmesh(),
     159                 :            :   m_boxelems(),
     160 [ +  - ][ +  - ]:       9394 :   m_shockmarker(m_u.nunk(), 1)
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
         [ +  - ][ +  - ]
                 [ +  - ]
     161                 :            : // *****************************************************************************
     162                 :            : //  Constructor
     163                 :            : //! \param[in] disc Discretization proxy
     164                 :            : //! \param[in] bface Boundary-faces mapped to side set ids
     165                 :            : //! \param[in] triinpoel Boundary-face connectivity
     166                 :            : // *****************************************************************************
     167                 :            : {
     168 [ +  + ][ +  + ]:       1304 :   if (g_inputdeck.get< tag::cmd, tag::chare >() ||
     169         [ +  + ]:        633 :       g_inputdeck.get< tag::cmd, tag::quiescence >())
     170 [ +  - ][ +  - ]:        446 :     stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
         [ +  - ][ +  - ]
                 [ +  - ]
     171                 :            :                                         "DG" );
     172                 :            : 
     173                 :            :   // assign number of dofs for each equation in all pde systems
     174 [ +  - ][ +  - ]:        671 :   g_dgpde[Disc()->MeshId()].numEquationDofs(m_numEqDof);
     175                 :            : 
     176                 :            :   // Allocate storage for the vector of nodal extrema
     177 [ +  - ][ +  - ]:        671 :   m_uNodalExtrm.resize( Disc()->Bid().size(),
     178                 :       1342 :     std::vector<tk::real>( 2 * m_ndof_NodalExtrm *
     179         [ +  - ]:        671 :     g_inputdeck.get< tag::ncomp >() ) );
     180 [ +  - ][ +  - ]:        671 :   m_pNodalExtrm.resize( Disc()->Bid().size(),
     181                 :       1342 :     std::vector<tk::real>( 2 * m_ndof_NodalExtrm *
     182         [ +  - ]:        671 :     m_p.nprop() / g_inputdeck.get< tag::rdof >() ) );
     183                 :            : 
     184                 :            :   // Initialization for the buffer vector of nodal extrema
     185         [ +  - ]:        671 :   resizeNodalExtremac();
     186                 :            : 
     187                 :        671 :   usesAtSync = true;    // enable migration at AtSync
     188                 :            : 
     189                 :        671 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
     190                 :            : 
     191                 :            :   // Enable SDAG wait for initially building the solution vector and limiting
     192         [ +  - ]:        671 :   if (m_initial) {
     193 [ +  - ][ +  - ]:        671 :     thisProxy[ thisIndex ].wait4sol();
     194 [ +  + ][ +  - ]:        671 :     if (pref) thisProxy[ thisIndex ].wait4refine();
                 [ +  - ]
     195 [ +  - ][ +  - ]:        671 :     thisProxy[ thisIndex ].wait4smooth();
     196 [ +  - ][ +  - ]:        671 :     thisProxy[ thisIndex ].wait4lim();
     197 [ +  - ][ +  - ]:        671 :     thisProxy[ thisIndex ].wait4nod();
     198 [ +  - ][ +  - ]:        671 :     thisProxy[ thisIndex ].wait4reco();
     199 [ +  - ][ +  - ]:        671 :     thisProxy[ thisIndex ].wait4nodalExtrema();
     200                 :            :   }
     201                 :            : 
     202 [ +  - ][ +  - ]:       1342 :   m_ghosts[thisIndex].insert(m_disc, bface, triinpoel, m_u.nunk(),
     203 [ +  - ][ +  - ]:       1342 :     CkCallback(CkIndex_DG::resizeSolVectors(), thisProxy[thisIndex]));
                 [ +  - ]
     204                 :            : 
     205                 :            :   // global-sync to call doneInserting on m_ghosts
     206         [ +  - ]:        671 :   auto meshid = Disc()->MeshId();
     207         [ +  - ]:        671 :   contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
     208 [ +  - ][ +  - ]:       1342 :     CkCallback(CkReductionTarget(Transporter,doneInsertingGhosts),
     209         [ +  - ]:        671 :     Disc()->Tr()) );
     210                 :        671 : }
     211                 :            : 
     212                 :            : void
     213                 :        529 : DG::registerReducers()
     214                 :            : // *****************************************************************************
     215                 :            : //  Configure Charm++ reduction types
     216                 :            : //! \details Since this is a [initnode] routine, the runtime system executes the
     217                 :            : //!   routine exactly once on every logical node early on in the Charm++ init
     218                 :            : //!   sequence. Must be static as it is called without an object. See also:
     219                 :            : //!   Section "Initializations at Program Startup" at in the Charm++ manual
     220                 :            : //!   http://charm.cs.illinois.edu/manuals/html/charm++/manual.html.
     221                 :            : // *****************************************************************************
     222                 :            : {
     223                 :        529 :   ElemDiagnostics::registerReducers();
     224                 :        529 : }
     225                 :            : 
     226                 :            : void
     227                 :      12084 : DG::ResumeFromSync()
     228                 :            : // *****************************************************************************
     229                 :            : //  Return from migration
     230                 :            : //! \details This is called when load balancing (LB) completes. The presence of
     231                 :            : //!   this function does not affect whether or not we block on LB.
     232                 :            : // *****************************************************************************
     233                 :            : {
     234 [ -  + ][ -  - ]:      12084 :   if (Disc()->It() == 0) Throw( "it = 0 in ResumeFromSync()" );
         [ -  - ][ -  - ]
     235                 :            : 
     236         [ +  - ]:      12084 :   if (!g_inputdeck.get< tag::cmd, tag::nonblocking >()) next();
     237                 :      12084 : }
     238                 :            : 
     239                 :            : void
     240                 :        671 : DG::resizeSolVectors()
     241                 :            : // *****************************************************************************
     242                 :            : // Resize solution vectors after extension due to Ghosts and continue with setup
     243                 :            : // *****************************************************************************
     244                 :            : {
     245                 :            :   // Resize solution vectors, lhs and rhs by the number of ghost tets
     246 [ +  - ][ +  - ]:        671 :   m_u.resize( myGhosts()->m_nunk );
     247 [ +  - ][ +  - ]:        671 :   m_un.resize( myGhosts()->m_nunk );
     248 [ +  - ][ +  - ]:        671 :   m_p.resize( myGhosts()->m_nunk );
     249 [ +  - ][ +  - ]:        671 :   m_rhs.resize( myGhosts()->m_nunk );
     250 [ +  - ][ +  - ]:        671 :   m_rhsprev.resize( myGhosts()->m_nunk );
     251 [ +  - ][ +  - ]:        671 :   m_stiffrhs.resize( myGhosts()->m_nunk );
     252 [ +  - ][ +  - ]:        671 :   m_stiffrhsprev.resize( myGhosts()->m_nunk );
     253                 :            : 
     254                 :            :   // Size communication buffer for solution and number of degrees of freedom
     255 [ +  + ][ +  - ]:       2684 :   for (auto& n : m_ndofc) n.resize( myGhosts()->m_bid.size() );
                 [ +  - ]
     256 [ +  + ][ +  - ]:       2684 :   for (auto& u : m_uc) u.resize( myGhosts()->m_bid.size() );
                 [ +  - ]
     257 [ +  + ][ +  - ]:       2684 :   for (auto& p : m_pc) p.resize( myGhosts()->m_bid.size() );
                 [ +  - ]
     258 [ +  + ][ +  - ]:       1342 :   for (auto& i : m_interfacec) i.resize( myGhosts()->m_bid.size() );
                 [ +  - ]
     259                 :            : 
     260                 :            :   // Initialize number of degrees of freedom in mesh elements
     261                 :        671 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
     262         [ +  + ]:        671 :   if( pref )
     263                 :            :   {
     264                 :        134 :     const auto ndofmax = g_inputdeck.get< tag::pref, tag::ndofmax >();
     265 [ +  - ][ +  - ]:        134 :     m_ndof.resize( myGhosts()->m_nunk, ndofmax );
     266                 :            :   }
     267                 :            :   else
     268                 :            :   {
     269                 :        537 :     const auto ndof = g_inputdeck.get< tag::ndof >();
     270 [ +  - ][ +  - ]:        537 :     m_ndof.resize( myGhosts()->m_nunk, ndof );
     271                 :            :   }
     272 [ +  - ][ +  - ]:        671 :   m_interface.resize( myGhosts()->m_nunk, 0 );
     273                 :            : 
     274                 :            :   // Ensure that we also have all the geometry and connectivity data
     275                 :            :   // (including those of ghosts)
     276 [ +  - ][ -  + ]:        671 :   Assert( myGhosts()->m_geoElem.nunk() == m_u.nunk(),
         [ -  - ][ -  - ]
                 [ -  - ]
     277                 :            :     "GeoElem unknowns size mismatch" );
     278                 :            : 
     279                 :            :   // Signal the runtime system that all workers have received their adjacency
     280 [ +  - ][ +  - ]:        671 :   std::vector< std::size_t > meshdata{ myGhosts()->m_initial, Disc()->MeshId() };
                 [ +  - ]
     281         [ +  - ]:        671 :   contribute( meshdata, CkReduction::sum_ulong,
     282 [ +  - ][ +  - ]:       1342 :     CkCallback(CkReductionTarget(Transporter,comfinal), Disc()->Tr()) );
                 [ +  - ]
     283                 :        671 : }
     284                 :            : 
     285                 :            : void
     286                 :        671 : DG::setup()
     287                 :            : // *****************************************************************************
     288                 :            : // Set initial conditions, generate lhs, output mesh
     289                 :            : // *****************************************************************************
     290                 :            : {
     291 [ +  + ][ +  + ]:       1304 :   if (g_inputdeck.get< tag::cmd, tag::chare >() ||
     292         [ +  + ]:        633 :       g_inputdeck.get< tag::cmd, tag::quiescence >())
     293 [ +  - ][ +  - ]:        446 :     stateProxy.ckLocalBranch()->insert( "DG", thisIndex, CkMyPe(), Disc()->It(),
         [ +  - ][ +  - ]
     294                 :            :                                         "setup" );
     295                 :            : 
     296                 :        671 :   auto d = Disc();
     297                 :            : 
     298                 :            :   // Compute left-hand side of discrete PDEs
     299                 :        671 :   lhs();
     300                 :            : 
     301                 :            :   // Determine elements inside user-defined IC box
     302                 :       1342 :   g_dgpde[d->MeshId()].IcBoxElems( myGhosts()->m_geoElem,
     303                 :        671 :     myGhosts()->m_fd.Esuel().size()/4, m_boxelems );
     304                 :            : 
     305                 :            :   // Compute volume of user-defined box IC
     306         [ +  - ]:        671 :   d->boxvol( {}, {}, 0 );      // punt for now
     307                 :            : 
     308                 :            :   // Query time history field output labels from all PDEs integrated
     309                 :        671 :   const auto& hist_points = g_inputdeck.get< tag::history_output, tag::point >();
     310         [ -  + ]:        671 :   if (!hist_points.empty()) {
     311                 :          0 :     std::vector< std::string > histnames;
     312         [ -  - ]:          0 :     auto n = g_dgpde[d->MeshId()].histNames();
     313         [ -  - ]:          0 :     histnames.insert( end(histnames), begin(n), end(n) );
     314         [ -  - ]:          0 :     d->histheader( std::move(histnames) );
     315                 :            :   }
     316                 :            : 
     317                 :            :   // If working with IMEX-RK, Store stiff equations into m_stiffEqIdx
     318         [ -  + ]:        671 :   if (g_inputdeck.get< tag::imex_runge_kutta >())
     319                 :            :   {
     320                 :          0 :     g_dgpde[Disc()->MeshId()].setStiffEqIdx(m_stiffEqIdx);
     321                 :          0 :     g_dgpde[Disc()->MeshId()].setNonStiffEqIdx(m_nonStiffEqIdx);
     322                 :            :   }
     323                 :        671 : }
     324                 :            : 
     325                 :            : void
     326                 :        671 : DG::box( tk::real v, const std::vector< tk::real >& )
     327                 :            : // *****************************************************************************
     328                 :            : // Receive total box IC volume and set conditions in box
     329                 :            : //! \param[in] v Total volume within user-specified box
     330                 :            : // *****************************************************************************
     331                 :            : {
     332                 :        671 :   auto d = Disc();
     333                 :            : 
     334                 :            :   // Store user-defined box IC volume
     335                 :        671 :   d->Boxvol() = v;
     336                 :            : 
     337                 :            :   // Set initial conditions for all PDEs
     338                 :       2013 :   g_dgpde[d->MeshId()].initialize( myGhosts()->m_geoElem, myGhosts()->m_inpoel,
     339                 :        671 :     myGhosts()->m_coord, m_boxelems, d->ElemBlockId(), m_u, d->T(),
     340                 :        671 :     myGhosts()->m_fd.Esuel().size()/4 );
     341                 :       1342 :   g_dgpde[d->MeshId()].updatePrimitives( m_u, myGhosts()->m_geoElem, m_p,
     342                 :        671 :     myGhosts()->m_fd.Esuel().size()/4, m_ndof );
     343                 :            : 
     344                 :        671 :   m_un = m_u;
     345                 :            : 
     346                 :            :   // Output initial conditions to file (regardless of whether it was requested)
     347 [ +  - ][ +  - ]:        671 :   startFieldOutput( CkCallback(CkIndex_DG::start(), thisProxy[thisIndex]) );
                 [ +  - ]
     348                 :        671 : }
     349                 :            : 
     350                 :            : void
     351                 :        671 : DG::start()
     352                 :            : // *****************************************************************************
     353                 :            : //  Start time stepping
     354                 :            : // *****************************************************************************
     355                 :            : {
     356                 :            :   // Free memory storing output mesh
     357                 :        671 :   m_outmesh.destroy();
     358                 :            : 
     359                 :            :   // Start timer measuring time stepping wall clock time
     360                 :        671 :   Disc()->Timer().zero();
     361                 :            :   // Zero grind-timer
     362                 :        671 :   Disc()->grindZero();
     363                 :            :   // Start time stepping by computing the size of the next time step)
     364                 :        671 :   next();
     365                 :        671 : }
     366                 :            : 
     367                 :            : void
     368                 :      13426 : DG::startFieldOutput( CkCallback c )
     369                 :            : // *****************************************************************************
     370                 :            : // Start preparing fields for output to file
     371                 :            : //! \param[in] c Function to continue with after the write
     372                 :            : // *****************************************************************************
     373                 :            : {
     374                 :            :   // No field output in benchmark mode or if field output frequency not hit
     375 [ +  + ][ +  + ]:      13426 :   if (g_inputdeck.get< tag::cmd, tag::benchmark >() || !fieldOutput()) {
                 [ +  + ]
     376                 :            : 
     377                 :      12101 :     c.send();
     378                 :            : 
     379                 :            :   } else {
     380                 :            : 
     381                 :            :     // Optionally refine mesh for field output
     382                 :       1325 :     auto d = Disc();
     383                 :            : 
     384         [ +  + ]:       1325 :     if (refinedOutput()) {
     385                 :            : 
     386 [ +  - ][ +  - ]:         33 :       const auto& tr = tk::remap( myGhosts()->m_fd.Triinpoel(), d->Gid() );
     387 [ +  - ][ +  - ]:         33 :       d->Ref()->outref( myGhosts()->m_fd.Bface(), {}, tr, c );
                 [ +  - ]
     388                 :            : 
     389                 :            :     } else {
     390                 :            : 
     391                 :            :       // cut off ghosts from mesh connectivity and coordinates
     392 [ +  - ][ +  - ]:       1292 :       const auto& tr = tk::remap( myGhosts()->m_fd.Triinpoel(), d->Gid() );
     393         [ +  - ]:       2584 :       extractFieldOutput( {}, d->Chunk(), d->Coord(), {}, {},
     394         [ +  - ]:       1292 :                           d->NodeCommMap(), myGhosts()->m_fd.Bface(), {}, tr, c );
     395                 :            : 
     396                 :            :     }
     397                 :            : 
     398                 :            :   }
     399                 :      13426 : }
     400                 :            : 
     401                 :            : void
     402                 :      38265 : DG::next()
     403                 :            : // *****************************************************************************
     404                 :            : // Advance equations to next time step
     405                 :            : // *****************************************************************************
     406                 :            : {
     407                 :      38265 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
     408                 :            : 
     409                 :      38265 :   auto d = Disc();
     410                 :            : 
     411 [ +  + ][ +  + ]:      38265 :   if (pref && m_stage == 0 && d->T() > 0)
         [ +  + ][ +  + ]
     412                 :       3272 :     g_dgpde[d->MeshId()].eval_ndof( myGhosts()->m_nunk, myGhosts()->m_coord,
     413                 :       1636 :                   myGhosts()->m_inpoel,
     414                 :       1636 :                   myGhosts()->m_fd, m_u, m_p,
     415                 :       1636 :                   g_inputdeck.get< tag::pref, tag::indicator >(),
     416                 :       1636 :                   g_inputdeck.get< tag::ndof >(),
     417                 :       1636 :                   g_inputdeck.get< tag::pref, tag::ndofmax >(),
     418                 :       1636 :                   g_inputdeck.get< tag::pref, tag::tolref >(),
     419                 :       1636 :                   m_ndof );
     420                 :            : 
     421                 :            :   // communicate solution ghost data (if any)
     422         [ +  + ]:      38265 :   if (myGhosts()->m_sendGhost.empty())
     423                 :       3390 :     comsol_complete();
     424                 :            :   else
     425 [ +  - ][ +  + ]:     419925 :     for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
     426         [ +  - ]:     770100 :       std::vector< std::size_t > tetid( ghostdata.size() );
     427         [ +  - ]:     770100 :       std::vector< std::vector< tk::real > > u( ghostdata.size() ),
     428         [ +  - ]:     770100 :                                              prim( ghostdata.size() );
     429         [ +  - ]:     770100 :       std::vector< std::size_t > interface( ghostdata.size() );
     430         [ +  - ]:     385050 :       std::vector< std::size_t > ndof( ghostdata.size() );
     431                 :     385050 :       std::size_t j = 0;
     432         [ +  + ]:    6854730 :       for(const auto& i : ghostdata) {
     433 [ +  - ][ -  + ]:    6469680 :         Assert( i < myGhosts()->m_fd.Esuel().size()/4,
         [ -  - ][ -  - ]
                 [ -  - ]
     434                 :            :           "Sending solution ghost data" );
     435                 :    6469680 :         tetid[j] = i;
     436         [ +  - ]:    6469680 :         u[j] = m_u[i];
     437         [ +  - ]:    6469680 :         prim[j] = m_p[i];
     438 [ +  + ][ +  + ]:    6469680 :         if (pref && m_stage == 0) {
     439                 :     395110 :           ndof[j] = m_ndof[i];
     440                 :     395110 :           interface[j] = m_interface[i];
     441                 :            :         }
     442                 :    6469680 :         ++j;
     443                 :            :       }
     444 [ +  - ][ +  - ]:     385050 :       thisProxy[ cid ].comsol( thisIndex, m_stage, tetid, u, prim, interface, ndof );
     445                 :            :     }
     446                 :            : 
     447                 :      38265 :   ownsol_complete();
     448                 :      38265 : }
     449                 :            : 
     450                 :            : void
     451                 :     385050 : DG::comsol( int fromch,
     452                 :            :             std::size_t fromstage,
     453                 :            :             const std::vector< std::size_t >& tetid,
     454                 :            :             const std::vector< std::vector< tk::real > >& u,
     455                 :            :             const std::vector< std::vector< tk::real > >& prim,
     456                 :            :             const std::vector< std::size_t >& interface,
     457                 :            :             const std::vector< std::size_t >& ndof )
     458                 :            : // *****************************************************************************
     459                 :            : //  Receive chare-boundary solution ghost data from neighboring chares
     460                 :            : //! \param[in] fromch Sender chare id
     461                 :            : //! \param[in] fromstage Sender chare time step stage
     462                 :            : //! \param[in] tetid Ghost tet ids we receive solution data for
     463                 :            : //! \param[in] u Solution ghost data
     464                 :            : //! \param[in] prim Primitive variables in ghost cells
     465                 :            : //! \param[in] interface Interface marker in ghost cells
     466                 :            : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
     467                 :            : //! \details This function receives contributions to the unlimited solution
     468                 :            : //!   from fellow chares.
     469                 :            : // *****************************************************************************
     470                 :            : {
     471 [ -  + ][ -  - ]:     385050 :   Assert( u.size() == tetid.size(), "Size mismatch in DG::comsol()" );
         [ -  - ][ -  - ]
     472 [ -  + ][ -  - ]:     385050 :   Assert( prim.size() == tetid.size(), "Size mismatch in DG::comsol()" );
         [ -  - ][ -  - ]
     473                 :            : 
     474                 :     385050 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
     475                 :            : 
     476 [ +  + ][ +  + ]:     385050 :   if (pref && fromstage == 0) {
     477 [ -  + ][ -  - ]:      13280 :     Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsol()" );
         [ -  - ][ -  - ]
     478 [ -  + ][ -  - ]:      13280 :     Assert( interface.size() == tetid.size(), "Size mismatch in DG::comsol()" );
         [ -  - ][ -  - ]
     479                 :            :   }
     480                 :            : 
     481                 :            :   // Find local-to-ghost tet id map for sender chare
     482                 :     385050 :   const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
     483                 :            : 
     484         [ +  + ]:    6854730 :   for (std::size_t i=0; i<tetid.size(); ++i) {
     485         [ +  - ]:    6469680 :     auto j = tk::cref_find( n, tetid[i] );
     486 [ +  - ][ -  + ]:    6469680 :     Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
         [ -  - ][ -  - ]
                 [ -  - ]
     487                 :            :       "Receiving solution non-ghost data" );
     488 [ +  - ][ +  - ]:    6469680 :     auto b = tk::cref_find( myGhosts()->m_bid, j );
     489 [ -  + ][ -  - ]:    6469680 :     Assert( b < m_uc[0].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
     490 [ -  + ][ -  - ]:    6469680 :     Assert( b < m_pc[0].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
     491         [ +  - ]:    6469680 :     m_uc[0][b] = u[i];
     492         [ +  - ]:    6469680 :     m_pc[0][b] = prim[i];
     493 [ +  + ][ +  + ]:    6469680 :     if (pref && fromstage == 0) {
     494 [ -  + ][ -  - ]:     395110 :       Assert( b < m_ndofc[0].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
     495                 :     395110 :       m_ndofc[0][b] = ndof[i];
     496 [ -  + ][ -  - ]:     395110 :       Assert( b < m_interfacec[0].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
     497                 :     395110 :       m_interfacec[0][b] = interface[i];
     498                 :            :     }
     499                 :            :   }
     500                 :            : 
     501                 :            :   // if we have received all solution ghost contributions from neighboring
     502                 :            :   // chares (chares we communicate along chare-boundary faces with), and
     503                 :            :   // contributed our solution to these neighbors, proceed to reconstructions
     504         [ +  + ]:     385050 :   if (++m_nsol == myGhosts()->m_sendGhost.size()) {
     505                 :      34875 :     m_nsol = 0;
     506                 :      34875 :     comsol_complete();
     507                 :            :   }
     508                 :     385050 : }
     509                 :            : 
     510                 :            : void
     511                 :       1325 : DG::extractFieldOutput(
     512                 :            :   const std::vector< std::size_t >& /*ginpoel*/,
     513                 :            :   const tk::UnsMesh::Chunk& chunk,
     514                 :            :   const tk::UnsMesh::Coords& coord,
     515                 :            :   const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
     516                 :            :   const std::unordered_map< std::size_t, std::size_t >& addedTets,
     517                 :            :   const tk::NodeCommMap& nodeCommMap,
     518                 :            :   const std::map< int, std::vector< std::size_t > >& bface,
     519                 :            :   const std::map< int, std::vector< std::size_t > >& /* bnode */,
     520                 :            :   const std::vector< std::size_t >& triinpoel,
     521                 :            :   CkCallback c )
     522                 :            : // *****************************************************************************
     523                 :            : // Extract field output going to file
     524                 :            : //! \param[in] chunk Field-output mesh chunk (connectivity and global<->local
     525                 :            : //!    id maps)
     526                 :            : //! \param[in] coord Field-output mesh node coordinates
     527                 :            : //! \param[in] addedTets Field-output mesh cells and their parents (local ids)
     528                 :            : //! \param[in] nodeCommMap Field-output mesh node communication map
     529                 :            : //! \param[in] bface Field-output meshndary-faces mapped to side set ids
     530                 :            : //! \param[in] triinpoel Field-output mesh boundary-face connectivity
     531                 :            : //! \param[in] c Function to continue with after the write
     532                 :            : // *****************************************************************************
     533                 :            : {
     534                 :       1325 :   m_outmesh.chunk = chunk;
     535                 :       1325 :   m_outmesh.coord = coord;
     536                 :       1325 :   m_outmesh.triinpoel = triinpoel;
     537                 :       1325 :   m_outmesh.bface = bface;
     538                 :       1325 :   m_outmesh.nodeCommMap = nodeCommMap;
     539                 :            : 
     540                 :       1325 :   const auto& inpoel = std::get< 0 >( chunk );
     541                 :            : 
     542                 :            :   // Evaluate element solution on incoming mesh
     543                 :       1325 :   evalSolution( *Disc(), inpoel, coord, addedTets, m_ndof, m_u, m_p,
     544                 :       1325 :     m_uElemfields, m_pElemfields, m_uNodefields, m_pNodefields );
     545                 :            : 
     546                 :            :   // Send node fields contributions to neighbor chares
     547         [ +  + ]:       1325 :   if (nodeCommMap.empty())
     548                 :        161 :     comnodeout_complete();
     549                 :            :   else {
     550                 :       1164 :     const auto& lid = std::get< 2 >( chunk );
     551         [ +  - ]:       2328 :     auto esup = tk::genEsup( inpoel, 4 );
     552         [ +  + ]:      11544 :     for(const auto& [ch,nodes] : nodeCommMap) {
     553                 :            :       // Pack node field data in chare boundary nodes
     554                 :            :       std::vector< std::vector< tk::real > >
     555 [ +  - ][ +  - ]:      31140 :         lu( m_uNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
     556                 :            :       std::vector< std::vector< tk::real > >
     557 [ +  - ][ +  - ]:      31140 :         lp( m_pNodefields.nprop(), std::vector< tk::real >( nodes.size() ) );
     558         [ +  + ]:      69172 :       for (std::size_t f=0; f<m_uNodefields.nprop(); ++f) {
     559                 :      58792 :         std::size_t j = 0;
     560         [ +  + ]:     603358 :         for (auto g : nodes)
     561 [ +  - ][ +  - ]:     544566 :           lu[f][j++] = m_uNodefields(tk::cref_find(lid,g),f);
     562                 :            :       }
     563         [ +  + ]:      21008 :       for (std::size_t f=0; f<m_pNodefields.nprop(); ++f) {
     564                 :      10628 :         std::size_t j = 0;
     565         [ +  + ]:     120108 :         for (auto g : nodes)
     566 [ +  - ][ +  - ]:     109480 :           lp[f][j++] = m_pNodefields(tk::cref_find(lid,g),f);
     567                 :            :       }
     568                 :            :       // Pack (partial) number of elements surrounding chare boundary nodes
     569         [ +  - ]:      10380 :       std::vector< std::size_t > nesup( nodes.size() );
     570                 :      10380 :       std::size_t j = 0;
     571         [ +  + ]:     107480 :       for (auto g : nodes) {
     572         [ +  - ]:      97100 :         auto i = tk::cref_find( lid, g );
     573                 :      97100 :         nesup[j++] = esup.second[i+1] - esup.second[i];
     574                 :            :       }
     575 [ +  - ][ +  - ]:      31140 :       thisProxy[ch].comnodeout(
     576         [ +  - ]:      20760 :         std::vector<std::size_t>(begin(nodes),end(nodes)), nesup, lu, lp );
     577                 :            :     }
     578                 :            :   }
     579                 :            : 
     580         [ +  - ]:       1325 :   ownnod_complete( c, addedTets );
     581                 :       1325 : }
     582                 :            : 
     583                 :      38265 : void DG::p_refine()
     584                 :            : // *****************************************************************************
     585                 :            : // Add the protective layer for ndof refinement
     586                 :            : // *****************************************************************************
     587                 :            : {
     588                 :      38265 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
     589                 :            : 
     590                 :            :   // Combine own and communicated contributions of unreconstructed solution and
     591                 :            :   // degrees of freedom in cells (if p-adaptive)
     592 [ +  - ][ +  + ]:    6507945 :   for (const auto& b : myGhosts()->m_bid) {
     593 [ -  + ][ -  - ]:    6469680 :     Assert( m_uc[0][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
         [ -  - ][ -  - ]
     594 [ -  + ][ -  - ]:    6469680 :     Assert( m_pc[0][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
         [ -  - ][ -  - ]
     595         [ +  + ]:  173340045 :     for (std::size_t c=0; c<m_u.nprop(); ++c) {
     596         [ +  - ]:  166870365 :       m_u(b.first,c) = m_uc[0][b.second][c];
     597                 :            :     }
     598         [ +  + ]:   34324605 :     for (std::size_t c=0; c<m_p.nprop(); ++c) {
     599         [ +  - ]:   27854925 :       m_p(b.first,c) = m_pc[0][b.second][c];
     600                 :            :     }
     601 [ +  + ][ +  + ]:    6469680 :     if (pref && m_stage == 0) {
     602                 :     395110 :       m_ndof[ b.first ] = m_ndofc[0][ b.second ];
     603                 :     395110 :       m_interface[ b.first ] = m_interfacec[0][ b.second ];
     604                 :            :     }
     605                 :            :   }
     606                 :            : 
     607 [ +  + ][ +  + ]:      38265 :   if (pref && m_stage==0) refine_ndof();
     608                 :            : 
     609         [ +  + ]:      38265 :   if (!pref) {
     610                 :            :     // if p-refinement is not configured, proceed directly to reconstructions
     611                 :      32955 :     reco();
     612                 :            :   }
     613                 :            :   else {
     614                 :            :     // if p-refinement is configured, do refine-smoothing before reconstruction
     615         [ +  + ]:       5310 :     if (myGhosts()->m_sendGhost.empty())
     616                 :        150 :       comrefine_complete();
     617                 :            :     else
     618 [ +  - ][ +  + ]:      45000 :       for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
     619         [ +  - ]:      79680 :         std::vector< std::size_t > tetid( ghostdata.size() );
     620         [ +  - ]:      79680 :         std::vector< std::vector< tk::real > > u( ghostdata.size() ),
     621         [ +  - ]:      79680 :                                                prim( ghostdata.size() );
     622         [ +  - ]:      39840 :         std::vector< std::size_t > ndof( ghostdata.size() );
     623                 :      39840 :         std::size_t j = 0;
     624         [ +  + ]:    1225170 :         for(const auto& i : ghostdata) {
     625 [ +  - ][ -  + ]:    1185330 :           Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending refined ndof  "
         [ -  - ][ -  - ]
                 [ -  - ]
     626                 :            :             "data" );
     627                 :    1185330 :           tetid[j] = i;
     628 [ +  - ][ +  + ]:    1185330 :           if (pref && m_stage == 0) ndof[j] = m_ndof[i];
     629                 :    1185330 :           ++j;
     630                 :            :         }
     631 [ +  - ][ +  - ]:      39840 :         thisProxy[ cid ].comrefine( thisIndex, tetid, ndof );
     632                 :            :       }
     633                 :            : 
     634                 :       5310 :     ownrefine_complete();
     635                 :            :   }
     636                 :      38265 : }
     637                 :            : 
     638                 :            : void
     639                 :      39840 : DG::comrefine( int fromch,
     640                 :            :                const std::vector< std::size_t >& tetid,
     641                 :            :                const std::vector< std::size_t >& ndof )
     642                 :            : // *****************************************************************************
     643                 :            : //  Receive chare-boundary ghost data from neighboring chares
     644                 :            : //! \param[in] fromch Sender chare id
     645                 :            : //! \param[in] tetid Ghost tet ids we receive solution data for
     646                 :            : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
     647                 :            : //! \details This function receives contributions to the refined ndof data
     648                 :            : //!   from fellow chares.
     649                 :            : // *****************************************************************************
     650                 :            : {
     651                 :      39840 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
     652                 :            : 
     653 [ +  - ][ +  + ]:      39840 :   if (pref && m_stage == 0)
     654 [ -  + ][ -  - ]:      13280 :     Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comrefine()" );
         [ -  - ][ -  - ]
     655                 :            : 
     656                 :            :   // Find local-to-ghost tet id map for sender chare
     657                 :      39840 :   const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
     658                 :            : 
     659         [ +  + ]:    1225170 :   for (std::size_t i=0; i<tetid.size(); ++i) {
     660         [ +  - ]:    1185330 :     auto j = tk::cref_find( n, tetid[i] );
     661 [ +  - ][ -  + ]:    1185330 :     Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
         [ -  - ][ -  - ]
                 [ -  - ]
     662                 :            :       "Receiving solution non-ghost data" );
     663 [ +  - ][ +  - ]:    1185330 :     auto b = tk::cref_find( myGhosts()->m_bid, j );
     664 [ +  - ][ +  + ]:    1185330 :     if (pref && m_stage == 0) {
     665 [ -  + ][ -  - ]:     395110 :       Assert( b < m_ndofc[1].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
     666                 :     395110 :       m_ndofc[1][b] = ndof[i];
     667                 :            :     }
     668                 :            :   }
     669                 :            : 
     670                 :            :   // if we have received all solution ghost contributions from neighboring
     671                 :            :   // chares (chares we communicate along chare-boundary faces with), and
     672                 :            :   // contributed our solution to these neighbors, proceed to limiting
     673         [ +  + ]:      39840 :   if (++m_nrefine == myGhosts()->m_sendGhost.size()) {
     674                 :       5160 :     m_nrefine = 0;
     675                 :       5160 :     comrefine_complete();
     676                 :            :   }
     677                 :      39840 : }
     678                 :            : 
     679                 :            : void
     680                 :       5310 : DG::smooth()
     681                 :            : // *****************************************************************************
     682                 :            : // Smooth the refined ndof distribution
     683                 :            : // *****************************************************************************
     684                 :            : {
     685                 :       5310 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
     686                 :            : 
     687 [ +  - ][ +  + ]:    1190640 :   for (const auto& b : myGhosts()->m_bid) {
     688 [ +  - ][ +  + ]:    1185330 :     if (pref && m_stage == 0)
     689                 :     395110 :       m_ndof[ b.first ] = m_ndofc[1][ b.second ];
     690                 :            :   }
     691                 :            : 
     692 [ +  - ][ +  + ]:       5310 :   if (pref && m_stage==0) smooth_ndof();
     693                 :            : 
     694         [ +  + ]:       5310 :   if (myGhosts()->m_sendGhost.empty())
     695                 :        150 :     comsmooth_complete();
     696                 :            :   else
     697 [ +  - ][ +  + ]:      45000 :     for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
     698         [ +  - ]:      79680 :       std::vector< std::size_t > tetid( ghostdata.size() );
     699                 :      39840 :       std::vector< std::size_t > ndof;
     700                 :      39840 :       std::size_t j = 0;
     701         [ +  + ]:    1225170 :       for(const auto& i : ghostdata) {
     702 [ +  - ][ -  + ]:    1185330 :         Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending ndof data" );
         [ -  - ][ -  - ]
                 [ -  - ]
     703                 :    1185330 :         tetid[j] = i;
     704 [ +  - ][ +  + ]:    1185330 :         if (pref && m_stage == 0) ndof.push_back( m_ndof[i] );
                 [ +  - ]
     705                 :    1185330 :         ++j;
     706                 :            :       }
     707 [ +  - ][ +  - ]:      39840 :       thisProxy[ cid ].comsmooth( thisIndex, tetid, ndof );
     708                 :            :     }
     709                 :            : 
     710                 :       5310 :   ownsmooth_complete();
     711                 :       5310 : }
     712                 :            : 
     713                 :            : void
     714                 :      39840 : DG::comsmooth( int fromch,
     715                 :            :                const std::vector< std::size_t >& tetid,
     716                 :            :                const std::vector< std::size_t >& ndof )
     717                 :            : // *****************************************************************************
     718                 :            : //  Receive chare-boundary ghost data from neighboring chares
     719                 :            : //! \param[in] fromch Sender chare id
     720                 :            : //! \param[in] tetid Ghost tet ids we receive solution data for
     721                 :            : //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
     722                 :            : //! \details This function receives contributions to the smoothed ndof data
     723                 :            : //!   from fellow chares.
     724                 :            : // *****************************************************************************
     725                 :            : {
     726                 :      39840 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
     727                 :            : 
     728 [ +  - ][ +  + ]:      39840 :   if (pref && m_stage == 0)
     729 [ -  + ][ -  - ]:      13280 :     Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsmooth()" );
         [ -  - ][ -  - ]
     730                 :            : 
     731                 :      39840 :   const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
     732                 :            : 
     733         [ +  + ]:    1225170 :   for (std::size_t i=0; i<tetid.size(); ++i) {
     734         [ +  - ]:    1185330 :     auto j = tk::cref_find( n, tetid[i] );
     735 [ +  - ][ -  + ]:    1185330 :     Assert( j >= myGhosts()->m_fd.Esuel().size()/4, "Receiving ndof data" );
         [ -  - ][ -  - ]
                 [ -  - ]
     736 [ +  - ][ +  - ]:    1185330 :     auto b = tk::cref_find( myGhosts()->m_bid, j );
     737 [ +  - ][ +  + ]:    1185330 :     if (pref && m_stage == 0) {
     738 [ -  + ][ -  - ]:     395110 :       Assert( b < m_ndofc[2].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
     739                 :     395110 :       m_ndofc[2][b] = ndof[i];
     740                 :            :     }
     741                 :            :   }
     742                 :            : 
     743         [ +  + ]:      39840 :   if (++m_nsmooth == myGhosts()->m_sendGhost.size()) {
     744                 :       5160 :     m_nsmooth = 0;
     745                 :       5160 :     comsmooth_complete();
     746                 :            :   }
     747                 :      39840 : }
     748                 :            : 
     749                 :            : void
     750                 :      38265 : DG::reco()
     751                 :            : // *****************************************************************************
     752                 :            : // Compute reconstructions
     753                 :            : // *****************************************************************************
     754                 :            : {
     755                 :      38265 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
     756                 :      38265 :   const auto rdof = g_inputdeck.get< tag::rdof >();
     757                 :            : 
     758                 :            :   // Combine own and communicated contributions of unreconstructed solution and
     759                 :            :   // degrees of freedom in cells (if p-adaptive)
     760         [ +  + ]:      38265 :   if (pref) {
     761 [ +  - ][ +  + ]:    1190640 :     for (const auto& b : myGhosts()->m_bid) {
     762         [ +  + ]:    1185330 :       if (m_stage == 0) {
     763                 :     395110 :         m_ndof[ b.first ] = m_ndofc[2][ b.second ];
     764                 :            :       }
     765                 :            :     }
     766                 :            :   }
     767                 :            : 
     768                 :      38265 :   auto d = Disc();
     769 [ +  + ][ +  + ]:      38265 :   if (pref && m_stage==0) {
     770                 :       1770 :     g_dgpde[d->MeshId()].resetAdapSol( myGhosts()->m_fd, m_u, m_p, m_ndof );
     771                 :            :   }
     772                 :            : 
     773         [ +  + ]:      38265 :   if (rdof > 1)
     774                 :            :     // Reconstruct second-order solution and primitive quantities
     775                 :      49680 :     g_dgpde[d->MeshId()].reconstruct( d->T(), myGhosts()->m_geoFace,
     776                 :      24840 :       myGhosts()->m_geoElem,
     777                 :      24840 :       myGhosts()->m_fd, myGhosts()->m_esup, myGhosts()->m_inpoel,
     778                 :      24840 :       myGhosts()->m_coord, m_u, m_p, pref, m_ndof );
     779                 :            : 
     780                 :            :   // Send reconstructed solution to neighboring chares
     781         [ +  + ]:      38265 :   if (myGhosts()->m_sendGhost.empty())
     782                 :       3390 :     comreco_complete();
     783                 :            :   else
     784 [ +  - ][ +  + ]:     419925 :     for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
     785         [ +  - ]:     770100 :       std::vector< std::size_t > tetid( ghostdata.size() );
     786         [ +  - ]:     770100 :       std::vector< std::vector< tk::real > > u( ghostdata.size() ),
     787         [ +  - ]:     385050 :                                              prim( ghostdata.size() );
     788                 :     385050 :       std::size_t j = 0;
     789         [ +  + ]:    6854730 :       for(const auto& i : ghostdata) {
     790 [ +  - ][ -  + ]:    6469680 :         Assert( i < myGhosts()->m_fd.Esuel().size()/4, "Sending reconstructed ghost "
         [ -  - ][ -  - ]
                 [ -  - ]
     791                 :            :           "data" );
     792                 :    6469680 :         tetid[j] = i;
     793         [ +  - ]:    6469680 :         u[j] = m_u[i];
     794         [ +  - ]:    6469680 :         prim[j] = m_p[i];
     795                 :    6469680 :         ++j;
     796                 :            :       }
     797 [ +  - ][ +  - ]:     385050 :       thisProxy[ cid ].comreco( thisIndex, tetid, u, prim );
     798                 :            :     }
     799                 :            : 
     800                 :      38265 :   ownreco_complete();
     801                 :      38265 : }
     802                 :            : 
     803                 :            : void
     804                 :     385050 : DG::comreco( int fromch,
     805                 :            :              const std::vector< std::size_t >& tetid,
     806                 :            :              const std::vector< std::vector< tk::real > >& u,
     807                 :            :              const std::vector< std::vector< tk::real > >& prim )
     808                 :            : // *****************************************************************************
     809                 :            : //  Receive chare-boundary reconstructed ghost data from neighboring chares
     810                 :            : //! \param[in] fromch Sender chare id
     811                 :            : //! \param[in] tetid Ghost tet ids we receive solution data for
     812                 :            : //! \param[in] u Reconstructed high-order solution
     813                 :            : //! \param[in] prim Limited high-order primitive quantities
     814                 :            : //! \details This function receives contributions to the reconstructed solution
     815                 :            : //!   from fellow chares.
     816                 :            : // *****************************************************************************
     817                 :            : {
     818 [ -  + ][ -  - ]:     385050 :   Assert( u.size() == tetid.size(), "Size mismatch in DG::comreco()" );
         [ -  - ][ -  - ]
     819 [ -  + ][ -  - ]:     385050 :   Assert( prim.size() == tetid.size(), "Size mismatch in DG::comreco()" );
         [ -  - ][ -  - ]
     820                 :            : 
     821                 :            :   // Find local-to-ghost tet id map for sender chare
     822                 :     385050 :   const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
     823                 :            : 
     824         [ +  + ]:    6854730 :   for (std::size_t i=0; i<tetid.size(); ++i) {
     825         [ +  - ]:    6469680 :     auto j = tk::cref_find( n, tetid[i] );
     826 [ +  - ][ -  + ]:    6469680 :     Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
         [ -  - ][ -  - ]
                 [ -  - ]
     827                 :            :       "Receiving solution non-ghost data" );
     828 [ +  - ][ +  - ]:    6469680 :     auto b = tk::cref_find( myGhosts()->m_bid, j );
     829 [ -  + ][ -  - ]:    6469680 :     Assert( b < m_uc[1].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
     830 [ -  + ][ -  - ]:    6469680 :     Assert( b < m_pc[1].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
     831         [ +  - ]:    6469680 :     m_uc[1][b] = u[i];
     832         [ +  - ]:    6469680 :     m_pc[1][b] = prim[i];
     833                 :            :   }
     834                 :            : 
     835                 :            :   // if we have received all solution ghost contributions from neighboring
     836                 :            :   // chares (chares we communicate along chare-boundary faces with), and
     837                 :            :   // contributed our solution to these neighbors, proceed to limiting
     838         [ +  + ]:     385050 :   if (++m_nreco == myGhosts()->m_sendGhost.size()) {
     839                 :      34875 :     m_nreco = 0;
     840                 :      34875 :     comreco_complete();
     841                 :            :   }
     842                 :     385050 : }
     843                 :            : 
     844                 :            : void
     845                 :      38265 : DG::nodalExtrema()
     846                 :            : // *****************************************************************************
     847                 :            : // Compute nodal extrema at chare-boundary nodes. Extrema at internal nodes
     848                 :            : // are calculated in limiter function.
     849                 :            : // *****************************************************************************
     850                 :            : {
     851         [ +  - ]:      38265 :   auto d = Disc();
     852         [ +  - ]:      76530 :   auto gid = d->Gid();
     853         [ +  - ]:      76530 :   auto bid = d->Bid();
     854                 :      38265 :   const auto rdof = g_inputdeck.get< tag::rdof >();
     855                 :      38265 :   const auto ncomp = m_u.nprop() / rdof;
     856                 :      38265 :   const auto nprim = m_p.nprop() / rdof;
     857                 :            : 
     858                 :            :   // Combine own and communicated contributions of unlimited solution, and
     859                 :            :   // if a p-adaptive algorithm is used, degrees of freedom in cells
     860 [ +  - ][ +  + ]:    6507945 :   for (const auto& [boundary, localtet] : myGhosts()->m_bid) {
     861 [ -  + ][ -  - ]:    6469680 :     Assert( m_uc[1][localtet].size() == m_u.nprop(), "ncomp size mismatch" );
         [ -  - ][ -  - ]
     862 [ -  + ][ -  - ]:    6469680 :     Assert( m_pc[1][localtet].size() == m_p.nprop(), "ncomp size mismatch" );
         [ -  - ][ -  - ]
     863         [ +  + ]:  173340045 :     for (std::size_t c=0; c<m_u.nprop(); ++c) {
     864         [ +  - ]:  166870365 :       m_u(boundary,c) = m_uc[1][localtet][c];
     865                 :            :     }
     866         [ +  + ]:   34324605 :     for (std::size_t c=0; c<m_p.nprop(); ++c) {
     867         [ +  - ]:   27854925 :       m_p(boundary,c) = m_pc[1][localtet][c];
     868                 :            :     }
     869                 :            :   }
     870                 :            : 
     871                 :            :   // Initialize nodal extrema vector
     872                 :      38265 :   auto large = std::numeric_limits< tk::real >::max();
     873         [ +  + ]:    1317765 :   for(std::size_t i = 0; i<bid.size(); i++)
     874                 :            :   {
     875         [ +  + ]:    9304560 :     for (std::size_t c=0; c<ncomp; ++c)
     876                 :            :     {
     877         [ +  + ]:   32100240 :       for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
     878                 :            :       {
     879                 :   24075180 :         auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
     880                 :   24075180 :         auto min_mark = max_mark + 1;
     881                 :   24075180 :         m_uNodalExtrm[i][max_mark] = -large;
     882                 :   24075180 :         m_uNodalExtrm[i][min_mark] =  large;
     883                 :            :       }
     884                 :            :     }
     885         [ +  + ]:    3460350 :     for (std::size_t c=0; c<nprim; ++c)
     886                 :            :     {
     887         [ +  + ]:    8723400 :       for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
     888                 :            :       {
     889                 :    6542550 :         auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
     890                 :    6542550 :         auto min_mark = max_mark + 1;
     891                 :    6542550 :         m_pNodalExtrm[i][max_mark] = -large;
     892                 :    6542550 :         m_pNodalExtrm[i][min_mark] =  large;
     893                 :            :       }
     894                 :            :     }
     895                 :            :   }
     896                 :            : 
     897                 :            :   // Evaluate the max/min value for the chare-boundary nodes
     898         [ +  + ]:      38265 :   if(rdof > 4) {
     899         [ +  - ]:       7110 :       evalNodalExtrmRefEl(ncomp, nprim, m_ndof_NodalExtrm, d->bndel(),
     900 [ +  - ][ +  - ]:       7110 :         myGhosts()->m_inpoel, gid, bid, m_u, m_p, m_uNodalExtrm, m_pNodalExtrm);
     901                 :            :   }
     902                 :            : 
     903                 :            :   // Communicate extrema at nodes to other chares on chare-boundary
     904         [ +  + ]:      38265 :   if (d->NodeCommMap().empty())        // in serial we are done
     905         [ +  - ]:       3390 :     comnodalExtrema_complete();
     906                 :            :   else  // send nodal extrema to chare-boundary nodes to fellow chares
     907                 :            :   {
     908         [ +  + ]:     419925 :     for (const auto& [c,n] : d->NodeCommMap()) {
     909 [ +  - ][ +  - ]:    1155150 :       std::vector< std::vector< tk::real > > g1( n.size() ), g2( n.size() );
     910                 :     385050 :       std::size_t j = 0;
     911         [ +  + ]:    3269040 :       for (auto i : n)
     912                 :            :       {
     913         [ +  - ]:    2883990 :         auto p = tk::cref_find(d->Bid(),i);
     914         [ +  - ]:    2883990 :         g1[ j   ] = m_uNodalExtrm[ p ];
     915         [ +  - ]:    2883990 :         g2[ j++ ] = m_pNodalExtrm[ p ];
     916                 :            :       }
     917 [ +  - ][ +  - ]:     385050 :       thisProxy[c].comnodalExtrema( std::vector<std::size_t>(begin(n),end(n)),
                 [ +  - ]
     918                 :            :         g1, g2 );
     919                 :            :     }
     920                 :            :   }
     921         [ +  - ]:      38265 :   ownnodalExtrema_complete();
     922                 :      38265 : }
     923                 :            : 
     924                 :            : void
     925                 :     385050 : DG::comnodalExtrema( const std::vector< std::size_t >& gid,
     926                 :            :                      const std::vector< std::vector< tk::real > >& G1,
     927                 :            :                      const std::vector< std::vector< tk::real > >& G2 )
     928                 :            : // *****************************************************************************
     929                 :            : //  Receive contributions to nodal extrema on chare-boundaries
     930                 :            : //! \param[in] gid Global mesh node IDs at which we receive grad contributions
     931                 :            : //! \param[in] G1 Partial contributions of extrema for conservative variables to
     932                 :            : //!   chare-boundary nodes
     933                 :            : //! \param[in] G2 Partial contributions of extrema for primitive variables to
     934                 :            : //!   chare-boundary nodes
     935                 :            : //! \details This function receives contributions to m_uNodalExtrm/m_pNodalExtrm
     936                 :            : //!   , which stores nodal extrems at mesh chare-boundary nodes. While
     937                 :            : //!   m_uNodalExtrm/m_pNodalExtrm stores own contributions, m_uNodalExtrmc
     938                 :            : //!   /m_pNodalExtrmc collects the neighbor chare contributions during
     939                 :            : //!   communication.
     940                 :            : // *****************************************************************************
     941                 :            : {
     942 [ -  + ][ -  - ]:     385050 :   Assert( G1.size() == gid.size(), "Size mismatch" );
         [ -  - ][ -  - ]
     943 [ -  + ][ -  - ]:     385050 :   Assert( G2.size() == gid.size(), "Size mismatch" );
         [ -  - ][ -  - ]
     944                 :            : 
     945                 :     385050 :   const auto rdof = g_inputdeck.get< tag::rdof >();
     946                 :     385050 :   const auto ncomp = m_u.nprop() / rdof;
     947                 :     385050 :   const auto nprim = m_p.nprop() / rdof;
     948                 :            : 
     949         [ +  + ]:    3269040 :   for (std::size_t i=0; i<gid.size(); ++i)
     950                 :            :   {
     951                 :    2883990 :     auto& u = m_uNodalExtrmc[gid[i]];
     952                 :    2883990 :     auto& p = m_pNodalExtrmc[gid[i]];
     953         [ +  + ]:   20524560 :     for (std::size_t c=0; c<ncomp; ++c)
     954                 :            :     {
     955         [ +  + ]:   70562280 :       for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
     956                 :            :       {
     957                 :   52921710 :         auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
     958                 :   52921710 :         auto min_mark = max_mark + 1;
     959                 :   52921710 :         u[max_mark] = std::max( G1[i][max_mark], u[max_mark] );
     960                 :   52921710 :         u[min_mark] = std::min( G1[i][min_mark], u[min_mark] );
     961                 :            :       }
     962                 :            :     }
     963         [ +  + ]:    7291440 :     for (std::size_t c=0; c<nprim; ++c)
     964                 :            :     {
     965         [ +  + ]:   17629800 :       for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
     966                 :            :       {
     967                 :   13222350 :         auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
     968                 :   13222350 :         auto min_mark = max_mark + 1;
     969                 :   13222350 :         p[max_mark] = std::max( G2[i][max_mark], p[max_mark] );
     970                 :   13222350 :         p[min_mark] = std::min( G2[i][min_mark], p[min_mark] );
     971                 :            :       }
     972                 :            :     }
     973                 :            :   }
     974                 :            : 
     975         [ +  + ]:     385050 :   if (++m_nnodalExtrema == Disc()->NodeCommMap().size())
     976                 :            :   {
     977                 :      34875 :     m_nnodalExtrema = 0;
     978                 :      34875 :     comnodalExtrema_complete();
     979                 :            :   }
     980                 :     385050 : }
     981                 :            : 
     982                 :      38936 : void DG::resizeNodalExtremac()
     983                 :            : // *****************************************************************************
     984                 :            : //  Resize the buffer vector of nodal extrema
     985                 :            : // *****************************************************************************
     986                 :            : {
     987                 :      38936 :   const auto rdof = g_inputdeck.get< tag::rdof >();
     988                 :      38936 :   const auto ncomp = m_u.nprop() / rdof;
     989                 :      38936 :   const auto nprim = m_p.nprop() / rdof;
     990                 :            : 
     991                 :      38936 :   auto large = std::numeric_limits< tk::real >::max();
     992 [ +  - ][ +  + ]:     431010 :   for (const auto& [c,n] : Disc()->NodeCommMap())
     993                 :            :   {
     994         [ +  + ]:    3322548 :     for (auto i : n) {
     995         [ +  - ]:    2930474 :       auto& u = m_uNodalExtrmc[i];
     996         [ +  - ]:    2930474 :       auto& p = m_pNodalExtrmc[i];
     997         [ +  - ]:    2930474 :       u.resize( 2*m_ndof_NodalExtrm*ncomp, large );
     998         [ +  - ]:    2930474 :       p.resize( 2*m_ndof_NodalExtrm*nprim, large );
     999                 :            : 
    1000                 :            :       // Initialize the minimum nodal extrema
    1001         [ +  + ]:   11721896 :       for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
    1002                 :            :       {
    1003         [ +  + ]:   62506722 :         for(std::size_t k = 0; k < ncomp; k++)
    1004                 :   53715300 :           u[2*k*m_ndof_NodalExtrm+2*idof] = -large;
    1005         [ +  + ]:   22160034 :         for(std::size_t k = 0; k < nprim; k++)
    1006                 :   13368612 :           p[2*k*m_ndof_NodalExtrm+2*idof] = -large;
    1007                 :            :       }
    1008                 :            :     }
    1009                 :            :   }
    1010                 :      38936 : }
    1011                 :            : 
    1012                 :       7110 : void DG::evalNodalExtrmRefEl(
    1013                 :            :   const std::size_t ncomp,
    1014                 :            :   const std::size_t nprim,
    1015                 :            :   const std::size_t ndof_NodalExtrm,
    1016                 :            :   const std::vector< std::size_t >& bndel,
    1017                 :            :   const std::vector< std::size_t >& inpoel,
    1018                 :            :   const std::vector< std::size_t >& gid,
    1019                 :            :   const std::unordered_map< std::size_t, std::size_t >& bid,
    1020                 :            :   const tk::Fields& U,
    1021                 :            :   const tk::Fields& P,
    1022                 :            :   std::vector< std::vector<tk::real> >& uNodalExtrm,
    1023                 :            :   std::vector< std::vector<tk::real> >& pNodalExtrm )
    1024                 :            : // *****************************************************************************
    1025                 :            : //  Compute the nodal extrema of ref el derivatives for chare-boundary nodes
    1026                 :            : //! \param[in] ncomp Number of conservative variables
    1027                 :            : //! \param[in] nprim Number of primitive variables
    1028                 :            : //! \param[in] ndof_NodalExtrm Degree of freedom for nodal extrema
    1029                 :            : //! \param[in] bndel List of elements contributing to chare-boundary nodes
    1030                 :            : //! \param[in] inpoel Element-node connectivity for element e
    1031                 :            : //! \param[in] gid Local->global node id map
    1032                 :            : //! \param[in] bid Local chare-boundary node ids (value) associated to
    1033                 :            : //!   global node ids (key)
    1034                 :            : //! \param[in] U Vector of conservative variables
    1035                 :            : //! \param[in] P Vector of primitive variables
    1036                 :            : //! \param[in,out] uNodalExtrm Chare-boundary nodal extrema for conservative
    1037                 :            : //!   variables
    1038                 :            : //! \param[in,out] pNodalExtrm Chare-boundary nodal extrema for primitive
    1039                 :            : //!   variables
    1040                 :            : // *****************************************************************************
    1041                 :            : {
    1042                 :       7110 :   const auto rdof = g_inputdeck.get< tag::rdof >();
    1043                 :            : 
    1044         [ +  + ]:     653190 :   for (auto e : bndel)
    1045                 :            :   {
    1046                 :            :     // access node IDs
    1047                 :            :     const std::vector<std::size_t> N
    1048         [ +  - ]:    1292160 :       { inpoel[e*4+0], inpoel[e*4+1], inpoel[e*4+2], inpoel[e*4+3] };
    1049                 :            : 
    1050                 :            :     // Loop over nodes of element e
    1051         [ +  + ]:    3230400 :     for(std::size_t ip=0; ip<4; ++ip)
    1052                 :            :     {
    1053         [ +  - ]:    2584320 :       auto i = bid.find( gid[N[ip]] );
    1054         [ +  + ]:    2584320 :       if (i != end(bid))      // If ip is the chare boundary point
    1055                 :            :       {
    1056                 :            :         // If DG(P2) is applied, find the nodal extrema of the gradients of
    1057                 :            :         // conservative/primitive variables in the reference element
    1058                 :            : 
    1059                 :            :         // Vector used to store the first order derivatives for both
    1060                 :            :         // conservative and primitive variables
    1061         [ +  - ]:    3509460 :         std::vector< std::array< tk::real, 3 > > gradc(ncomp, {0.0, 0.0, 0.0});
    1062         [ +  - ]:    3509460 :         std::vector< std::array< tk::real, 3 > > gradp(ncomp, {0.0, 0.0, 0.0});
    1063                 :            : 
    1064                 :            :         // Derivatives of the Dubiner basis
    1065                 :    1754730 :         std::array< tk::real, 3 > center {{0.25, 0.25, 0.25}};
    1066         [ +  - ]:    3509460 :         auto dBdxi = tk::eval_dBdxi(rdof, center);
    1067                 :            : 
    1068                 :            :         // Evaluate the first order derivative
    1069         [ +  + ]:   10528380 :         for(std::size_t icomp = 0; icomp < ncomp; icomp++)
    1070                 :            :         {
    1071                 :    8773650 :           auto mark = icomp * rdof;
    1072         [ +  + ]:   35094600 :           for(std::size_t idir = 0; idir < 3; idir++)
    1073                 :            :           {
    1074                 :   26320950 :             gradc[icomp][idir] = 0;
    1075         [ +  + ]:  263209500 :             for(std::size_t idof = 1; idof < rdof; idof++)
    1076         [ +  - ]:  236888550 :               gradc[icomp][idir] += U(e, mark+idof) * dBdxi[idir][idof];
    1077                 :            :           }
    1078                 :            :         }
    1079         [ -  + ]:    1754730 :         for(std::size_t icomp = 0; icomp < nprim; icomp++)
    1080                 :            :         {
    1081                 :          0 :           auto mark = icomp * rdof;
    1082         [ -  - ]:          0 :           for(std::size_t idir = 0; idir < 3; idir++)
    1083                 :            :           {
    1084                 :          0 :             gradp[icomp][idir] = 0;
    1085         [ -  - ]:          0 :             for(std::size_t idof = 1; idof < rdof; idof++)
    1086         [ -  - ]:          0 :               gradp[icomp][idir] += P(e, mark+idof) * dBdxi[idir][idof];
    1087                 :            :           }
    1088                 :            :         }
    1089                 :            : 
    1090                 :            :         // Store the extrema for the gradients
    1091         [ +  + ]:   10528380 :         for (std::size_t c=0; c<ncomp; ++c)
    1092                 :            :         {
    1093         [ +  + ]:   35094600 :           for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
    1094                 :            :           {
    1095                 :   26320950 :             auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
    1096                 :   26320950 :             auto min_mark = max_mark + 1;
    1097                 :   26320950 :             auto& ex = uNodalExtrm[i->second];
    1098                 :   26320950 :             ex[max_mark] = std::max(ex[max_mark], gradc[c][idof]);
    1099                 :   26320950 :             ex[min_mark] = std::min(ex[min_mark], gradc[c][idof]);
    1100                 :            :           }
    1101                 :            :         }
    1102         [ -  + ]:    1754730 :         for (std::size_t c=0; c<nprim; ++c)
    1103                 :            :         {
    1104         [ -  - ]:          0 :           for (std::size_t idof = 0; idof < ndof_NodalExtrm; idof++)
    1105                 :            :           {
    1106                 :          0 :             auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
    1107                 :          0 :             auto min_mark = max_mark + 1;
    1108                 :          0 :             auto& ex = pNodalExtrm[i->second];
    1109                 :          0 :             ex[max_mark] = std::max(ex[max_mark], gradp[c][idof]);
    1110                 :          0 :             ex[min_mark] = std::min(ex[min_mark], gradp[c][idof]);
    1111                 :            :           }
    1112                 :            :         }
    1113                 :            :       }
    1114                 :            :     }
    1115                 :            :   }
    1116                 :       7110 : }
    1117                 :            : 
    1118                 :            : void
    1119                 :      38265 : DG::lim()
    1120                 :            : // *****************************************************************************
    1121                 :            : // Compute limiter function
    1122                 :            : // *****************************************************************************
    1123                 :            : {
    1124                 :      38265 :   auto d = Disc();
    1125                 :      38265 :   const auto rdof = g_inputdeck.get< tag::rdof >();
    1126                 :      38265 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
    1127                 :      38265 :   const auto ncomp = m_u.nprop() / rdof;
    1128                 :      38265 :   const auto nprim = m_p.nprop() / rdof;
    1129                 :            : 
    1130                 :            :   // Combine own and communicated contributions to nodal extrema
    1131         [ +  + ]:    1317765 :   for (const auto& [gid,g] : m_uNodalExtrmc) {
    1132         [ +  - ]:    1279500 :     auto bid = tk::cref_find( d->Bid(), gid );
    1133         [ +  + ]:    9304560 :     for (ncomp_t c=0; c<ncomp; ++c)
    1134                 :            :     {
    1135         [ +  + ]:   32100240 :       for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
    1136                 :            :       {
    1137                 :   24075180 :         auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
    1138                 :   24075180 :         auto min_mark = max_mark + 1;
    1139                 :   24075180 :         m_uNodalExtrm[bid][max_mark] =
    1140                 :   24075180 :           std::max(g[max_mark], m_uNodalExtrm[bid][max_mark]);
    1141                 :   24075180 :         m_uNodalExtrm[bid][min_mark] =
    1142                 :   24075180 :           std::min(g[min_mark], m_uNodalExtrm[bid][min_mark]);
    1143                 :            :       }
    1144                 :            :     }
    1145                 :            :   }
    1146         [ +  + ]:    1317765 :   for (const auto& [gid,g] : m_pNodalExtrmc) {
    1147         [ +  - ]:    1279500 :     auto bid = tk::cref_find( d->Bid(), gid );
    1148         [ +  + ]:    3460350 :     for (ncomp_t c=0; c<nprim; ++c)
    1149                 :            :     {
    1150         [ +  + ]:    8723400 :       for(std::size_t idof=0; idof<m_ndof_NodalExtrm; idof++)
    1151                 :            :       {
    1152                 :    6542550 :         auto max_mark = 2*c*m_ndof_NodalExtrm + 2*idof;
    1153                 :    6542550 :         auto min_mark = max_mark + 1;
    1154                 :    6542550 :         m_pNodalExtrm[bid][max_mark] =
    1155                 :    6542550 :           std::max(g[max_mark], m_pNodalExtrm[bid][max_mark]);
    1156                 :    6542550 :         m_pNodalExtrm[bid][min_mark] =
    1157                 :    6542550 :           std::min(g[min_mark], m_pNodalExtrm[bid][min_mark]);
    1158                 :            :       }
    1159                 :            :     }
    1160                 :            :   }
    1161                 :            : 
    1162                 :            :   // clear gradients receive buffer
    1163                 :      38265 :   tk::destroy(m_uNodalExtrmc);
    1164                 :      38265 :   tk::destroy(m_pNodalExtrmc);
    1165                 :            : 
    1166         [ +  + ]:      38265 :   if (rdof > 1) {
    1167                 :      74520 :     g_dgpde[d->MeshId()].limit( d->T(), pref, myGhosts()->m_geoFace,
    1168                 :      24840 :               myGhosts()->m_geoElem, myGhosts()->m_fd, myGhosts()->m_esup,
    1169                 :      24840 :               myGhosts()->m_inpoel, myGhosts()->m_coord, m_ndof, d->Gid(),
    1170                 :      24840 :               d->Bid(), m_uNodalExtrm, m_pNodalExtrm, m_mtInv, m_u, m_p,
    1171                 :      24840 :               m_shockmarker );
    1172                 :            : 
    1173         [ +  + ]:      24840 :     if (g_inputdeck.get< tag::limsol_projection >())
    1174                 :      41580 :       g_dgpde[d->MeshId()].CPL(m_p, myGhosts()->m_geoElem,
    1175                 :      20790 :         myGhosts()->m_inpoel, myGhosts()->m_coord, m_u,
    1176                 :      20790 :         myGhosts()->m_fd.Esuel().size()/4);
    1177                 :            :   }
    1178                 :            : 
    1179                 :            :   // Send limited solution to neighboring chares
    1180         [ +  + ]:      38265 :   if (myGhosts()->m_sendGhost.empty())
    1181                 :       3390 :     comlim_complete();
    1182                 :            :   else
    1183 [ +  - ][ +  + ]:     419925 :     for(const auto& [cid, ghostdata] : myGhosts()->m_sendGhost) {
    1184         [ +  - ]:     770100 :       std::vector< std::size_t > tetid( ghostdata.size() );
    1185         [ +  - ]:     770100 :       std::vector< std::vector< tk::real > > u( ghostdata.size() ),
    1186         [ +  - ]:     770100 :                                              prim( ghostdata.size() );
    1187                 :     385050 :       std::vector< std::size_t > ndof;
    1188                 :     385050 :       std::size_t j = 0;
    1189         [ +  + ]:    6854730 :       for(const auto& i : ghostdata) {
    1190 [ +  - ][ -  + ]:    6469680 :         Assert( i < myGhosts()->m_fd.Esuel().size()/4,
         [ -  - ][ -  - ]
                 [ -  - ]
    1191                 :            :           "Sending limiter ghost data" );
    1192                 :    6469680 :         tetid[j] = i;
    1193         [ +  - ]:    6469680 :         u[j] = m_u[i];
    1194         [ +  - ]:    6469680 :         prim[j] = m_p[i];
    1195                 :    6469680 :         ++j;
    1196                 :            :       }
    1197 [ +  - ][ +  - ]:     385050 :       thisProxy[ cid ].comlim( thisIndex, tetid, u, prim );
    1198                 :            :     }
    1199                 :            : 
    1200                 :      38265 :   ownlim_complete();
    1201                 :      38265 : }
    1202                 :            : 
    1203                 :            : void
    1204                 :       1770 : DG::refine_ndof()
    1205                 :            : // *****************************************************************************
    1206                 :            : //  p-refine all elements that are adjacent to p-refined elements
    1207                 :            : //! \details This function p-refines all the neighbors of an element that has
    1208                 :            : //!   been p-refined as a result of an error indicator.
    1209                 :            : // *****************************************************************************
    1210                 :            : {
    1211         [ +  - ]:       1770 :   auto d = Disc();
    1212                 :       1770 :   const auto& coord = d->Coord();
    1213                 :       1770 :   const auto& inpoel = d->Inpoel();
    1214                 :       1770 :   const auto npoin = coord[0].size();
    1215         [ +  - ]:       1770 :   const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
    1216         [ +  - ]:       3540 :   std::vector<std::size_t> node_ndof(npoin, 1);
    1217                 :            : 
    1218                 :            :   // Mark the max ndof for each node and store in node_ndof
    1219         [ +  + ]:     188540 :   for(std::size_t ip=0; ip<npoin; ip++)
    1220                 :            :   {
    1221 [ +  - ][ +  - ]:     186770 :     const auto& pesup = tk::cref_find(myGhosts()->m_esup, ip);
    1222         [ +  + ]:    2641810 :     for(auto er : pesup)
    1223                 :    2455040 :       node_ndof[ip] = std::max(m_ndof[er], node_ndof[ip]);
    1224                 :            :   }
    1225                 :            : 
    1226         [ +  + ]:     414490 :   for(std::size_t e = 0; e < nelem; e++)
    1227                 :            :   {
    1228                 :            :     // Find if any node of this element has p1/p2 ndofs
    1229                 :     412720 :     std::size_t counter_p2(0);
    1230                 :     412720 :     std::size_t counter_p1(0);
    1231         [ +  + ]:    2063600 :     for(std::size_t inode = 0; inode < 4; inode++)
    1232                 :            :     {
    1233                 :    1650880 :       auto node = inpoel[4*e+inode];
    1234         [ +  + ]:    1650880 :       if(node_ndof[node] == 10)
    1235                 :     320808 :         counter_p2++;
    1236         [ +  + ]:    1330072 :       else if (node_ndof[node] == 4)
    1237                 :     180140 :         counter_p1++;
    1238                 :            :     }
    1239                 :            : 
    1240                 :            :     // If there is at least one node with p1/p2 ndofs, all of the elements
    1241                 :            :     // around this node are refined to p1/p2.
    1242 [ +  + ][ +  + ]:     412720 :     if(counter_p2 > 0 && m_ndof[e] < 10)
                 [ +  + ]
    1243                 :            :     {
    1244         [ +  + ]:      16717 :       if(m_ndof[e] == 4)
    1245                 :      15693 :         m_ndof[e] = 10;
    1246         [ +  + ]:      16717 :       if(m_ndof[e] == 1)
    1247                 :       1024 :         m_ndof[e] = 4;
    1248                 :            :     }
    1249 [ +  + ][ +  + ]:     396003 :     else if(counter_p1 > 0 && m_ndof[e] < 4)
                 [ +  + ]
    1250                 :      13452 :       m_ndof[e] = 4;
    1251                 :            :   }
    1252                 :       1770 : }
    1253                 :            : 
    1254                 :       1770 : void DG::smooth_ndof()
    1255                 :            : // *****************************************************************************
    1256                 :            : //  Smooth the refined ndof distribution to avoid zigzag refinement
    1257                 :            : // *****************************************************************************
    1258                 :            : {
    1259         [ +  - ]:       1770 :   auto d = Disc();
    1260                 :       1770 :   const auto& inpoel = d->Inpoel();
    1261                 :       1770 :   const auto& coord = d->Coord();
    1262                 :       1770 :   const auto npoin = coord[0].size();
    1263         [ +  - ]:       1770 :   const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
    1264         [ +  - ]:       3540 :   std::vector<std::size_t> node_ndof(npoin, 1);
    1265                 :            : 
    1266                 :            :   // Mark the max ndof for each node and store in node_ndof
    1267         [ +  + ]:     188540 :   for(std::size_t ip=0; ip<npoin; ip++)
    1268                 :            :   {
    1269 [ +  - ][ +  - ]:     186770 :     const auto& pesup = tk::cref_find(myGhosts()->m_esup, ip);
    1270         [ +  + ]:    2641810 :     for(auto er : pesup)
    1271                 :    2455040 :       node_ndof[ip] = std::max(m_ndof[er], node_ndof[ip]);
    1272                 :            :   }
    1273                 :            : 
    1274         [ +  + ]:     414490 :   for(std::size_t e = 0; e < nelem; e++)
    1275                 :            :   {
    1276                 :            :     // Find if any node of this element has p1/p2 ndofs
    1277                 :     412720 :     std::size_t counter_p2(0);
    1278                 :     412720 :     std::size_t counter_p1(0);
    1279         [ +  + ]:    2063600 :     for(std::size_t inode = 0; inode < 4; inode++)
    1280                 :            :     {
    1281                 :    1650880 :       auto node = inpoel[4*e+inode];
    1282         [ +  + ]:    1650880 :       if(node_ndof[node] == 10)
    1283                 :     382503 :         counter_p2++;
    1284         [ +  + ]:    1268377 :       else if (node_ndof[node] == 4)
    1285                 :     182564 :         counter_p1++;
    1286                 :            :     }
    1287                 :            : 
    1288                 :            :     // If all the nodes in the element are p1/p2, this element is refined to
    1289                 :            :     // p1/p2.
    1290 [ +  + ][ +  + ]:     412720 :     if(counter_p2 == 4 && m_ndof[e] == 4)
                 [ +  + ]
    1291                 :       1469 :       m_ndof[e] = 10;
    1292 [ +  + ][ +  + ]:     411251 :     else if(counter_p1 == 4 && m_ndof[e] == 1)
                 [ +  + ]
    1293                 :       1553 :       m_ndof[e] = 4;
    1294                 :            :   }
    1295                 :       1770 : }
    1296                 :            : 
    1297                 :            : void
    1298                 :     385050 : DG::comlim( int fromch,
    1299                 :            :             const std::vector< std::size_t >& tetid,
    1300                 :            :             const std::vector< std::vector< tk::real > >& u,
    1301                 :            :             const std::vector< std::vector< tk::real > >& prim )
    1302                 :            : // *****************************************************************************
    1303                 :            : //  Receive chare-boundary limiter ghost data from neighboring chares
    1304                 :            : //! \param[in] fromch Sender chare id
    1305                 :            : //! \param[in] tetid Ghost tet ids we receive solution data for
    1306                 :            : //! \param[in] u Limited high-order solution
    1307                 :            : //! \param[in] prim Limited high-order primitive quantities
    1308                 :            : //! \details This function receives contributions to the limited solution from
    1309                 :            : //!   fellow chares.
    1310                 :            : // *****************************************************************************
    1311                 :            : {
    1312 [ -  + ][ -  - ]:     385050 :   Assert( u.size() == tetid.size(), "Size mismatch in DG::comlim()" );
         [ -  - ][ -  - ]
    1313 [ -  + ][ -  - ]:     385050 :   Assert( prim.size() == tetid.size(), "Size mismatch in DG::comlim()" );
         [ -  - ][ -  - ]
    1314                 :            : 
    1315                 :            :   // Find local-to-ghost tet id map for sender chare
    1316                 :     385050 :   const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
    1317                 :            : 
    1318         [ +  + ]:    6854730 :   for (std::size_t i=0; i<tetid.size(); ++i) {
    1319         [ +  - ]:    6469680 :     auto j = tk::cref_find( n, tetid[i] );
    1320 [ +  - ][ -  + ]:    6469680 :     Assert( j >= myGhosts()->m_fd.Esuel().size()/4,
         [ -  - ][ -  - ]
                 [ -  - ]
    1321                 :            :       "Receiving solution non-ghost data" );
    1322 [ +  - ][ +  - ]:    6469680 :     auto b = tk::cref_find( myGhosts()->m_bid, j );
    1323 [ -  + ][ -  - ]:    6469680 :     Assert( b < m_uc[2].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
    1324 [ -  + ][ -  - ]:    6469680 :     Assert( b < m_pc[2].size(), "Indexing out of bounds" );
         [ -  - ][ -  - ]
    1325         [ +  - ]:    6469680 :     m_uc[2][b] = u[i];
    1326         [ +  - ]:    6469680 :     m_pc[2][b] = prim[i];
    1327                 :            :   }
    1328                 :            : 
    1329                 :            :   // if we have received all solution ghost contributions from neighboring
    1330                 :            :   // chares (chares we communicate along chare-boundary faces with), and
    1331                 :            :   // contributed our solution to these neighbors, proceed to limiting
    1332         [ +  + ]:     385050 :   if (++m_nlim == myGhosts()->m_sendGhost.size()) {
    1333                 :      34875 :     m_nlim = 0;
    1334                 :      34875 :     comlim_complete();
    1335                 :            :   }
    1336                 :     385050 : }
    1337                 :            : 
    1338                 :            : void
    1339                 :      38265 : DG::dt()
    1340                 :            : // *****************************************************************************
    1341                 :            : // Compute time step size
    1342                 :            : // *****************************************************************************
    1343                 :            : {
    1344         [ +  - ]:      38265 :   auto d = Disc();
    1345                 :            : 
    1346                 :            :   // Combine own and communicated contributions of limited solution and degrees
    1347                 :            :   // of freedom in cells (if p-adaptive)
    1348 [ +  - ][ +  + ]:    6507945 :   for (const auto& b : myGhosts()->m_bid) {
    1349 [ -  + ][ -  - ]:    6469680 :     Assert( m_uc[2][b.second].size() == m_u.nprop(), "ncomp size mismatch" );
         [ -  - ][ -  - ]
    1350 [ -  + ][ -  - ]:    6469680 :     Assert( m_pc[2][b.second].size() == m_p.nprop(), "ncomp size mismatch" );
         [ -  - ][ -  - ]
    1351         [ +  + ]:  173340045 :     for (std::size_t c=0; c<m_u.nprop(); ++c) {
    1352         [ +  - ]:  166870365 :       m_u(b.first,c) = m_uc[2][b.second][c];
    1353                 :            :     }
    1354         [ +  + ]:   34324605 :     for (std::size_t c=0; c<m_p.nprop(); ++c) {
    1355         [ +  - ]:   27854925 :       m_p(b.first,c) = m_pc[2][b.second][c];
    1356                 :            :     }
    1357                 :            :   }
    1358                 :            : 
    1359                 :      38265 :   auto mindt = std::numeric_limits< tk::real >::max();
    1360                 :            : 
    1361         [ +  + ]:      38265 :   if (m_stage == 0)
    1362                 :            :   {
    1363                 :      12755 :     auto const_dt = g_inputdeck.get< tag::dt >();
    1364                 :      12755 :     auto eps = std::numeric_limits< tk::real >::epsilon();
    1365                 :            : 
    1366                 :            :     // use constant dt if configured
    1367         [ +  + ]:      12755 :     if (std::abs(const_dt) > eps) {
    1368                 :            : 
    1369                 :       9925 :       mindt = const_dt;
    1370                 :            : 
    1371                 :            :     } else {      // compute dt based on CFL
    1372                 :            : 
    1373                 :            :       // find the minimum dt across all PDEs integrated
    1374                 :            :       auto eqdt =
    1375 [ +  - ][ +  - ]:       5660 :         g_dgpde[d->MeshId()].dt( myGhosts()->m_coord, myGhosts()->m_inpoel,
    1376         [ +  - ]:       2830 :           myGhosts()->m_fd,
    1377 [ +  - ][ +  - ]:       2830 :           myGhosts()->m_geoFace, myGhosts()->m_geoElem, m_ndof, m_u, m_p,
                 [ +  - ]
    1378         [ +  - ]:       2830 :           myGhosts()->m_fd.Esuel().size()/4 );
    1379         [ +  - ]:       2830 :       if (eqdt < mindt) mindt = eqdt;
    1380                 :            : 
    1381                 :            :       // time-step suppression for unsteady problems
    1382                 :       2830 :       tk::real coeff(1.0);
    1383 [ -  + ][ -  - ]:       2830 :       if (g_inputdeck.get< tag::cfl_ramping >() && d->It() < 100) coeff = 0.01 * static_cast< tk::real >(d->It()+1);
                 [ -  + ]
    1384                 :            : 
    1385                 :       2830 :       mindt *= coeff * g_inputdeck.get< tag::cfl >();
    1386                 :            :     }
    1387                 :            :   }
    1388                 :            :   else
    1389                 :            :   {
    1390                 :      25510 :     mindt = d->Dt();
    1391                 :            :   }
    1392                 :            : 
    1393                 :            :   // Resize the buffer vector of nodal extrema
    1394         [ +  - ]:      38265 :   resizeNodalExtremac();
    1395                 :            : 
    1396                 :            :   // Contribute to minimum dt across all chares then advance to next step
    1397         [ +  - ]:      38265 :   contribute( sizeof(tk::real), &mindt, CkReduction::min_double,
    1398 [ +  - ][ +  - ]:      76530 :               CkCallback(CkReductionTarget(DG,solve), thisProxy) );
    1399                 :      38265 : }
    1400                 :            : 
    1401                 :            : void
    1402                 :      38265 : DG::solve( tk::real newdt )
    1403                 :            : // *****************************************************************************
    1404                 :            : // Compute right-hand side of discrete transport equations
    1405                 :            : //! \param[in] newdt Size of this new time step
    1406                 :            : // *****************************************************************************
    1407                 :            : {
    1408                 :      38265 :   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
    1409                 :            : 
    1410                 :            :   // Enable SDAG wait for building the solution vector during the next stage
    1411         [ +  - ]:      38265 :   thisProxy[ thisIndex ].wait4sol();
    1412 [ +  + ][ +  - ]:      38265 :   if (pref) thisProxy[ thisIndex ].wait4refine();
    1413         [ +  - ]:      38265 :   thisProxy[ thisIndex ].wait4smooth();
    1414         [ +  - ]:      38265 :   thisProxy[ thisIndex ].wait4reco();
    1415         [ +  - ]:      38265 :   thisProxy[ thisIndex ].wait4nodalExtrema();
    1416         [ +  - ]:      38265 :   thisProxy[ thisIndex ].wait4lim();
    1417         [ +  - ]:      38265 :   thisProxy[ thisIndex ].wait4nod();
    1418                 :            : 
    1419                 :      38265 :   auto d = Disc();
    1420                 :      38265 :   const auto rdof = g_inputdeck.get< tag::rdof >();
    1421                 :      38265 :   const auto ndof = g_inputdeck.get< tag::ndof >();
    1422                 :      38265 :   const auto neq = m_u.nprop()/rdof;
    1423                 :            : 
    1424                 :            :   // Set new time step size
    1425         [ +  + ]:      38265 :   if (m_stage == 0) d->setdt( newdt );
    1426                 :            : 
    1427                 :            :   // Update Un
    1428         [ +  + ]:      38265 :   if (m_stage == 0) m_un = m_u;
    1429                 :            : 
    1430                 :            :   // Explicit or IMEX
    1431                 :      38265 :   const auto imex_runge_kutta = g_inputdeck.get< tag::imex_runge_kutta >();
    1432                 :            : 
    1433                 :            :   // physical time at time-stage for computing exact source terms
    1434                 :      38265 :   tk::real physT(d->T());
    1435         [ +  + ]:      38265 :   if (m_stage == 1) {
    1436                 :      12755 :     physT += d->Dt();
    1437                 :            :   }
    1438         [ +  + ]:      25510 :   else if (m_stage == 2) {
    1439                 :      12755 :     physT += 0.5*d->Dt();
    1440                 :            :   }
    1441                 :            : 
    1442         [ -  + ]:      38265 :   if (imex_runge_kutta) {
    1443         [ -  - ]:          0 :     if (m_stage == 0)
    1444                 :            :     {
    1445                 :            :       // Save previous rhs
    1446                 :          0 :       m_rhsprev = m_rhs;
    1447                 :            :       // Initialize m_stiffrhs to zero
    1448                 :          0 :       m_stiffrhs.fill(0.0);
    1449                 :          0 :       m_stiffrhsprev.fill(0.0);
    1450                 :            :     }
    1451                 :            :   }
    1452                 :            : 
    1453                 :      76530 :   g_dgpde[d->MeshId()].rhs( physT, pref, myGhosts()->m_geoFace,
    1454                 :      38265 :     myGhosts()->m_geoElem, myGhosts()->m_fd, myGhosts()->m_inpoel, m_boxelems,
    1455                 :      38265 :     myGhosts()->m_coord, m_u, m_p, m_ndof, d->Dt(), m_rhs );
    1456                 :            : 
    1457         [ +  - ]:      38265 :   if (!imex_runge_kutta) {
    1458                 :            :     // Explicit time-stepping using RK3 to discretize time-derivative
    1459         [ +  + ]:   16351845 :     for(std::size_t e=0; e<myGhosts()->m_nunk; ++e) {
    1460                 :   16313580 :       auto vole = myGhosts()->m_geoElem(e,0);
    1461         [ +  + ]:  111412515 :       for(std::size_t c=0; c<neq; ++c)
    1462                 :            :       {
    1463         [ +  + ]:  373799400 :         for (std::size_t k=0; k<m_numEqDof[c]; ++k)
    1464                 :            :         {
    1465         [ +  + ]:  278700465 :           if(k < m_ndof[e]) {
    1466                 :  193948005 :             auto rmark = c*rdof+k;
    1467                 :  193948005 :             auto mark = c*ndof+k;
    1468                 :  193948005 :             m_u(e, rmark) =  rkcoef[0][m_stage] * m_un(e, rmark)
    1469                 :  193948005 :               + rkcoef[1][m_stage] * ( m_u(e, rmark)
    1470                 :  193948005 :                 + d->Dt() * m_rhs(e, mark)/ (vole*mass_dubiner[k]));
    1471         [ +  + ]:  193948005 :             if(fabs(m_u(e, rmark)) < 1e-16)
    1472                 :   31380279 :               m_u(e, rmark) = 0;
    1473                 :            :           }
    1474                 :            :         }
    1475                 :            :       }
    1476                 :            :     }
    1477                 :            :   }
    1478                 :            :   else {
    1479                 :            :     // Implicit-Explicit time-stepping using RK3 to discretize time-derivative
    1480                 :          0 :     DG::imex_integrate();
    1481                 :            :   }
    1482                 :            : 
    1483         [ +  + ]:   16351845 :   for(std::size_t e=0; e<myGhosts()->m_nunk; ++e)
    1484         [ +  + ]:  111412515 :     for(std::size_t c=0; c<neq; ++c)
    1485                 :            :     {
    1486                 :            :       // zero out unused/reconstructed dofs of equations using reduced dofs
    1487                 :            :       // (see DGMultiMat::numEquationDofs())
    1488         [ +  + ]:   95098935 :       if (m_numEqDof[c] < rdof) {
    1489         [ +  + ]:   96926640 :         for (std::size_t k=m_numEqDof[c]; k<rdof; ++k)
    1490                 :            :         {
    1491                 :   72694980 :           auto rmark = c*rdof+k;
    1492                 :   72694980 :           m_u(e, rmark) = 0.0;
    1493                 :            :         }
    1494                 :            :       }
    1495                 :            :     }
    1496                 :            : 
    1497                 :            :   // Update primitives based on the evolved solution
    1498                 :      76530 :   g_dgpde[d->MeshId()].updateInterfaceCells( m_u,
    1499                 :      38265 :     myGhosts()->m_fd.Esuel().size()/4, m_ndof, m_interface );
    1500                 :      76530 :   g_dgpde[d->MeshId()].updatePrimitives( m_u, myGhosts()->m_geoElem, m_p,
    1501                 :      38265 :     myGhosts()->m_fd.Esuel().size()/4, m_ndof );
    1502         [ +  - ]:      38265 :   if (!g_inputdeck.get< tag::accuracy_test >()) {
    1503                 :      76530 :     g_dgpde[d->MeshId()].cleanTraceMaterial( physT, myGhosts()->m_geoElem, m_u,
    1504                 :      38265 :       m_p, myGhosts()->m_fd.Esuel().size()/4 );
    1505                 :            :   }
    1506                 :            : 
    1507         [ +  + ]:      38265 :   if (m_stage < m_nstage-1) {
    1508                 :            : 
    1509                 :            :     // continue with next time step stage
    1510                 :      25510 :     stage();
    1511                 :            : 
    1512                 :            :   } else {
    1513                 :            : 
    1514                 :            :     // Increase number of iterations and physical time
    1515                 :      12755 :     d->next();
    1516                 :            : 
    1517                 :            :     // Compute diagnostics, e.g., residuals
    1518                 :      38265 :     auto diag_computed = m_diag.compute( *d,
    1519                 :      12755 :       m_u.nunk()-myGhosts()->m_fd.Esuel().size()/4, myGhosts()->m_geoElem,
    1520                 :      12755 :       m_ndof, m_u, m_un );
    1521                 :            : 
    1522                 :            :     // Continue to mesh refinement (if configured)
    1523 [ +  + ][ +  - ]:      12755 :     if (!diag_computed) refine( std::vector< tk::real >( m_u.nprop(), 0.0 ) );
                 [ +  - ]
    1524                 :            : 
    1525                 :            :   }
    1526                 :      38265 : }
    1527                 :            : 
    1528                 :            : void
    1529                 :      12755 : DG::refine( [[maybe_unused]] const std::vector< tk::real >& l2res )
    1530                 :            : // *****************************************************************************
    1531                 :            : // Optionally refine/derefine mesh
    1532                 :            : //! \param[in] l2res L2-norms of the residual for each scalar component
    1533                 :            : //!   computed across the whole problem
    1534                 :            : // *****************************************************************************
    1535                 :            : {
    1536                 :      12755 :   auto d = Disc();
    1537                 :            : 
    1538                 :      12755 :   auto dtref = g_inputdeck.get< tag::amr, tag::dtref >();
    1539                 :      12755 :   auto dtfreq = g_inputdeck.get< tag::amr, tag::dtfreq >();
    1540                 :            : 
    1541                 :            :   // if t>0 refinement enabled and we hit the dtref frequency
    1542 [ -  + ][ -  - ]:      12755 :   if (dtref && !(d->It() % dtfreq)) {   // refine
                 [ -  + ]
    1543                 :            : 
    1544                 :          0 :     d->startvol();
    1545 [ -  - ][ -  - ]:          0 :     d->Ref()->dtref( myGhosts()->m_fd.Bface(), {},
    1546                 :          0 :       tk::remap(myGhosts()->m_fd.Triinpoel(),d->Gid()) );
    1547                 :          0 :     d->refined() = 1;
    1548                 :            : 
    1549                 :            :   } else {      // do not refine
    1550                 :            : 
    1551                 :      12755 :     d->refined() = 0;
    1552                 :      12755 :     stage();
    1553                 :            : 
    1554                 :            :   }
    1555                 :      12755 : }
    1556                 :            : 
    1557                 :            : void
    1558                 :          0 : DG::resizePostAMR(
    1559                 :            :   const std::vector< std::size_t >& /*ginpoel*/,
    1560                 :            :   const tk::UnsMesh::Chunk& chunk,
    1561                 :            :   const tk::UnsMesh::Coords& coord,
    1562                 :            :   const std::unordered_map< std::size_t, tk::UnsMesh::Edge >& /*addedNodes*/,
    1563                 :            :   const std::unordered_map< std::size_t, std::size_t >& addedTets,
    1564                 :            :   const std::set< std::size_t >& removedNodes,
    1565                 :            :   const std::unordered_map< std::size_t, std::size_t >& amrNodeMap,
    1566                 :            :   const tk::NodeCommMap& nodeCommMap,
    1567                 :            :   const std::map< int, std::vector< std::size_t > >& bface,
    1568                 :            :   const std::map< int, std::vector< std::size_t > >& /* bnode */,
    1569                 :            :   const std::vector< std::size_t >& triinpoel,
    1570                 :            :   const std::unordered_map< std::size_t, std::set< std::size_t > >& elemblockid )
    1571                 :            : // *****************************************************************************
    1572                 :            : //  Receive new mesh from Refiner
    1573                 :            : //! \param[in] chunk New mesh chunk (connectivity and global<->local id maps)
    1574                 :            : //! \param[in] coord New mesh node coordinates
    1575                 :            : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
    1576                 :            : //! \param[in] removedNodes Newly removed mesh node local ids
    1577                 :            : //! \param[in] amrNodeMap Node id map after amr (local ids)
    1578                 :            : //! \param[in] nodeCommMap New node communication map
    1579                 :            : //! \param[in] bface Boundary-faces mapped to side set ids
    1580                 :            : //! \param[in] triinpoel Boundary-face connectivity
    1581                 :            : //! \param[in] elemblockid Local tet ids associated with mesh block ids
    1582                 :            : // *****************************************************************************
    1583                 :            : {
    1584         [ -  - ]:          0 :   auto d = Disc();
    1585                 :            : 
    1586                 :            :   // Set flag that indicates that we are during time stepping
    1587                 :          0 :   m_initial = 0;
    1588         [ -  - ]:          0 :   myGhosts()->m_initial = 0;
    1589                 :            : 
    1590                 :            :   // Zero field output iteration count between two mesh refinement steps
    1591                 :          0 :   d->Itf() = 0;
    1592                 :            : 
    1593                 :            :   // Increase number of iterations with mesh refinement
    1594                 :          0 :   ++d->Itr();
    1595                 :            : 
    1596                 :            :   // Save old number of elements
    1597         [ -  - ]:          0 :   [[maybe_unused]] auto old_nelem = myGhosts()->m_inpoel.size()/4;
    1598                 :            : 
    1599                 :            :   // Resize mesh data structures
    1600         [ -  - ]:          0 :   d->resizePostAMR( chunk, coord, amrNodeMap, nodeCommMap, removedNodes,
    1601                 :            :     elemblockid );
    1602                 :            : 
    1603                 :            :   // Update state
    1604 [ -  - ][ -  - ]:          0 :   myGhosts()->m_inpoel = d->Inpoel();
    1605 [ -  - ][ -  - ]:          0 :   myGhosts()->m_coord = d->Coord();
    1606         [ -  - ]:          0 :   auto nelem = myGhosts()->m_inpoel.size()/4;
    1607         [ -  - ]:          0 :   m_p.resize( nelem );
    1608         [ -  - ]:          0 :   m_u.resize( nelem );
    1609         [ -  - ]:          0 :   m_un.resize( nelem );
    1610         [ -  - ]:          0 :   m_rhs.resize( nelem );
    1611         [ -  - ]:          0 :   m_rhsprev.resize( nelem );
    1612         [ -  - ]:          0 :   m_stiffrhs.resize( nelem );
    1613         [ -  - ]:          0 :   m_stiffrhsprev.resize( nelem );
    1614 [ -  - ][ -  - ]:          0 :   m_uNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
    1615         [ -  - ]:          0 :     m_ndof_NodalExtrm*g_inputdeck.get< tag::ncomp >() ) );
    1616 [ -  - ][ -  - ]:          0 :   m_pNodalExtrm.resize( Disc()->Bid().size(), std::vector<tk::real>( 2*
    1617         [ -  - ]:          0 :     m_ndof_NodalExtrm*m_p.nprop()/g_inputdeck.get< tag::rdof >()));
    1618                 :            : 
    1619                 :            :   // Resize the buffer vector of nodal extrema
    1620         [ -  - ]:          0 :   resizeNodalExtremac();
    1621                 :            : 
    1622 [ -  - ][ -  - ]:          0 :   myGhosts()->m_fd = FaceData( myGhosts()->m_inpoel, bface,
                 [ -  - ]
    1623         [ -  - ]:          0 :     tk::remap(triinpoel,d->Lid()) );
    1624                 :            : 
    1625         [ -  - ]:          0 :   myGhosts()->m_geoFace =
    1626 [ -  - ][ -  - ]:          0 :     tk::Fields( tk::genGeoFaceTri( myGhosts()->m_fd.Nipfac(),
    1627         [ -  - ]:          0 :     myGhosts()->m_fd.Inpofa(), coord ) );
    1628 [ -  - ][ -  - ]:          0 :   myGhosts()->m_geoElem = tk::Fields( tk::genGeoElemTet( myGhosts()->m_inpoel,
                 [ -  - ]
    1629                 :          0 :     coord ) );
    1630                 :            : 
    1631 [ -  - ][ -  - ]:          0 :   myGhosts()->m_nfac = myGhosts()->m_fd.Inpofa().size()/3;
    1632         [ -  - ]:          0 :   myGhosts()->m_nunk = nelem;
    1633                 :          0 :   m_npoin = coord[0].size();
    1634         [ -  - ]:          0 :   myGhosts()->m_bndFace.clear();
    1635         [ -  - ]:          0 :   myGhosts()->m_exptGhost.clear();
    1636         [ -  - ]:          0 :   myGhosts()->m_sendGhost.clear();
    1637         [ -  - ]:          0 :   myGhosts()->m_ghost.clear();
    1638         [ -  - ]:          0 :   myGhosts()->m_esup.clear();
    1639                 :            : 
    1640                 :            :   // Update solution on new mesh, P0 (cell center value) only for now
    1641         [ -  - ]:          0 :   m_un = m_u;
    1642         [ -  - ]:          0 :   auto pn = m_p;
    1643                 :          0 :   auto unprop = m_u.nprop();
    1644                 :          0 :   auto pnprop = m_p.nprop();
    1645         [ -  - ]:          0 :   for (const auto& [child,parent] : addedTets) {
    1646 [ -  - ][ -  - ]:          0 :     Assert( child < nelem, "Indexing out of new solution vector" );
         [ -  - ][ -  - ]
    1647 [ -  - ][ -  - ]:          0 :     Assert( parent < old_nelem, "Indexing out of old solution vector" );
         [ -  - ][ -  - ]
    1648 [ -  - ][ -  - ]:          0 :     for (std::size_t i=0; i<unprop; ++i) m_u(child,i) = m_un(parent,i);
                 [ -  - ]
    1649 [ -  - ][ -  - ]:          0 :     for (std::size_t i=0; i<pnprop; ++i) m_p(child,i) = pn(parent,i);
                 [ -  - ]
    1650                 :            :   }
    1651         [ -  - ]:          0 :   m_un = m_u;
    1652                 :            : 
    1653                 :            :   // Resize communication buffers
    1654 [ -  - ][ -  - ]:          0 :   m_ghosts[thisIndex].resizeComm();
    1655                 :          0 : }
    1656                 :            : 
    1657                 :            : bool
    1658                 :      10214 : DG::fieldOutput() const
    1659                 :            : // *****************************************************************************
    1660                 :            : // Decide wether to output field data
    1661                 :            : //! \return True if field data is output in this step
    1662                 :            : // *****************************************************************************
    1663                 :            : {
    1664                 :      10214 :   auto d = Disc();
    1665                 :            : 
    1666                 :            :   // Output field data
    1667 [ +  + ][ +  - ]:      10214 :   return d->fielditer() or d->fieldtime() or d->fieldrange() or d->finished();
         [ +  - ][ +  + ]
    1668                 :            : }
    1669                 :            : 
    1670                 :            : bool
    1671                 :       1325 : DG::refinedOutput() const
    1672                 :            : // *****************************************************************************
    1673                 :            : // Decide if we write field output using a refined mesh
    1674                 :            : //! \return True if field output will use a refined mesh
    1675                 :            : // *****************************************************************************
    1676                 :            : {
    1677         [ +  + ]:       1358 :   return g_inputdeck.get< tag::field_output, tag::refined >() &&
    1678         [ +  - ]:       1358 :          g_inputdeck.get< tag::scheme >() != ctr::SchemeType::DGP0;
    1679                 :            : }
    1680                 :            : 
    1681                 :            : void
    1682                 :       1325 : DG::writeFields(
    1683                 :            :   CkCallback c,
    1684                 :            :   const std::unordered_map< std::size_t, std::size_t >& addedTets )
    1685                 :            : // *****************************************************************************
    1686                 :            : // Output mesh field data
    1687                 :            : //! \param[in] c Function to continue with after the write
    1688                 :            : //! \param[in] addedTets Newly added mesh cells and their parents (local ids)
    1689                 :            : // *****************************************************************************
    1690                 :            : {
    1691         [ +  - ]:       1325 :   auto d = Disc();
    1692                 :            : 
    1693                 :       1325 :   const auto& inpoel = std::get< 0 >( m_outmesh.chunk );
    1694         [ +  - ]:       2650 :   auto esup = tk::genEsup( inpoel, 4 );
    1695                 :       1325 :   auto nelem = inpoel.size() / 4;
    1696                 :            : 
    1697                 :            :   // Combine own and communicated contributions and finish averaging of node
    1698                 :            :   // field output in chare boundary nodes
    1699                 :       1325 :   const auto& lid = std::get< 2 >( m_outmesh.chunk );
    1700         [ +  + ]:      51485 :   for (const auto& [g,f] : m_uNodefieldsc) {
    1701 [ -  + ][ -  - ]:      50160 :     Assert( m_uNodefields.nprop() == f.first.size(), "Size mismatch" );
         [ -  - ][ -  - ]
    1702         [ +  - ]:      50160 :     auto p = tk::cref_find( lid, g );
    1703         [ +  + ]:     319474 :     for (std::size_t i=0; i<f.first.size(); ++i) {
    1704         [ +  - ]:     269314 :       m_uNodefields(p,i) += f.first[i];
    1705                 :     269314 :       m_uNodefields(p,i) /= static_cast< tk::real >(
    1706         [ +  - ]:     269314 :                               esup.second[p+1] - esup.second[p] + f.second );
    1707                 :            :     }
    1708                 :            :   }
    1709                 :       1325 :   tk::destroy( m_uNodefieldsc );
    1710         [ +  + ]:      51485 :   for (const auto& [g,f] : m_pNodefieldsc) {
    1711 [ -  + ][ -  - ]:      50160 :     Assert( m_pNodefields.nprop() == f.first.size(), "Size mismatch" );
         [ -  - ][ -  - ]
    1712         [ +  - ]:      50160 :     auto p = tk::cref_find( lid, g );
    1713         [ +  + ]:     100760 :     for (std::size_t i=0; i<f.first.size(); ++i) {
    1714         [ +  - ]:      50600 :       m_pNodefields(p,i) += f.first[i];
    1715                 :      50600 :       m_pNodefields(p,i) /= static_cast< tk::real >(
    1716         [ +  - ]:      50600 :                               esup.second[p+1] - esup.second[p] + f.second );
    1717                 :            :     }
    1718                 :            :   }
    1719                 :       1325 :   tk::destroy( m_pNodefieldsc );
    1720                 :            : 
    1721                 :            :   // Lambda to decide if a node (global id) is on a chare boundary of the field
    1722                 :            :   // output mesh. p - global node id, return true if node is on the chare
    1723                 :            :   // boundary.
    1724                 :     525866 :   auto chbnd = [ this ]( std::size_t p ) {
    1725                 :            :     return
    1726                 :     262933 :       std::any_of( m_outmesh.nodeCommMap.cbegin(), m_outmesh.nodeCommMap.cend(),
    1727         [ +  - ]:     670654 :         [&](const auto& s) { return s.second.find(p) != s.second.cend(); } );
    1728                 :       1325 :   };
    1729                 :            : 
    1730                 :            :   // Finish computing node field output averages in internal nodes
    1731                 :       1325 :   auto npoin = m_outmesh.coord[0].size();
    1732                 :       1325 :   auto& gid = std::get< 1 >( m_outmesh.chunk );
    1733         [ +  + ]:     264258 :   for (std::size_t p=0; p<npoin; ++p) {
    1734 [ +  - ][ +  + ]:     262933 :     if (!chbnd(gid[p])) {
    1735                 :     212773 :       auto n = static_cast< tk::real >( esup.second[p+1] - esup.second[p] );
    1736         [ +  + ]:    1060697 :       for (std::size_t i=0; i<m_uNodefields.nprop(); ++i)
    1737         [ +  - ]:     847924 :         m_uNodefields(p,i) /= n;
    1738         [ +  + ]:     320807 :       for (std::size_t i=0; i<m_pNodefields.nprop(); ++i)
    1739         [ +  - ]:     108034 :         m_pNodefields(p,i) /= n;
    1740                 :            :     }
    1741                 :            :   }
    1742                 :            : 
    1743                 :            :   // Collect field output from numerical solution requested by user
    1744                 :       1325 :   auto elemfields = numericFieldOutput( m_uElemfields, tk::Centering::ELEM,
    1745 [ +  - ][ +  - ]:       2650 :     g_dgpde[Disc()->MeshId()].OutVarFn(), m_pElemfields );
                 [ +  - ]
    1746                 :       1325 :   auto nodefields = numericFieldOutput( m_uNodefields, tk::Centering::NODE,
    1747 [ +  - ][ +  - ]:       2650 :     g_dgpde[Disc()->MeshId()].OutVarFn(), m_pNodefields );
                 [ +  - ]
    1748                 :            : 
    1749                 :            :   // Collect field output from analytical solutions (if exist)
    1750                 :       1325 :   const auto& coord = m_outmesh.coord;
    1751         [ +  - ]:       2650 :   auto geoElem = tk::genGeoElemTet( inpoel, coord );
    1752         [ +  - ]:       1325 :   auto t = Disc()->T();
    1753         [ +  - ]:       1325 :   analyticFieldOutput( g_dgpde[d->MeshId()], tk::Centering::ELEM,
    1754 [ +  - ][ +  - ]:       2650 :     geoElem.extract_comp(1), geoElem.extract_comp(2), geoElem.extract_comp(3),
                 [ +  - ]
    1755                 :            :     t, elemfields );
    1756         [ +  - ]:       1325 :   analyticFieldOutput( g_dgpde[d->MeshId()], tk::Centering::NODE, coord[0],
    1757                 :       1325 :     coord[1], coord[2], t, nodefields );
    1758                 :            : 
    1759                 :            :   // Add adaptive indicator array to element-centered field output
    1760         [ +  + ]:       1325 :   if (g_inputdeck.get< tag::pref, tag::pref >()) {
    1761         [ +  - ]:        804 :     std::vector< tk::real > ndof( begin(m_ndof), end(m_ndof) );
    1762         [ +  - ]:        402 :     ndof.resize( nelem );
    1763         [ +  + ]:     182529 :     for(std::size_t k = 0; k < nelem; k++) {
    1764                 :            :       // Mark the cell with THINC reconstruction as 0 for output
    1765         [ +  + ]:     182127 :       if(m_interface[k] == 1) ndof[k] = 0;
    1766                 :            :     }
    1767         [ +  + ]:      87834 :     for (const auto& [child,parent] : addedTets)
    1768                 :      87432 :       ndof[child] = static_cast< tk::real >( m_ndof[parent] );
    1769         [ +  - ]:        402 :     elemfields.push_back( ndof );
    1770                 :            :   }
    1771                 :            : 
    1772                 :            :   // Add shock detection marker array to element-centered field output
    1773         [ +  - ]:       2650 :   std::vector< tk::real > shockmarker( begin(m_shockmarker), end(m_shockmarker) );
    1774                 :            :   // Here m_shockmarker has a size of m_u.nunk() which is the number of the
    1775                 :            :   // elements within this partition (nelem) plus the ghost partition cells. In
    1776                 :            :   // terms of output purpose, we only need the solution data within this
    1777                 :            :   // partition. Therefore, resizing it to nelem removes the extra partition
    1778                 :            :   // boundary allocations in the shockmarker vector. Since the code assumes that
    1779                 :            :   // the boundary elements are on the top, the resize operation keeps the lower
    1780                 :            :   // portion.
    1781         [ +  - ]:       1325 :   shockmarker.resize( nelem );
    1782         [ +  + ]:     158837 :   for (const auto& [child,parent] : addedTets)
    1783                 :     157512 :     shockmarker[child] = static_cast< tk::real >(m_shockmarker[parent]);
    1784         [ +  - ]:       1325 :   elemfields.push_back( shockmarker );
    1785                 :            : 
    1786                 :            :   // Add rho0*det(g)/rho to make sure it is staying close to 1,
    1787                 :            :   // averaged for all materials
    1788         [ +  - ]:       2650 :   std::vector< tk::real > densityConstr(nelem);
    1789         [ +  - ]:       1325 :   g_dgpde[d->MeshId()].computeDensityConstr(nelem, m_u, densityConstr);
    1790         [ +  + ]:     158837 :   for (const auto& [child,parent] : addedTets)
    1791                 :     157512 :     densityConstr[child] = 0.0;
    1792 [ +  + ][ +  - ]:       1325 :   if (densityConstr.size() > 0) elemfields.push_back( densityConstr );
    1793                 :            : 
    1794                 :            :   // Query fields names requested by user
    1795         [ +  - ]:       2650 :   auto elemfieldnames = numericFieldNames( tk::Centering::ELEM );
    1796         [ +  - ]:       2650 :   auto nodefieldnames = numericFieldNames( tk::Centering::NODE );
    1797                 :            : 
    1798                 :            :   // Collect field output names for analytical solutions
    1799         [ +  - ]:       1325 :   analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::ELEM, elemfieldnames );
    1800         [ +  - ]:       1325 :   analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::NODE, nodefieldnames );
    1801                 :            : 
    1802         [ +  + ]:       1325 :   if (g_inputdeck.get< tag::pref, tag::pref >()) {
    1803 [ +  - ][ +  - ]:        402 :     elemfieldnames.push_back( "NDOF" );
    1804                 :            :   }
    1805                 :            : 
    1806 [ +  - ][ +  - ]:       1325 :   elemfieldnames.push_back( "shock_marker" );
    1807                 :            : 
    1808         [ +  + ]:       1325 :   if (densityConstr.size() > 0)
    1809 [ +  - ][ +  - ]:        179 :     elemfieldnames.push_back( "density_constraint" );
    1810                 :            : 
    1811 [ -  + ][ -  - ]:       1325 :   Assert( elemfieldnames.size() == elemfields.size(), "Size mismatch" );
         [ -  - ][ -  - ]
    1812 [ -  + ][ -  - ]:       1325 :   Assert( nodefieldnames.size() == nodefields.size(), "Size mismatch" );
         [ -  - ][ -  - ]
    1813                 :            : 
    1814                 :            :   // Collect surface output names
    1815         [ +  - ]:       2650 :   auto surfnames = g_dgpde[d->MeshId()].surfNames();
    1816                 :            : 
    1817                 :            :   // Collect surface field solution
    1818         [ +  - ]:       1325 :   const auto& fd = myGhosts()->m_fd;
    1819         [ +  - ]:       1325 :   auto elemsurfs = g_dgpde[d->MeshId()].surfOutput(fd, m_u, m_p);
    1820                 :            : 
    1821                 :            :   // Output chare mesh and fields metadata to file
    1822                 :       1325 :   const auto& triinpoel = m_outmesh.triinpoel;
    1823         [ +  - ]:       3975 :   d->write( inpoel, m_outmesh.coord, m_outmesh.bface, {},
    1824         [ +  - ]:       2650 :             tk::remap( triinpoel, lid ), elemfieldnames, nodefieldnames,
    1825                 :            :             surfnames, {}, elemfields, nodefields, elemsurfs, {}, c );
    1826                 :       1325 : }
    1827                 :            : 
    1828                 :            : void
    1829                 :      10380 : DG::comnodeout( const std::vector< std::size_t >& gid,
    1830                 :            :                 const std::vector< std::size_t >& nesup,
    1831                 :            :                 const std::vector< std::vector< tk::real > >& Lu,
    1832                 :            :                 const std::vector< std::vector< tk::real > >& Lp )
    1833                 :            : // *****************************************************************************
    1834                 :            : //  Receive chare-boundary nodal solution (for field output) contributions from
    1835                 :            : //  neighboring chares
    1836                 :            : //! \param[in] gid Global mesh node IDs at which we receive contributions
    1837                 :            : //! \param[in] nesup Number of elements surrounding points
    1838                 :            : //! \param[in] Lu Partial contributions of solution nodal fields to
    1839                 :            : //!   chare-boundary nodes
    1840                 :            : //! \param[in] Lp Partial contributions of primitive quantity nodal fields to
    1841                 :            : //!   chare-boundary nodes
    1842                 :            : // *****************************************************************************
    1843                 :            : {
    1844 [ -  + ][ -  - ]:      10380 :   Assert( gid.size() == nesup.size(), "Size mismatch" );
         [ -  - ][ -  - ]
    1845 [ -  + ][ -  - ]:      10380 :   Assert(Lu.size() == m_uNodefields.nprop(), "Fields size mismatch");
         [ -  - ][ -  - ]
    1846 [ -  + ][ -  - ]:      10380 :   Assert(Lp.size() == m_pNodefields.nprop(), "Fields size mismatch");
         [ -  - ][ -  - ]
    1847         [ +  + ]:      69172 :   for (std::size_t f=0; f<Lu.size(); ++f)
    1848 [ -  + ][ -  - ]:      58792 :     Assert( gid.size() == Lu[f].size(), "Size mismatch" );
         [ -  - ][ -  - ]
    1849         [ +  + ]:      21008 :   for (std::size_t f=0; f<Lp.size(); ++f)
    1850 [ -  + ][ -  - ]:      10628 :     Assert( gid.size() == Lp[f].size(), "Size mismatch" );
         [ -  - ][ -  - ]
    1851                 :            : 
    1852         [ +  + ]:     107480 :   for (std::size_t i=0; i<gid.size(); ++i) {
    1853                 :      97100 :     auto& nfu = m_uNodefieldsc[ gid[i] ];
    1854                 :      97100 :     nfu.first.resize( Lu.size() );
    1855         [ +  + ]:     641666 :     for (std::size_t f=0; f<Lu.size(); ++f) nfu.first[f] += Lu[f][i];
    1856                 :      97100 :     nfu.second += nesup[i];
    1857                 :      97100 :     auto& nfp = m_pNodefieldsc[ gid[i] ];
    1858                 :      97100 :     nfp.first.resize( Lp.size() );
    1859         [ +  + ]:     206580 :     for (std::size_t f=0; f<Lp.size(); ++f) nfp.first[f] += Lp[f][i];
    1860                 :      97100 :     nfp.second += nesup[i];
    1861                 :            :   }
    1862                 :            : 
    1863                 :            :   // When we have heard from all chares we communicate with, this chare is done
    1864         [ +  + ]:      10380 :   if (++m_nnod == Disc()->NodeCommMap().size()) {
    1865                 :       1164 :     m_nnod = 0;
    1866                 :       1164 :     comnodeout_complete();
    1867                 :            :   }
    1868                 :      10380 : }
    1869                 :            : 
    1870                 :            : void
    1871                 :      38265 : DG::stage()
    1872                 :            : // *****************************************************************************
    1873                 :            : // Evaluate whether to continue with next time step stage
    1874                 :            : // *****************************************************************************
    1875                 :            : {
    1876                 :            :   // Increment Runge-Kutta stage counter
    1877                 :      38265 :   ++m_stage;
    1878                 :            : 
    1879                 :            :   // if not all Runge-Kutta stages complete, continue to next time stage,
    1880                 :            :   // otherwise prepare for nodal field output
    1881         [ +  + ]:      38265 :   if (m_stage < m_nstage)
    1882                 :      25510 :     next();
    1883                 :            :   else
    1884 [ +  - ][ +  - ]:      12755 :     startFieldOutput( CkCallback(CkIndex_DG::step(), thisProxy[thisIndex]) );
                 [ +  - ]
    1885                 :      38265 : }
    1886                 :            : 
    1887                 :            : void
    1888                 :      12084 : DG::evalLB( int nrestart )
    1889                 :            : // *****************************************************************************
    1890                 :            : // Evaluate whether to do load balancing
    1891                 :            : //! \param[in] nrestart Number of times restarted
    1892                 :            : // *****************************************************************************
    1893                 :            : {
    1894                 :      12084 :   auto d = Disc();
    1895                 :            : 
    1896                 :            :   // Detect if just returned from a checkpoint and if so, zero timers
    1897                 :      12084 :   d->restarted( nrestart );
    1898                 :            : 
    1899                 :      12084 :   const auto lbfreq = g_inputdeck.get< tag::cmd, tag::lbfreq >();
    1900                 :      12084 :   const auto nonblocking = g_inputdeck.get< tag::cmd, tag::nonblocking >();
    1901                 :            : 
    1902                 :            :   // Load balancing if user frequency is reached or after the second time-step
    1903 [ -  + ][ -  - ]:      12084 :   if ( (d->It()) % lbfreq == 0 || d->It() == 2 ) {
                 [ +  - ]
    1904                 :            : 
    1905                 :      12084 :     AtSync();
    1906         [ -  + ]:      12084 :     if (nonblocking) next();
    1907                 :            : 
    1908                 :            :   } else {
    1909                 :            : 
    1910                 :          0 :     next();
    1911                 :            : 
    1912                 :            :   }
    1913                 :      12084 : }
    1914                 :            : 
    1915                 :            : void
    1916                 :      12084 : DG::evalRestart()
    1917                 :            : // *****************************************************************************
    1918                 :            : // Evaluate whether to save checkpoint/restart
    1919                 :            : // *****************************************************************************
    1920                 :            : {
    1921                 :      12084 :   auto d = Disc();
    1922                 :            : 
    1923                 :      12084 :   const auto rsfreq = g_inputdeck.get< tag::cmd, tag::rsfreq >();
    1924                 :      12084 :   const auto benchmark = g_inputdeck.get< tag::cmd, tag::benchmark >();
    1925                 :            : 
    1926 [ +  + ][ -  + ]:      12084 :   if (not benchmark and not (d->It() % rsfreq)) {
                 [ -  + ]
    1927                 :            : 
    1928         [ -  - ]:          0 :     std::vector< std::size_t > meshdata{ /* finished = */ 0, d->MeshId() };
    1929         [ -  - ]:          0 :     contribute( meshdata, CkReduction::nop,
    1930 [ -  - ][ -  - ]:          0 :       CkCallback(CkReductionTarget(Transporter,checkpoint), d->Tr()) );
    1931                 :            : 
    1932                 :            :   } else {
    1933                 :            : 
    1934                 :      12084 :     evalLB( /* nrestart = */ -1 );
    1935                 :            : 
    1936                 :            :   }
    1937                 :      12084 : }
    1938                 :            : 
    1939                 :            : void
    1940                 :      12755 : DG::step()
    1941                 :            : // *****************************************************************************
    1942                 :            : // Evaluate wether to continue with next time step
    1943                 :            : // *****************************************************************************
    1944                 :            : {
    1945                 :      12755 :   auto d = Disc();
    1946                 :            : 
    1947                 :            :   // Output time history
    1948 [ +  - ][ +  - ]:      12755 :   if (d->histiter() or d->histtime() or d->histrange()) {
         [ -  + ][ -  + ]
    1949                 :          0 :     std::vector< std::vector< tk::real > > hist;
    1950                 :          0 :     auto h = g_dgpde[d->MeshId()].histOutput( d->Hist(), myGhosts()->m_inpoel,
    1951 [ -  - ][ -  - ]:          0 :       myGhosts()->m_coord, m_u, m_p );
                 [ -  - ]
    1952         [ -  - ]:          0 :     hist.insert( end(hist), begin(h), end(h) );
    1953         [ -  - ]:          0 :     d->history( std::move(hist) );
    1954                 :            :   }
    1955                 :            : 
    1956                 :            :   // Free memory storing output mesh
    1957                 :      12755 :   m_outmesh.destroy();
    1958                 :            : 
    1959                 :            :   // Output one-liner status report to screen
    1960                 :      12755 :   d->status();
    1961                 :            :   // Reset Runge-Kutta stage counter
    1962                 :      12755 :   m_stage = 0;
    1963                 :            : 
    1964                 :      12755 :   const auto term = g_inputdeck.get< tag::term >();
    1965                 :      12755 :   const auto nstep = g_inputdeck.get< tag::nstep >();
    1966                 :      12755 :   const auto eps = std::numeric_limits< tk::real >::epsilon();
    1967                 :            : 
    1968                 :            :   // If neither max iterations nor max time reached, continue, otherwise finish
    1969 [ +  - ][ +  + ]:      12755 :   if (std::fabs(d->T()-term) > eps && d->It() < nstep) {
                 [ +  + ]
    1970                 :            : 
    1971                 :      12084 :     evalRestart();
    1972                 :            :  
    1973                 :            :   } else {
    1974                 :            : 
    1975                 :        671 :     auto meshid = d->MeshId();
    1976         [ +  - ]:        671 :     d->contribute( sizeof(std::size_t), &meshid, CkReduction::nop,
    1977 [ +  - ][ +  - ]:       1342 :                    CkCallback(CkReductionTarget(Transporter,finish), d->Tr()) );
    1978                 :            : 
    1979                 :            :   }
    1980                 :      12755 : }
    1981                 :            : 
    1982                 :            : void
    1983                 :          0 : DG::imex_integrate()
    1984                 :            : // *****************************************************************************
    1985                 :            : //  Perform the Implicit-Explicit Runge-Kutta stage update
    1986                 :            : //
    1987                 :            : //!  \details
    1988                 :            : //!    Performs the Implicit-Explicit Runge-Kutta step. Scheme taken from
    1989                 :            : //!    Cavaglieri, D., & Bewley, T. (2015). Low-storage implicit/explicit
    1990                 :            : //!    Runge–Kutta schemes for the simulation of stiff high-dimensional ODE
    1991                 :            : //!    systems. Journal of Computational Physics, 286, 172-193.
    1992                 :            : //!
    1993                 :            : //!    Scheme given by equations (25a,b):
    1994                 :            : //!
    1995                 :            : //!    u[0] = u[n] + dt * (expl_rkcoef[1,0]*R_ex(u[n])+impl_rkcoef[1,1]*R_im(u[0]))
    1996                 :            : //!
    1997                 :            : //!    u[1] = u[n] + dt * (expl_rkcoef[2,1]*R_ex(u[0])+impl_rkcoef[2,1]*R_im(u[0])
    1998                 :            : //!                                                   +impl_rkcoef[2,2]*R_im(u[1]))
    1999                 :            : //!
    2000                 :            : //!    u[n+1] = u[n] + dt * (expl_rkcoef[3,1]*R_ex(u[0])+impl_rkcoef[3,1]*R_im(u[0])
    2001                 :            : //!                          expl_rkcoef[3,2]*R_ex(u[1])+impl_rkcoef[3,2]*R_im(u[1]))
    2002                 :            : //!
    2003                 :            : //!    In order to solve the first two equations we need to solve a series of
    2004                 :            : //!    systems of non-linear equations:
    2005                 :            : //!
    2006                 :            : //!    F1(u[0]) = B1 + R1(u[0]) = 0, and
    2007                 :            : //!    F2(u[1]) = B2 + R2(u[1]) = 0,
    2008                 :            : //!
    2009                 :            : //!    where
    2010                 :            : //!
    2011                 :            : //!    B1 = u[n] + dt * expl_rkcoef[1,0]*R_ex(u[n]),
    2012                 :            : //!    R1 = dt * impl_rkcoef[1,1]*R_im(u[0]) - u([0]),
    2013                 :            : //!    B2 = u[n] + dt * (expl_rkcoef[2,1]*R_ex(u[0])+impl_rkcoef[2,1]*R_im(u[0])),
    2014                 :            : //!    R2 = dt * impl_rkcoef[2,2]*R_im(u[1]) - u([1]).
    2015                 :            : //!
    2016                 :            : //!    In order to solve the non-linear system F(U) = 0, we employ:
    2017                 :            : //!    First, Broyden's method.
    2018                 :            : //!    If that fails, Newton's method (with FD approximation for jacobian).
    2019                 :            : //!
    2020                 :            : //!    Broyden's method:
    2021                 :            : //!    ----------------
    2022                 :            : //!
    2023                 :            : //!    Taken from https://en.wikipedia.org/wiki/Broyden%27s_method.
    2024                 :            : //!    The method consists in obtaining an approximation for the inverse of the
    2025                 :            : //!    Jacobian H = J^(-1) and advancing in a quasi-newton step:
    2026                 :            : //!
    2027                 :            : //!    U[k+1] = U[k] - H[k]*F(U[k]),
    2028                 :            : //!
    2029                 :            : //!    until F(U) is close enough to zero.
    2030                 :            : //!
    2031                 :            : //!    The approximation H[k] is improved at every iteration following
    2032                 :            : //!
    2033                 :            : //!    H[k] = H[k-1] + (DU[k]-H[k-1]*DF[k])/(DU[k]^T*H[k-1]*DF[k]) * DU[k]^T*H[k-1],
    2034                 :            : //!
    2035                 :            : //!    where DU[k] = U[k] - U[k-1] and DF[k] = F(U[k]) - F(U[k-1)).
    2036                 :            : //!
    2037                 :            : //!    Newton's method:
    2038                 :            : //!    ----------------
    2039                 :            : //!
    2040                 :            : //!    Taken from https://en.wikipedia.org/wiki/Newton%27s_method
    2041                 :            : //!    The method consists in inverting the jacobian
    2042                 :            : //!    Jacobian J and advancing in a Newton step:
    2043                 :            : //!
    2044                 :            : //!    U[k+1] = U[k] - J^(-1)[k]*F(U[k]),
    2045                 :            : //!
    2046                 :            : //!    until F(U) is close enough to zero.
    2047                 :            : //!
    2048                 :            : //!
    2049                 :            : //!    This function performs the following main algorithmic steps:
    2050                 :            : //!    - If stage == 0 or stage == 1:
    2051                 :            : //!      - Take Initial value:
    2052                 :            : //!        U[0] = U[n] + dt * expl_rkcoef[1,0]*R_ex(U[n]) (for stage 0)
    2053                 :            : //!        U[1] = U[n] + dt * (expl_rkcoef[2,1]*R_ex(U[0])
    2054                 :            : //!                           +impl_rkcoef[2,1]*R_im(U[0])) (for stage 1)
    2055                 :            : //!      - Loop over the Elements (e++)
    2056                 :            : //!        - Broyden steps:
    2057                 :            : //!        - Initialize Jacobian inverse approximation using FD
    2058                 :            : //!        - Compute implicit right-hand-side (F_im) with current U
    2059                 :            : //!        - Iterate for the solution (iter++)
    2060                 :            : //!          - Perform line search prior to solution update
    2061                 :            : //!          - Compute new solution U[k+1] = U[k] - H[k]*F(U[k])
    2062                 :            : //!          - Compute implicit right-hand-side (F_im) with current U
    2063                 :            : //!          - Compute DU and DF
    2064                 :            : //!          - Update inverse Jacobian approximation by:
    2065                 :            : //!            - Compute V1 = H[k-1]*DF[k] and V2 = DU[k]^T*H[k-1]
    2066                 :            : //!            - Compute d = DU[k]^T*V1 and V3 = DU[k]-V1
    2067                 :            : //!            - Compute V4 = V3/d
    2068                 :            : //!            - Update H[k] = H[k-1] + V4*V2
    2069                 :            : //!          - Save old U and F
    2070                 :            : //!          - Compute absolute and relative errors
    2071                 :            : //!          - Break iterations if error < tol or iter == max_iter
    2072                 :            : //!        - Newton steps:
    2073                 :            : //!          - Initialize Jacobian using FD approximation.
    2074                 :            : //!          - Compute implicit right-hand-side (F_im) with current U
    2075                 :            : //!          - Iterate for the solution (iter++)
    2076                 :            : //!          - Perform line search prior to solution update
    2077                 :            : //!          - Compute new solution U[k+1] = U[k] - J^(-1)[k]*F(U[k])
    2078                 :            : //!          - Compute implicit right-hand-side (F_im) with current U
    2079                 :            : //!          - Compute DU and DF
    2080                 :            : //!          - Save old U and F
    2081                 :            : //!          - Compute absolute and relative errors
    2082                 :            : //!          - Break iterations if error < tol or iter == max_iter
    2083                 :            : //!       - Update explicit equations using only the explicit terms.
    2084                 :            : //!    - Else if stage == 2:
    2085                 :            : //!       - Update explicit equations using only the explicit terms.
    2086                 :            : //!       - Update implicit equations using:
    2087                 :            : //!       u[n+1] = u[n]+dt*(expl_rkcoef[3,1]*R_ex(u[0])+impl_rkcoef[3,1]*R_im(u[0])
    2088                 :            : //!                         expl_rkcoef[3,2]*R_ex(u[1])+impl_rkcoef[3,2]*R_im(u[1]))
    2089                 :            : // *****************************************************************************
    2090                 :            : {
    2091                 :          0 :   auto d = Disc();
    2092                 :          0 :   const auto rdof = g_inputdeck.get< tag::rdof >();
    2093                 :          0 :   const auto ndof = g_inputdeck.get< tag::ndof >();
    2094         [ -  - ]:          0 :   if (m_stage < m_nstage-1) {
    2095                 :            :     // Save previous stiff_rhs
    2096                 :          0 :     m_stiffrhsprev = m_stiffrhs;
    2097                 :            : 
    2098                 :            :     // Compute the imex update
    2099                 :            : 
    2100                 :            :     // Integrate explicitly on the imex equations
    2101                 :            :     // (To use as initial values)
    2102         [ -  - ]:          0 :     for (std::size_t e=0; e<myGhosts()->m_nunk; ++e) {
    2103                 :          0 :       auto vole = myGhosts()->m_geoElem(e,0);
    2104         [ -  - ]:          0 :       for (std::size_t c=0; c<m_nstiffeq; ++c)
    2105                 :            :       {
    2106         [ -  - ]:          0 :         for (std::size_t k=0; k<m_numEqDof[c]; ++k)
    2107                 :            :         {
    2108                 :          0 :           auto rmark = m_stiffEqIdx[c]*rdof+k;
    2109                 :          0 :           auto mark = m_stiffEqIdx[c]*ndof+k;
    2110                 :          0 :           m_u(e, rmark) =  m_un(e, rmark) + d->Dt() * (
    2111                 :          0 :             expl_rkcoef[0][m_stage] * m_rhsprev(e, mark)/(vole*mass_dubiner[k])
    2112                 :          0 :             + expl_rkcoef[1][m_stage] * m_rhs(e, mark)/(vole*mass_dubiner[k])
    2113                 :          0 :             + impl_rkcoef[0][m_stage]
    2114                 :          0 :             * m_stiffrhsprev(e,c*ndof+k)/(vole*mass_dubiner[k]) );
    2115         [ -  - ]:          0 :           if(fabs(m_u(e, rmark)) < 1e-16)
    2116                 :          0 :             m_u(e, rmark) = 0;
    2117                 :            :         }
    2118                 :            :       }
    2119                 :            :     }
    2120                 :            : 
    2121                 :            :     // Solve for implicit-explicit equations
    2122                 :          0 :     const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
    2123         [ -  - ]:          0 :     for (std::size_t e=0; e<nelem; ++e)
    2124                 :            :     {
    2125                 :            : 
    2126                 :            :       // Non-linear solver solves for x.
    2127                 :            :       // Copy the relevant variables from the state array into x.
    2128         [ -  - ]:          0 :       std::vector< tk::real > x(m_nstiffeq*ndof, 0.0);
    2129         [ -  - ]:          0 :       for (size_t ieq=0; ieq<m_nstiffeq; ++ieq)
    2130         [ -  - ]:          0 :         for (size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
    2131                 :            :         {
    2132                 :          0 :           auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
    2133         [ -  - ]:          0 :           x[ieq*ndof+idof] = m_u(e, stiffrmark);
    2134                 :            :         }
    2135                 :            : 
    2136                 :            :       // Save all the values of m_u at stiffEqIdx as x_star,
    2137                 :            :       // They will serve to balance the energy exchange
    2138                 :            :       // from the implicit step
    2139         [ -  - ]:          0 :       auto x_star = x;
    2140                 :            : 
    2141                 :            :       // Solve nonlinear system, first try broyden
    2142                 :          0 :       bool solver_failed = false;
    2143 [ -  - ][ -  - ]:          0 :       x = DG::nonlinear_broyden(e, x, solver_failed);
    2144                 :            : 
    2145                 :            :       // If solver_failed, do newton
    2146         [ -  - ]:          0 :       if (solver_failed) {
    2147                 :          0 :         solver_failed = false;
    2148 [ -  - ][ -  - ]:          0 :         x = DG::nonlinear_newton(e, x, solver_failed);
    2149                 :            :       }
    2150                 :            : 
    2151                 :            :       // If newton failed, crash
    2152         [ -  - ]:          0 :       if (solver_failed)
    2153 [ -  - ][ -  - ]:          0 :         Throw("At element " + std::to_string(e) +
         [ -  - ][ -  - ]
                 [ -  - ]
    2154                 :            :               " nonlinear solvers was not able to converge");
    2155                 :            : 
    2156                 :            :       // Balance energy
    2157 [ -  - ][ -  - ]:          0 :       g_dgpde[d->MeshId()].balance_plastic_energy(e, x_star, x, m_un);
                 [ -  - ]
    2158                 :            : 
    2159                 :            :       // Update the state u with the converged vector x.
    2160         [ -  - ]:          0 :       for (size_t ieq=0; ieq<m_nstiffeq; ++ieq)
    2161         [ -  - ]:          0 :         for (size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
    2162                 :            :         {
    2163                 :          0 :           auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
    2164         [ -  - ]:          0 :           m_u(e, stiffrmark) = x[ieq*ndof+idof];
    2165                 :            :         }
    2166                 :            : 
    2167                 :            :     }
    2168                 :            : 
    2169                 :            :     // Then, integrate explicitly on the remaining equations
    2170         [ -  - ]:          0 :     for (std::size_t e=0; e<nelem; ++e) {
    2171                 :          0 :       auto vole = myGhosts()->m_geoElem(e,0);
    2172         [ -  - ]:          0 :       for (std::size_t c=0; c<m_nnonstiffeq; ++c)
    2173                 :            :       {
    2174         [ -  - ]:          0 :         for (std::size_t k=0; k<m_numEqDof[c]; ++k)
    2175                 :            :         {
    2176                 :          0 :           auto rmark = m_nonStiffEqIdx[c]*rdof+k;
    2177                 :          0 :           auto mark = m_nonStiffEqIdx[c]*ndof+k;
    2178                 :          0 :           m_u(e, rmark) =  m_un(e, rmark) + d->Dt() * (
    2179                 :          0 :             expl_rkcoef[0][m_stage] * m_rhsprev(e, mark)/(vole*mass_dubiner[k])
    2180                 :          0 :             + expl_rkcoef[1][m_stage] * m_rhs(e, mark)/(vole*mass_dubiner[k]));
    2181         [ -  - ]:          0 :           if(fabs(m_u(e, rmark)) < 1e-16)
    2182                 :          0 :             m_u(e, rmark) = 0;
    2183                 :            :         }
    2184                 :            :       }
    2185                 :            :     }
    2186                 :            :   }
    2187                 :            :   else {
    2188                 :            :     // For last stage just use all previously computed stages
    2189                 :          0 :     const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
    2190         [ -  - ]:          0 :     for (std::size_t e=0; e<nelem; ++e)
    2191                 :            :     {
    2192                 :          0 :       auto vole = myGhosts()->m_geoElem(e,0);
    2193                 :            :       // First integrate explicitly on nonstiff equations
    2194         [ -  - ]:          0 :       for (std::size_t c=0; c<m_nnonstiffeq; ++c)
    2195                 :            :       {
    2196         [ -  - ]:          0 :         for (std::size_t k=0; k<m_numEqDof[c]; ++k)
    2197                 :            :         {
    2198                 :          0 :           auto rmark = m_nonStiffEqIdx[c]*rdof+k;
    2199                 :          0 :           auto mark = m_nonStiffEqIdx[c]*ndof+k;
    2200                 :          0 :           m_u(e, rmark) =  m_un(e, rmark) + d->Dt() * (
    2201                 :          0 :             expl_rkcoef[0][m_stage] * m_rhsprev(e, mark)/(vole*mass_dubiner[k])
    2202                 :          0 :             + expl_rkcoef[1][m_stage] * m_rhs(e, mark)/(vole*mass_dubiner[k]));
    2203         [ -  - ]:          0 :           if(fabs(m_u(e, rmark)) < 1e-16)
    2204                 :          0 :             m_u(e, rmark) = 0;
    2205                 :            :         }
    2206                 :            :       }
    2207                 :            :       // Then, integrate the imex-equations
    2208         [ -  - ]:          0 :       for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
    2209         [ -  - ]:          0 :         for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
    2210                 :            :         {
    2211                 :          0 :           auto rmark = m_stiffEqIdx[ieq]*rdof+idof;
    2212                 :          0 :           auto mark = m_stiffEqIdx[ieq]*ndof+idof;
    2213                 :          0 :           m_u(e, rmark) = m_un(e, rmark)
    2214                 :          0 :             + d->Dt() * (expl_rkcoef[0][m_stage]
    2215                 :          0 :                          * m_rhsprev(e,mark)/(vole*mass_dubiner[idof])
    2216                 :          0 :                          + expl_rkcoef[1][m_stage]
    2217                 :          0 :                          * m_rhs(e,mark)/(vole*mass_dubiner[idof])
    2218                 :          0 :                          + impl_rkcoef[0][m_stage]
    2219                 :          0 :                          * m_stiffrhsprev(e,ieq*ndof+idof)/(vole*mass_dubiner[idof])
    2220                 :          0 :                          + impl_rkcoef[1][m_stage]
    2221                 :          0 :                          * m_stiffrhs(e,ieq*ndof+idof)/(vole*mass_dubiner[idof]) );
    2222         [ -  - ]:          0 :           if(fabs(m_u(e, rmark)) < 1e-16)
    2223                 :          0 :             m_u(e, rmark) = 0;
    2224                 :            :         }
    2225                 :            :     }
    2226                 :            :   }
    2227                 :          0 : }
    2228                 :            : 
    2229                 :          0 : std::vector< tk::real > DG::nonlinear_func(std::size_t e,
    2230                 :            :                                            std::vector< tk::real > x)
    2231                 :            : // *****************************************************************************
    2232                 :            : // Evaluate the stiff RHS and stiff equations f = b - A(x)
    2233                 :            : //! \param[in] e Element number
    2234                 :            : //! \param[in,out] x Array of unknowns to solve for
    2235                 :            : //! \details
    2236                 :            : //!   Defines the F(x) function that the non-linear solvers
    2237                 :            : //!   look to minimize. Deals with properly calling the stiff
    2238                 :            : //!   RHS functions.
    2239                 :            : // *****************************************************************************
    2240                 :            : {
    2241         [ -  - ]:          0 :   auto d = Disc();
    2242                 :          0 :   const auto rdof = g_inputdeck.get< tag::rdof >();
    2243                 :          0 :   const auto ndof = g_inputdeck.get< tag::ndof >();
    2244                 :          0 :   std::size_t n = x.size();
    2245                 :            : 
    2246                 :            :   // m_u <- x
    2247         [ -  - ]:          0 :   for (size_t ieq=0; ieq<m_nstiffeq; ++ieq)
    2248         [ -  - ]:          0 :     for (size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
    2249                 :            :     {
    2250                 :          0 :       auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
    2251         [ -  - ]:          0 :       m_u(e, stiffrmark) = x[ieq*ndof+idof];
    2252                 :            :     }
    2253                 :            : 
    2254 [ -  - ][ -  - ]:          0 :   auto vole = myGhosts()->m_geoElem(e,0);
    2255                 :            : 
    2256                 :            :   // Compute explicit terms (Should be computed once)
    2257         [ -  - ]:          0 :   std::vector< tk::real > expl_terms(n, 0.0);
    2258         [ -  - ]:          0 :   for (size_t ieq=0; ieq<m_nstiffeq; ++ieq)
    2259         [ -  - ]:          0 :     for (size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
    2260                 :            :     {
    2261                 :          0 :       auto stiffmark = m_stiffEqIdx[ieq]*ndof+idof;
    2262                 :          0 :       auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
    2263         [ -  - ]:          0 :       expl_terms[ieq*ndof+idof] = m_un(e, stiffrmark)
    2264                 :          0 :         + d->Dt() * ( expl_rkcoef[0][m_stage]
    2265         [ -  - ]:          0 :         * m_rhsprev(e,stiffmark)/(vole*mass_dubiner[idof])
    2266                 :          0 :         + expl_rkcoef[1][m_stage]
    2267         [ -  - ]:          0 :         * m_rhs(e,stiffmark)/(vole*mass_dubiner[idof])
    2268                 :          0 :         + impl_rkcoef[0][m_stage]
    2269         [ -  - ]:          0 :         * m_stiffrhsprev(e,ieq*ndof+idof)/(vole*mass_dubiner[idof]) );
    2270                 :            :     }
    2271                 :            : 
    2272                 :            :   // Compute stiff_rhs
    2273         [ -  - ]:          0 :   g_dgpde[d->MeshId()].stiff_rhs( e, myGhosts()->m_geoElem,
    2274 [ -  - ][ -  - ]:          0 :     myGhosts()->m_inpoel, myGhosts()->m_coord,
    2275         [ -  - ]:          0 :     m_u, m_p, m_ndof, m_stiffrhs );
    2276                 :            : 
    2277                 :            :   // Store f
    2278         [ -  - ]:          0 :   std::vector< tk::real > f(n, 0.0);
    2279         [ -  - ]:          0 :   for (std::size_t ieq=0; ieq<m_nstiffeq; ++ieq)
    2280         [ -  - ]:          0 :     for (std::size_t idof=0; idof<m_numEqDof[ieq]; ++idof)
    2281                 :            :     {
    2282                 :          0 :       auto stiffrmark = m_stiffEqIdx[ieq]*rdof+idof;
    2283                 :          0 :       f[ieq*ndof+idof] = expl_terms[ieq*ndof+idof]
    2284                 :          0 :         + d->Dt() * impl_rkcoef[1][m_stage]
    2285         [ -  - ]:          0 :         * m_stiffrhs(e,ieq*ndof+idof)/(vole*mass_dubiner[idof])
    2286         [ -  - ]:          0 :         - m_u(e, stiffrmark);
    2287                 :            :     }
    2288                 :            : 
    2289                 :          0 :   return f;
    2290                 :            : }
    2291                 :            : 
    2292                 :          0 : std::vector< tk::real > DG::nonlinear_broyden(std::size_t e,
    2293                 :            :                                               std::vector< tk::real > x,
    2294                 :            :                                               bool solver_failed )
    2295                 :            : // *****************************************************************************
    2296                 :            : // Performs Broyden's method to solve a non-linear system on
    2297                 :            : // element e.
    2298                 :            : //! \param[in] e Element number
    2299                 :            : //! \param[in,out] x Array of unknowns to solve for
    2300                 :            : //! \param[out] solver_failed Returns true if solver did not converge
    2301                 :            : //! \details
    2302                 :            : //!    Taken from https://en.wikipedia.org/wiki/Broyden%27s_method.
    2303                 :            : //!    The method consists in obtaining an approximation for the inverse of the
    2304                 :            : //!    Jacobian H = J^(-1) and advancing in a quasi-newton step:
    2305                 :            : //!
    2306                 :            : //!    U[k+1] = U[k] - H[k]*F(U[k]),
    2307                 :            : //!
    2308                 :            : //!    until F(U) is close enough to zero.
    2309                 :            : // *****************************************************************************
    2310                 :            : {
    2311                 :            :   // Broyden's method
    2312                 :            :   // Control parameters
    2313                 :          0 :   std::size_t max_iter = g_inputdeck.get< tag::imex_maxiter >();
    2314                 :          0 :   tk::real rel_tol = g_inputdeck.get< tag::imex_reltol >();
    2315                 :          0 :   tk::real abs_tol = g_inputdeck.get< tag::imex_abstol >();
    2316                 :          0 :   tk::real rel_err = rel_tol+1;
    2317                 :          0 :   tk::real abs_err = abs_tol+1;
    2318                 :          0 :   std::size_t n = x.size();
    2319                 :            : 
    2320                 :            :   // Compute f with initial guess
    2321 [ -  - ][ -  - ]:          0 :   std::vector< tk::real > f = DG::nonlinear_func(e, x);
    2322                 :            : 
    2323                 :            :   // Initialize x_old and f_old
    2324 [ -  - ][ -  - ]:          0 :   std::vector< tk::real > x_old(n, 0.0), f_old(n, 0.0);
    2325         [ -  - ]:          0 :   for (std::size_t i=0; i<n; ++i)
    2326                 :            :   {
    2327                 :          0 :     x_old[i] = x[i];
    2328                 :          0 :     f_old[i] = f[i];
    2329                 :            :   }
    2330                 :            : 
    2331                 :            :   // Initialize delta_x and delta_f
    2332 [ -  - ][ -  - ]:          0 :   std::vector< tk::real > delta_x(n, 0.0), delta_f(n, 0.0);
    2333                 :            : 
    2334                 :            :   // Store the norm of f initially, for relative error measure
    2335                 :          0 :   tk::real err0 = 0.0;
    2336         [ -  - ]:          0 :   for (std::size_t i=0; i<n; ++i)
    2337                 :          0 :     err0 += f[i]*f[i];
    2338                 :          0 :   err0 = std::sqrt(err0);
    2339                 :            : 
    2340                 :            :   // Iterate for the solution if err0 > 0
    2341                 :          0 :   solver_failed = false;
    2342                 :          0 :   tk::real alpha_jacob = 1.0;
    2343         [ -  - ]:          0 :   if (err0 > abs_tol) {
    2344                 :            : 
    2345                 :            :     // Evaluate finite difference based jacobian
    2346         [ -  - ]:          0 :     std::vector< double > jacob(n*n);
    2347                 :          0 :     tk::real dx = 0.0;
    2348         [ -  - ]:          0 :     for (std::size_t i=0; i<n; ++i)
    2349         [ -  - ]:          0 :       for (std::size_t j=0; j<n; ++j)
    2350                 :            :       {
    2351                 :            :         // Set dx in the order 1% of the unknown
    2352                 :          0 :         dx = std::max(std::abs(0.1*x[j]),1.0e-06);
    2353                 :            :         // Derivative of f[i] with respect to x[j]
    2354         [ -  - ]:          0 :         auto x_perturb = x;
    2355                 :          0 :         x_perturb[j] += dx;
    2356 [ -  - ][ -  - ]:          0 :         auto f_perturb = DG::nonlinear_func(e, x_perturb);
    2357                 :          0 :         jacob[i*n+j] = (f_perturb[i]-f[i])/dx;
    2358                 :            :       }
    2359                 :            : 
    2360                 :            :     // Initialize Jacobian to be the inverse of this jacobian
    2361                 :          0 :     lapack_int ln = static_cast< lapack_int >(n);
    2362         [ -  - ]:          0 :     std::vector< lapack_int > ipiv(n);
    2363                 :            : 
    2364                 :            :     #ifndef NDEBUG
    2365                 :            :     lapack_int ierr =
    2366                 :            :     #endif
    2367         [ -  - ]:          0 :       LAPACKE_dgetrf(LAPACK_ROW_MAJOR, ln, ln, jacob.data(), ln, ipiv.data());
    2368 [ -  - ][ -  - ]:          0 :     Assert(ierr==0, "Lapack error in LU factorization of FD Jacobian");
         [ -  - ][ -  - ]
    2369                 :            : 
    2370                 :            :     #ifndef NDEBUG
    2371                 :            :     lapack_int jerr =
    2372                 :            :     #endif
    2373         [ -  - ]:          0 :       LAPACKE_dgetri(LAPACK_ROW_MAJOR, ln, jacob.data(), ln, ipiv.data());
    2374 [ -  - ][ -  - ]:          0 :     Assert(jerr==0, "Lapack error in inverting FD Jacobian");
         [ -  - ][ -  - ]
    2375                 :            : 
    2376                 :            :     std::vector< std::vector< tk::real > >
    2377 [ -  - ][ -  - ]:          0 :       approx_jacob(n, std::vector< tk::real >(n, 0.0));
    2378         [ -  - ]:          0 :     for (std::size_t i=0; i<n; ++i)
    2379         [ -  - ]:          0 :       for (std::size_t j=0; j<n; ++j)
    2380                 :          0 :         approx_jacob[i][j] = jacob[i*n+j];
    2381                 :            : 
    2382         [ -  - ]:          0 :     for (size_t iter=0; iter<max_iter; ++iter)
    2383                 :            :     {
    2384                 :            : 
    2385                 :            :       // Scale inverse of jacobian if things are not going well
    2386         [ -  - ]:          0 :       for (std::size_t i=0; i<n; ++i)
    2387         [ -  - ]:          0 :         for (std::size_t j=0; j<n; ++j)
    2388                 :          0 :           approx_jacob[i][j] *= alpha_jacob;
    2389                 :            : 
    2390                 :            :       // Compute new solution
    2391         [ -  - ]:          0 :       std::vector < tk::real > delta(n, 0.0);
    2392         [ -  - ]:          0 :       for (std::size_t i=0; i<n; ++i)
    2393                 :            :       {
    2394                 :          0 :         delta[i] = 0.0;
    2395         [ -  - ]:          0 :         for (std::size_t j=0; j<n; ++j)
    2396                 :          0 :           delta[i] -= approx_jacob[i][j] * f[j];
    2397                 :            :       }
    2398                 :            : 
    2399                 :            :       // Update x using line search
    2400                 :          0 :       bool ls_failed = false;
    2401                 :          0 :       tk::real alpha_ls = 1.0E+00;
    2402                 :          0 :       std::size_t nline = 25;
    2403         [ -  - ]:          0 :       auto xtest = x;
    2404         [ -  - ]:          0 :       for (std::size_t iline = 0; iline<nline; ++iline)
    2405                 :            :       {
    2406                 :            :         // Evaluate xtest
    2407         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2408                 :          0 :           xtest[i] = x[i] + alpha_ls*delta[i];
    2409                 :            : 
    2410                 :            :         // Compute new f(x)
    2411 [ -  - ][ -  - ]:          0 :         f = DG::nonlinear_func(e, xtest);
    2412                 :            : 
    2413                 :          0 :         tk::real err = 0.0;
    2414         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2415                 :          0 :           err += f[i]*f[i];
    2416                 :          0 :         abs_err = std::sqrt(err);
    2417                 :            : 
    2418                 :            :         // If 1. The error went up
    2419                 :            :         // or 2. The function f flipped in sign
    2420                 :            :         // Reduce the factor alpha_ls
    2421                 :          0 :         bool flipped_sign = false;
    2422         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2423         [ -  - ]:          0 :           if (f_old[i]*f[i] < 0.0) {
    2424                 :          0 :             flipped_sign = true;
    2425                 :          0 :             break;
    2426                 :            :           }
    2427                 :            : 
    2428         [ -  - ]:          0 :         if (!flipped_sign)
    2429                 :            :         {
    2430                 :          0 :           break;
    2431                 :            :         }
    2432                 :            :         else
    2433                 :            :         {
    2434                 :          0 :           alpha_ls *= 0.5;
    2435                 :            :         }
    2436         [ -  - ]:          0 :         if (iline == nline-1) {
    2437                 :            :           // Try again by reducing the jacobian,
    2438                 :            :           // but only a few times, otherwise give up
    2439                 :          0 :           alpha_jacob *= 0.5;
    2440         [ -  - ]:          0 :           if (alpha_jacob < 1.0E-03)
    2441                 :          0 :             solver_failed = true;
    2442                 :            :           else
    2443                 :          0 :             ls_failed = true;
    2444                 :            :         }
    2445                 :            :       }
    2446                 :            : 
    2447         [ -  - ]:          0 :       if (solver_failed) {
    2448                 :          0 :         break;
    2449                 :            :       }
    2450                 :            : 
    2451         [ -  - ]:          0 :       if (!ls_failed) {
    2452                 :            :         // Save x
    2453         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2454                 :          0 :           x[i] = xtest[i];
    2455                 :            : 
    2456                 :            :         // Compute delta_x and delta_f
    2457         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2458                 :            :         {
    2459                 :          0 :           delta_x[i] = x[i] - x_old[i];
    2460                 :          0 :           delta_f[i] = f[i] - f_old[i];
    2461                 :            :         }
    2462                 :            : 
    2463                 :            :         // Update inverse Jacobian approximation
    2464                 :            : 
    2465                 :            :         // 1. Compute approx_jacob*delta_f and delta_x*jacob_approx
    2466                 :            :         tk::real sum1, sum2;
    2467 [ -  - ][ -  - ]:          0 :         std::vector< tk::real > auxvec1(n, 0.0), auxvec2(n, 0.0);
    2468         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2469                 :            :         {
    2470                 :          0 :           sum1 = 0.0;
    2471                 :          0 :           sum2 = 0.0;
    2472         [ -  - ]:          0 :           for (std::size_t j=0; j<n; ++j)
    2473                 :            :           {
    2474                 :          0 :             sum1 += approx_jacob[i][j]*delta_f[j];
    2475                 :          0 :             sum2 += delta_x[j]*approx_jacob[j][i];
    2476                 :            :           }
    2477                 :          0 :           auxvec1[i] = sum1;
    2478                 :          0 :           auxvec2[i] = sum2;
    2479                 :            :         }
    2480                 :            : 
    2481                 :            :         // 2. Compute delta_x*approx_jacob*delta_f
    2482                 :            :         // and delta_x-approx_jacob*delta_f
    2483                 :          0 :         tk::real denom = 0.0;
    2484         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2485                 :            :         {
    2486                 :          0 :           denom += delta_x[i]*auxvec1[i];
    2487                 :          0 :           auxvec1[i] = delta_x[i]-auxvec1[i];
    2488                 :            :         }
    2489                 :            : 
    2490                 :            :         // 3. Divide delta_u+approx_jacob*delta_f
    2491                 :            :         // by delta_x*(approx_jacob*delta_f)
    2492         [ -  - ]:          0 :         if (std::abs(denom) < 1.0e-18)
    2493                 :            :         {
    2494         [ -  - ]:          0 :           if (denom < 0.0)
    2495                 :            :           {
    2496         [ -  - ]:          0 :             for (std::size_t i=0; i<n; ++i)
    2497                 :          0 :               auxvec1[i] /= -1.0e-18;
    2498                 :            :           }
    2499                 :            :           else
    2500                 :            :           {
    2501         [ -  - ]:          0 :             for (std::size_t i=0; i<n; ++i)
    2502                 :          0 :               auxvec1[i] /= 1.0e-18;
    2503                 :            :           }
    2504                 :            :         }
    2505                 :            :         else
    2506                 :            :         {
    2507         [ -  - ]:          0 :           for (std::size_t i=0; i<n; ++i)
    2508                 :          0 :             auxvec1[i] /= denom;
    2509                 :            :         }
    2510                 :            : 
    2511                 :            :         // 4. Perform outter product between the two arrays and
    2512                 :            :         // add that quantity to the new jacobian approximation
    2513         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2514         [ -  - ]:          0 :           for (std::size_t j=0; j<n; ++j)
    2515                 :          0 :             approx_jacob[i][j] += auxvec1[i] * auxvec2[j];
    2516                 :            : 
    2517                 :            :         // Compute a measure of error, use norm of f
    2518                 :          0 :         tk::real err = 0.0;
    2519         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2520                 :          0 :           err += f[i]*f[i];
    2521                 :          0 :         abs_err = std::sqrt(err);
    2522                 :          0 :         rel_err = abs_err/err0;
    2523                 :            : 
    2524                 :            :         // Save solution and f
    2525         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2526                 :            :         {
    2527                 :          0 :           x_old[i] = x[i];
    2528                 :          0 :           f_old[i] = f[i];
    2529                 :            :         }
    2530                 :            : 
    2531                 :            :         // check if error condition is met and loop back
    2532 [ -  - ][ -  - ]:          0 :         if (rel_err < rel_tol || abs_err < abs_tol)
    2533                 :            :           break;
    2534                 :            : 
    2535                 :            :         // If we did not converge, print a message and keep going
    2536         [ -  - ]:          0 :         if (iter == max_iter-1)
    2537                 :            :         {
    2538                 :          0 :           solver_failed = true;
    2539                 :            :         }
    2540                 :            :       }
    2541                 :            :     }
    2542                 :            :   }
    2543                 :            : 
    2544                 :          0 :   return x;
    2545                 :            : }
    2546                 :            : 
    2547                 :          0 : std::vector< tk::real > DG::nonlinear_newton(std::size_t e,
    2548                 :            :                                              std::vector< tk::real > x,
    2549                 :            :                                              bool solver_failed )
    2550                 :            : // *****************************************************************************
    2551                 :            : // Performs Newton's method to solve a non-linear system on
    2552                 :            : // element e.
    2553                 :            : //! \param[in] e Element number
    2554                 :            : //! \param[in,out] x Array of unknowns to solve for
    2555                 :            : //! \param[out] solver_failed Returns true if solver did not converge
    2556                 :            : //! \details
    2557                 :            : //!    Taken from https://en.wikipedia.org/wiki/Newton%27s_method
    2558                 :            : //!    The method consists in inverting the jacobian
    2559                 :            : //!    Jacobian J and advancing in a Newton step:
    2560                 :            : //!
    2561                 :            : //!    U[k+1] = U[k] - J^(-1)[k]*F(U[k]),
    2562                 :            : //!
    2563                 :            : //!    until F(U) is close enough to zero.
    2564                 :            : // *****************************************************************************
    2565                 :            : {
    2566                 :            :   // Newton's method
    2567                 :            :   // Control parameters
    2568                 :          0 :   std::size_t max_iter = g_inputdeck.get< tag::imex_maxiter >();
    2569                 :          0 :   tk::real rel_tol = g_inputdeck.get< tag::imex_reltol >();
    2570                 :          0 :   tk::real abs_tol = g_inputdeck.get< tag::imex_abstol >();
    2571                 :          0 :   tk::real rel_err = rel_tol+1;
    2572                 :          0 :   tk::real abs_err = abs_tol+1;
    2573                 :          0 :   std::size_t n = x.size();
    2574                 :            : 
    2575                 :            :   // Define jacobian
    2576         [ -  - ]:          0 :   std::vector< double > jacob(n*n);
    2577                 :            : 
    2578                 :            :   // Compute f with initial guess
    2579 [ -  - ][ -  - ]:          0 :   std::vector< tk::real > f = DG::nonlinear_func(e, x);
    2580                 :            : 
    2581                 :            :   // Store the norm of f initially, for relative error measure
    2582                 :          0 :   tk::real err0 = 0.0;
    2583         [ -  - ]:          0 :   for (std::size_t i=0; i<n; ++i)
    2584                 :          0 :     err0 += f[i]*f[i];
    2585                 :          0 :   err0 = std::sqrt(err0);
    2586                 :          0 :   auto abs_err_old = err0;
    2587                 :            : 
    2588                 :            :   // Iterate for the solution if err0 > 0
    2589                 :          0 :   solver_failed = false;
    2590                 :          0 :   tk::real alpha_jacob = 1.0;
    2591         [ -  - ]:          0 :   if (err0 > abs_tol)
    2592         [ -  - ]:          0 :     for (std::size_t iter=0; iter<max_iter; ++iter)
    2593                 :            :     {
    2594                 :            : 
    2595                 :            :       // Evaluate jacobian
    2596                 :          0 :       tk::real dx = 0.0;
    2597         [ -  - ]:          0 :       for (std::size_t i=0; i<n; ++i)
    2598         [ -  - ]:          0 :         for (std::size_t j=0; j<n; ++j)
    2599                 :            :         {
    2600                 :            :           // Set dx in the order 1% of the unknown
    2601                 :          0 :           dx = alpha_jacob*std::max(std::abs(0.1*x[j]),1.0e-06);
    2602                 :            :           // Derivative of f[i] with respect to x[j]
    2603         [ -  - ]:          0 :           auto x_perturb = x;
    2604                 :          0 :           x_perturb[j] += dx;
    2605 [ -  - ][ -  - ]:          0 :           auto f_perturb = DG::nonlinear_func(e, x_perturb);
    2606                 :          0 :           jacob[i*n+j] = (f_perturb[i]-f[i])/dx;
    2607                 :            :         }
    2608                 :            : 
    2609                 :            :       // Compute new solution by solving linear system J*dx = -f
    2610                 :          0 :       lapack_int ln = static_cast< lapack_int >(n);
    2611         [ -  - ]:          0 :       std::vector< double > delta(n);
    2612         [ -  - ]:          0 :       for (std::size_t i=0; i<n; ++i)
    2613                 :          0 :         delta[i] = -f[i];
    2614                 :            :       lapack_int info;
    2615         [ -  - ]:          0 :       std::vector< lapack_int > ipiv(n);
    2616         [ -  - ]:          0 :       info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, ln, 1, jacob.data(), ln, ipiv.data(), delta.data(), 1);
    2617                 :            : 
    2618         [ -  - ]:          0 :       if (info != 0) {
    2619         [ -  - ]:          0 :         printf("Failed with info: %ld\n", info);
    2620                 :            :       }
    2621                 :            : 
    2622                 :            :       // Save f as fold
    2623         [ -  - ]:          0 :       std::vector< tk::real > fold(n);
    2624         [ -  - ]:          0 :       for (std::size_t i=0; i<n; ++i)
    2625                 :          0 :         fold[i] = f[i];
    2626                 :            : 
    2627                 :            :       // Update x using line search
    2628                 :          0 :       bool ls_failed = false;
    2629                 :          0 :       tk::real alpha_ls = 1.0E+00;
    2630                 :          0 :       std::size_t nline = 25;
    2631         [ -  - ]:          0 :       auto xtest = x;
    2632         [ -  - ]:          0 :       for (std::size_t iline = 0; iline<nline; ++iline)
    2633                 :            :       {
    2634                 :            :         // Evaluate xtest
    2635         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2636                 :          0 :           xtest[i] = x[i] + alpha_ls*delta[i];
    2637                 :            : 
    2638                 :            :         // Compute new f(x)
    2639 [ -  - ][ -  - ]:          0 :         f = DG::nonlinear_func(e, xtest);
    2640                 :            : 
    2641                 :          0 :         tk::real err = 0.0;
    2642         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2643                 :          0 :           err += f[i]*f[i];
    2644                 :          0 :         abs_err = std::sqrt(err);
    2645                 :            : 
    2646                 :            :         // If 1. The error went up
    2647                 :            :         // or 2. The function f flipped in sign
    2648                 :            :         // Reduce the factor alpha_ls
    2649                 :          0 :         bool flipped_sign = false;
    2650         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2651         [ -  - ]:          0 :           if (fold[i]*f[i] < 0.0) {
    2652                 :          0 :             flipped_sign = true;
    2653                 :          0 :             break;
    2654                 :            :           }
    2655                 :            : 
    2656 [ -  - ][ -  - ]:          0 :         if (abs_err < abs_err_old && !flipped_sign)
    2657                 :            :         {
    2658                 :            :           break;
    2659                 :            :         }
    2660                 :            :         else
    2661                 :            :         {
    2662                 :          0 :           alpha_ls *= 0.5;
    2663                 :            :         }
    2664         [ -  - ]:          0 :         if (iline == nline-1) {
    2665                 :            :           //printf("Line search failed to decrease f\n");
    2666                 :            :           // Try again by reducing the jacobian,
    2667                 :            :           // but only a few times, otherwise give up
    2668                 :          0 :           alpha_jacob *= 0.5;
    2669         [ -  - ]:          0 :           if (alpha_jacob < 1.0E-03)
    2670                 :          0 :             solver_failed = true;
    2671                 :            :           else
    2672                 :          0 :             ls_failed = true;
    2673                 :            :         }
    2674                 :            :       }
    2675                 :            : 
    2676         [ -  - ]:          0 :       if (solver_failed) {
    2677 [ -  - ][ -  - ]:          0 :         f = DG::nonlinear_func(e, x);
    2678         [ -  - ]:          0 :         printf("\nIMEX-RK: Non-linear solver did not converge in %lu iterations\n", iter+1);
    2679         [ -  - ]:          0 :         printf("Element #%lu\n", e);
    2680         [ -  - ]:          0 :         printf("Relative error: %e\n", rel_err);
    2681         [ -  - ]:          0 :         printf("Absolute error: %e\n\n", abs_err);
    2682                 :          0 :         break;
    2683                 :            :       }
    2684                 :            : 
    2685         [ -  - ]:          0 :       if (!ls_failed) {
    2686                 :            :         // Save x
    2687         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2688                 :          0 :           x[i] = xtest[i];
    2689                 :            : 
    2690                 :            :         // Compute a measure of error, use norm of f
    2691                 :          0 :         tk::real err = 0.0;
    2692         [ -  - ]:          0 :         for (std::size_t i=0; i<n; ++i)
    2693                 :          0 :           err += f[i]*f[i];
    2694                 :          0 :         abs_err = std::sqrt(err);
    2695                 :          0 :         rel_err = abs_err/err0;
    2696                 :            : 
    2697                 :            :         // check if error condition is met and loop back
    2698 [ -  - ][ -  - ]:          0 :         if (rel_err < rel_tol || abs_err < abs_tol) {
    2699                 :            :           break;
    2700                 :            :         }
    2701                 :            : 
    2702                 :            :         // If we did not converge, print a message and keep going
    2703         [ -  - ]:          0 :         if (iter == max_iter-1)
    2704                 :            :         {
    2705         [ -  - ]:          0 :           printf("\nIMEX-RK: Non-linear solver did not converge in %lu iterations\n", max_iter);
    2706         [ -  - ]:          0 :           printf("Element #%lu\n", e);
    2707         [ -  - ]:          0 :           printf("Relative error: %e\n", rel_err);
    2708         [ -  - ]:          0 :           printf("Absolute error: %e\n\n", abs_err);
    2709                 :            :         }
    2710                 :            :       }
    2711                 :            :     }
    2712                 :            : 
    2713                 :          0 :   return x;
    2714                 :            : 
    2715                 :            : }
    2716                 :            : 
    2717                 :            : #include "NoWarning/dg.def.h"

Generated by: LCOV version 1.14