p-adaptation finalized for multi-material

src/Inciter/DG.cpp

@@ -117,10 +117,12 @@ DG::DG( const CProxy_Discretization& disc,
   m_diag(),
   m_stage( 0 ),
   m_ndof(),
+  m_interface(),
   m_numEqDof(),
   m_uc(),
   m_pc(),
   m_ndofc(),
+  m_interfacec(),
   m_initial( 1 ),
   m_uElemfields( m_u.nunk(),
                  g_inputdeck.get< tag::ncomp >() ),
@@ -231,6 +233,7 @@ DG::resizeSolVectors()
   for (auto& n : m_ndofc) n.resize( myGhosts()->m_bid.size() );
   for (auto& u : m_uc) u.resize( myGhosts()->m_bid.size() );
   for (auto& p : m_pc) p.resize( myGhosts()->m_bid.size() );
+  for (auto& i : m_interfacec) i.resize( myGhosts()->m_bid.size() );
 
   // Initialize number of degrees of freedom in mesh elements
   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
@@ -244,6 +247,7 @@ DG::resizeSolVectors()
     const auto ndof = g_inputdeck.get< tag::ndof >();
     m_ndof.resize( myGhosts()->m_nunk, ndof );
   }
+  m_interface.resize( myGhosts()->m_nunk, 0 );
 
   // Ensure that we also have all the geometry and connectivity data
   // (including those of ghosts)
@@ -317,7 +321,7 @@ DG::box( tk::real v, const std::vector< tk::real >& )
     myGhosts()->m_coord, m_boxelems, d->ElemBlockId(), m_u, d->T(),
     myGhosts()->m_fd.Esuel().size()/4 );
   g_dgpde[d->MeshId()].updatePrimitives( m_u, m_lhs, myGhosts()->m_geoElem, m_p,
-    myGhosts()->m_fd.Esuel().size()/4 );
+    myGhosts()->m_fd.Esuel().size()/4, m_ndof );
 
   m_un = m_u;
 
@@ -404,18 +408,22 @@ DG::next()
       std::vector< std::size_t > tetid( ghostdata.size() );
       std::vector< std::vector< tk::real > > u( ghostdata.size() ),
                                              prim( ghostdata.size() );
-      std::vector< std::size_t > ndof;
+      std::vector< std::size_t > interface( ghostdata.size() );
+      std::vector< std::size_t > ndof( ghostdata.size() );
       std::size_t j = 0;
       for(const auto& i : ghostdata) {
         Assert( i < myGhosts()->m_fd.Esuel().size()/4,
           "Sending solution ghost data" );
         tetid[j] = i;
         u[j] = m_u[i];
         prim[j] = m_p[i];
-        if (pref && m_stage == 0) ndof.push_back( m_ndof[i] );
+        if (pref && m_stage == 0) {
+          ndof[j] = m_ndof[i];
+          interface[j] = m_interface[i];
+        }
         ++j;
       }
-      thisProxy[ cid ].comsol( thisIndex, m_stage, tetid, u, prim, ndof );
+      thisProxy[ cid ].comsol( thisIndex, m_stage, tetid, u, prim, interface, ndof );
     }
 
   ownsol_complete();
@@ -427,6 +435,7 @@ DG::comsol( int fromch,
             const std::vector< std::size_t >& tetid,
             const std::vector< std::vector< tk::real > >& u,
             const std::vector< std::vector< tk::real > >& prim,
+            const std::vector< std::size_t >& interface,
             const std::vector< std::size_t >& ndof )
 // *****************************************************************************
 //  Receive chare-boundary solution ghost data from neighboring chares
@@ -435,6 +444,7 @@ DG::comsol( int fromch,
 //! \param[in] tetid Ghost tet ids we receive solution data for
 //! \param[in] u Solution ghost data
 //! \param[in] prim Primitive variables in ghost cells
+//! \param[in] interface Interface marker in ghost cells
 //! \param[in] ndof Number of degrees of freedom for chare-boundary elements
 //! \details This function receives contributions to the unlimited solution
 //!   from fellow chares.
@@ -445,8 +455,10 @@ DG::comsol( int fromch,
 
   const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
 
-  if (pref && fromstage == 0)
+  if (pref && fromstage == 0) {
     Assert( ndof.size() == tetid.size(), "Size mismatch in DG::comsol()" );
+    Assert( interface.size() == tetid.size(), "Size mismatch in DG::comsol()" );
+  }
 
   // Find local-to-ghost tet id map for sender chare
   const auto& n = tk::cref_find( myGhosts()->m_ghost, fromch );
@@ -462,6 +474,8 @@ DG::comsol( int fromch,
     if (pref && fromstage == 0) {
       Assert( b < m_ndofc[0].size(), "Indexing out of bounds" );
       m_ndofc[0][b] = ndof[i];
+      Assert( b < m_interfacec[0].size(), "Indexing out of bounds" );
+      m_interfacec[0][b] = interface[i];
     }
   }
 
@@ -578,6 +592,7 @@ void DG::refine()
     }
     if (pref && m_stage == 0) {
       m_ndof[ b.first ] = m_ndofc[0][ b.second ];
+      m_interface[ b.first ] = m_interfacec[0][ b.second ];
     }
   }
 
@@ -734,12 +749,16 @@ DG::reco()
   }
 
   auto d = Disc();
+  if (pref && m_stage==0) {
+    g_dgpde[d->MeshId()].resetAdapSol( myGhosts()->m_fd, m_u, m_p, m_ndof );
+  }
+
   if (rdof > 1)
     // Reconstruct second-order solution and primitive quantities
     g_dgpde[d->MeshId()].reconstruct( d->T(), myGhosts()->m_geoFace,
       myGhosts()->m_geoElem,
       myGhosts()->m_fd, myGhosts()->m_esup, myGhosts()->m_inpoel,
-      myGhosts()->m_coord, m_u, m_p );
+      myGhosts()->m_coord, m_u, m_p, pref, m_ndof );
 
   // Send reconstructed solution to neighboring chares
   if (myGhosts()->m_sendGhost.empty())
@@ -1087,6 +1106,7 @@ DG::lim()
 {
   auto d = Disc();
   const auto rdof = g_inputdeck.get< tag::rdof >();
+  const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
   const auto ncomp = m_u.nprop() / rdof;
   const auto nprim = m_p.nprop() / rdof;
 
@@ -1127,10 +1147,11 @@ DG::lim()
   tk::destroy(m_pNodalExtrmc);
 
   if (rdof > 1) {
-    g_dgpde[d->MeshId()].limit( d->T(), myGhosts()->m_geoFace, myGhosts()->m_geoElem,
-              myGhosts()->m_fd, myGhosts()->m_esup, myGhosts()->m_inpoel,
-              myGhosts()->m_coord, m_ndof, d->Gid(), d->Bid(), m_uNodalExtrm,
-              m_pNodalExtrm, m_mtInv, m_u, m_p, m_shockmarker );
+    g_dgpde[d->MeshId()].limit( d->T(), pref, myGhosts()->m_geoFace,
+              myGhosts()->m_geoElem, myGhosts()->m_fd, myGhosts()->m_esup,
+              myGhosts()->m_inpoel, myGhosts()->m_coord, m_ndof, d->Gid(),
+              d->Bid(), m_uNodalExtrm, m_pNodalExtrm, m_mtInv, m_u, m_p,
+              m_shockmarker );
 
     if (g_inputdeck.get< tag::limsol_projection >())
       g_dgpde[d->MeshId()].CPL(m_p, myGhosts()->m_geoElem,
@@ -1171,8 +1192,9 @@ DG::refine_ndof()
 // *****************************************************************************
 {
   auto d = Disc();
+  const auto& coord = d->Coord();
   const auto& inpoel = d->Inpoel();
-  const auto npoin = d->Coord()[0].size();
+  const auto npoin = coord[0].size();
   const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
   std::vector<std::size_t> node_ndof(npoin, 1);
 
@@ -1219,7 +1241,8 @@ void DG::smooth_ndof()
 {
   auto d = Disc();
   const auto& inpoel = d->Inpoel();
-  const auto npoin = d->Coord()[0].size();
+  const auto& coord = d->Coord();
+  const auto npoin = coord[0].size();
   const auto nelem = myGhosts()->m_fd.Esuel().size()/4;
   std::vector<std::size_t> node_ndof(npoin, 1);
 
@@ -1374,42 +1397,11 @@ DG::solve( tk::real newdt )
   const auto rdof = g_inputdeck.get< tag::rdof >();
   const auto ndof = g_inputdeck.get< tag::ndof >();
   const auto neq = m_u.nprop()/rdof;
+  const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
 
   // Set new time step size
   if (m_stage == 0) d->setdt( newdt );
 
-  const auto pref = g_inputdeck.get< tag::pref, tag::pref >();
-  if (pref && m_stage == 0)
-  {
-    // When the element are coarsened, high order terms should be zero
-    for(std::size_t e = 0; e < myGhosts()->m_nunk; e++)
-    {
-      const auto ncomp= m_u.nprop()/rdof;
-      if(m_ndof[e] == 1)
-      {
-        for (std::size_t c=0; c<ncomp; ++c)
-        {
-          auto mark = c*rdof;
-          m_u(e, mark+1) = 0.0;
-          m_u(e, mark+2) = 0.0;
-          m_u(e, mark+3) = 0.0;
-        }
-      } else if(m_ndof[e] == 4)
-      {
-        for (std::size_t c=0; c<ncomp; ++c)
-        {
-          auto mark = c*ndof;
-          m_u(e, mark+4) = 0.0;
-          m_u(e, mark+5) = 0.0;
-          m_u(e, mark+6) = 0.0;
-          m_u(e, mark+7) = 0.0;
-          m_u(e, mark+8) = 0.0;
-          m_u(e, mark+9) = 0.0;
-        }
-      }
-    }
-  }
-
   // Update Un
   if (m_stage == 0) m_un = m_u;
 
@@ -1436,9 +1428,9 @@ DG::solve( tk::real newdt )
     }
   }
 
-  g_dgpde[d->MeshId()].rhs( physT, myGhosts()->m_geoFace, myGhosts()->m_geoElem,
-    myGhosts()->m_fd, myGhosts()->m_inpoel, m_boxelems, myGhosts()->m_coord,
-    m_u, m_p, m_ndof, d->Dt(), m_rhs );
+  g_dgpde[d->MeshId()].rhs( physT, pref, myGhosts()->m_geoFace,
+    myGhosts()->m_geoElem, myGhosts()->m_fd, myGhosts()->m_inpoel, m_boxelems,
+    myGhosts()->m_coord, m_u, m_p, m_ndof, d->Dt(), m_rhs );
 
   if (!imex_runge_kutta) {
     // Explicit time-stepping using RK3 to discretize time-derivative
@@ -1447,13 +1439,15 @@ DG::solve( tk::real newdt )
       {
         for (std::size_t k=0; k<m_numEqDof[c]; ++k)
         {
-          auto rmark = c*rdof+k;
-          auto mark = c*ndof+k;
-          m_u(e, rmark) =  rkcoef[0][m_stage] * m_un(e, rmark)
-            + rkcoef[1][m_stage] * ( m_u(e, rmark)
-              + d->Dt() * m_rhs(e, mark)/m_lhs(e, mark) );
-          if(fabs(m_u(e, rmark)) < 1e-16)
-            m_u(e, rmark) = 0;
+          if(k < m_ndof[e]) {
+            auto rmark = c*rdof+k;
+            auto mark = c*ndof+k;
+            m_u(e, rmark) =  rkcoef[0][m_stage] * m_un(e, rmark)
+              + rkcoef[1][m_stage] * ( m_u(e, rmark)
+                + d->Dt() * m_rhs(e, mark)/m_lhs(e, mark) );
+            if(fabs(m_u(e, rmark)) < 1e-16)
+              m_u(e, rmark) = 0;
+          }
         }
       }
   }
@@ -1478,9 +1472,9 @@ DG::solve( tk::real newdt )
 
   // Update primitives based on the evolved solution
   g_dgpde[d->MeshId()].updateInterfaceCells( m_u,
-    myGhosts()->m_fd.Esuel().size()/4, m_ndof );
+    myGhosts()->m_fd.Esuel().size()/4, m_ndof, m_interface );
   g_dgpde[d->MeshId()].updatePrimitives( m_u, m_lhs, myGhosts()->m_geoElem, m_p,
-    myGhosts()->m_fd.Esuel().size()/4 );
+    myGhosts()->m_fd.Esuel().size()/4, m_ndof );
   if (!g_inputdeck.get< tag::accuracy_test >()) {
     g_dgpde[d->MeshId()].cleanTraceMaterial( physT, myGhosts()->m_geoElem, m_u,
       m_p, myGhosts()->m_fd.Esuel().size()/4 );
@@ -1743,6 +1737,10 @@ DG::writeFields(
   if (g_inputdeck.get< tag::pref, tag::pref >()) {
     std::vector< tk::real > ndof( begin(m_ndof), end(m_ndof) );
     ndof.resize( nelem );
+    for(std::size_t k = 0; k < nelem; k++) {
+      // Mark the cell with THINC reconstruction as 0 for output
+      if(m_interface[k] == 1) ndof[k] = 0;
+    }
     for (const auto& [child,parent] : addedTets)
       ndof[child] = static_cast< tk::real >( m_ndof[parent] );
     elemfields.push_back( ndof );
@@ -1778,8 +1776,9 @@ DG::writeFields(
   analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::ELEM, elemfieldnames );
   analyticFieldNames( g_dgpde[d->MeshId()], tk::Centering::NODE, nodefieldnames );
 
-  if (g_inputdeck.get< tag::pref, tag::pref >())
+  if (g_inputdeck.get< tag::pref, tag::pref >()) {
     elemfieldnames.push_back( "NDOF" );
+  }
 
   elemfieldnames.push_back( "shock_marker" );
 

src/Inciter/DG.hpp

@@ -135,6 +135,7 @@ class DG : public CBase_DG {
                  const std::vector< std::size_t >& tetid,
                  const std::vector< std::vector< tk::real > >& u,
                  const std::vector< std::vector< tk::real > >& prim,
+                 const std::vector< std::size_t >& interface,
                  const std::vector< std::size_t >& ndof );
 
     //! Receive contributions to nodal gradients on chare-boundaries
@@ -259,10 +260,12 @@ class DG : public CBase_DG {
       p | m_diag;
       p | m_stage;
       p | m_ndof;
+      p | m_interface;
       p | m_numEqDof;
       p | m_uc;
       p | m_pc;
       p | m_ndofc;
+      p | m_interfacec;
       p | m_initial;
       p | m_uElemfields;
       p | m_pElemfields;
@@ -350,6 +353,8 @@ class DG : public CBase_DG {
     std::size_t m_stage;
     //! Vector of local number of degrees of freedom for each element
     std::vector< std::size_t > m_ndof;
+    //! Interface marker for field output
+    std::vector< std::size_t > m_interface;
     //! Vector of number of degrees of freedom for each PDE equation/component
     std::vector< std::size_t > m_numEqDof;
     //! Solution receive buffers for ghosts only
@@ -359,6 +364,8 @@ class DG : public CBase_DG {
     //! \brief Number of degrees of freedom (for p-adaptive) receive buffers
     //!   for ghosts only
     std::array< std::vector< std::size_t >, 3 > m_ndofc;
+    //! Interface marker receive buffers for ghost only
+    std::array< std::vector< std::size_t >, 1 > m_interfacec;
     //! 1 if starting time stepping, 0 if during time stepping
     std::size_t m_initial;
     //! Solution elem output fields

src/Inciter/dg.ci

@@ -57,6 +57,7 @@ module dg {
                          const std::vector< std::size_t >& tetid,
                          const std::vector< std::vector< tk::real > >& u,
                          const std::vector< std::vector< tk::real > >& prim,
+                         const std::vector< std::size_t >& interface,
                          const std::vector< std::size_t >& ndof );
       entry void refine( const std::vector< tk::real >& l2ref );
       entry [reductiontarget] void solve( tk::real newdt );