Quinoa all test code coverage report
Current view: top level - Inciter - DG.cpp (source / functions) Hit Total Coverage
Commit: -128-NOTFOUND Lines: 824 1219 67.6 %
Date: 2025-10-02 14:31:17 Functions: 39 45 86.7 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 815 2064 39.5 %

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

Generated by: LCOV version 1.14