Header path

fem/nonlinearform_ext.cpp

@@ -15,6 +15,9 @@
 #include "nonlinearform.hpp"
 #include "ceed/interface/util.hpp"
 
+#define NVTX_COLOR ::gpu::nvtx::color_names::kLime
+#include "../general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -40,13 +43,17 @@ PANonlinearFormExtension::PANonlinearFormExtension(const NonlinearForm *nlf):
 
 real_t PANonlinearFormExtension::GetGridFunctionEnergy(const Vector &x) const
 {
+   // dbg("x: {}",x*x);
    real_t energy = 0.0;
 
    elemR->Mult(x, xe);
+   // dbg("xe: {}",xe*xe);
+
    for (int i = 0; i < dnfi.Size(); i++)
    {
       energy += dnfi[i]->GetLocalStateEnergyPA(xe);
    }
+   dbg("energy: {}",energy);
    return energy;
 }
 

fem/tmop.cpp

@@ -15,6 +15,9 @@
 #include "tmop_tools.hpp"
 #include "../general/forall.hpp"
 
+#define NVTX_COLOR ::gpu::nvtx::color_names::kMoccasin
+#include "../general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -1565,6 +1568,8 @@ void TargetConstructor::ComputeAllElementTargets_Fallback(
    const Vector &xe,
    DenseTensor &Jtr) const
 {
+   dbg("xe:{} {}",xe*xe, xe.Norml2());
+
    // Fallback to the 1-element method, ComputeElementTargets()
 
    // When UsesPhysicalCoordinates() == true, we assume 'xe' uses
@@ -1584,8 +1589,17 @@ void TargetConstructor::ComputeAllElementTargets_Fallback(
    MFEM_VERIFY(!UsesPhysicalCoordinates() ||
                xe.Size() == NE*nvdofs, "invalid input Vector 'xe'!");
    const int NQ = ir.GetNPoints();
-   const Array<int> *dof_map = nullptr;
-   if (UsesPhysicalCoordinates())
+   // const Array<int> *dof_map = nullptr;
+
+   dbg("UsesPhysicalCoordinates:{}",UsesPhysicalCoordinates());
+   /*{
+      assert(sdim==2); assert(dim==2);
+      assert(Jtr.TotalSize() == 2*2*NQ*NE);
+      Vector jtr(Jtr.HostReadWrite(), 2* 2*NQ*NE);
+      dbg("Jtr:{}",jtr*jtr);
+   }*/
+
+   /*if (UsesPhysicalCoordinates())
    {
       const TensorBasisElement *tfe =
          dynamic_cast<const TensorBasisElement *>(&fe);
@@ -1594,23 +1608,25 @@ void TargetConstructor::ComputeAllElementTargets_Fallback(
          dof_map = &tfe->GetDofMap();
          if (dof_map->Size() == 0) { dof_map = nullptr; }
       }
-   }
+   }*/
 
    Vector elfun_lex, elfun_nat;
    DenseTensor J;
    xe.HostRead();
-   Jtr.HostWrite();
-   if (UsesPhysicalCoordinates() && dof_map != nullptr)
+   Jtr.HostReadWrite();
+   /*if (UsesPhysicalCoordinates() && dof_map != nullptr)
    {
       elfun_nat.SetSize(nvdofs);
-   }
+   }*/
    for (int e = 0; e < NE; e++)
    {
-      if (UsesPhysicalCoordinates())
+      /*if (UsesPhysicalCoordinates())
       {
+         assert(false);
          if (!dof_map)
          {
             elfun_nat.SetDataAndSize(xe.GetData()+e*nvdofs, nvdofs);
+            dbg("elfun_nat!:{}",elfun_nat*elfun_nat);
          }
          else
          {
@@ -1624,11 +1640,22 @@ void TargetConstructor::ComputeAllElementTargets_Fallback(
                      elfun_lex[i_lex+d*ndofs];
                }
             }
+            dbg("elfun_nat:{}",elfun_nat*elfun_nat);
          }
-      }
+      }*/
+
+#warning 🔥
       J.UseExternalData(Jtr(e*NQ).Data(), sdim, dim, NQ);
       ComputeElementTargets(e, fe, ir, elfun_nat, J);
+      // Jtr.Read();
+
+      // Vector j(J.HostReadWrite(), sdim* dim*NQ);
+      // dbg("\x1B[32m j:{}",j*j);
    }
+
+   // const auto J = Reshape(J.Read(), 2, 2, Q1D, Q1D, NE);
+   Vector jtr(Jtr.HostReadWrite(), 2* 2* NQ* NE);
+   dbg("\x1B[32m Jtr:{}",jtr*jtr);
 }
 
 bool TargetConstructor::ContainsVolumeInfo() const
@@ -1650,6 +1677,7 @@ void TargetConstructor::ComputeElementTargets(int e_id, const FiniteElement &fe,
                                               const Vector &elfun,
                                               DenseTensor &Jtr) const
 {
+   // dbg("e_id:{} target_type:{}",e_id, static_cast<uint32_t>(target_type));
    MFEM_CONTRACT_VAR(elfun);
    MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != NULL, "");
 
@@ -1744,6 +1772,7 @@ void AnalyticAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
                                             const Vector &elfun,
                                             DenseTensor &Jtr) const
 {
+   dbg("e_id:{} target_type:{}",e_id, static_cast<uint32_t>(target_type));
    DenseMatrix point_mat;
    point_mat.UseExternalData(elfun.GetData(), fe.GetDof(), fe.GetDim());
 
@@ -1822,6 +1851,7 @@ static inline void device_copy(real_t *d_dest, const real_t *d_src, int size)
 #ifdef MFEM_USE_MPI
 void DiscreteAdaptTC::FinalizeParDiscreteTargetSpec(const ParGridFunction &t)
 {
+   dbg();
    MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
    MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
 
@@ -1847,6 +1877,7 @@ void DiscreteAdaptTC::FinalizeParDiscreteTargetSpec(const ParGridFunction &t)
 
 void DiscreteAdaptTC::ParUpdateAfterMeshTopologyChange()
 {
+   dbg();
    ptspec_fesv->Update();
    if (tspec_fesv)
    {
@@ -1865,6 +1896,7 @@ void DiscreteAdaptTC::ParUpdateAfterMeshTopologyChange()
 
 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const ParGridFunction &tspec_)
 {
+   dbg();
    const int vdim = tspec_.FESpace()->GetVDim(),
              ndof = tspec_.FESpace()->GetNDofs();
    MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTspecAtIndex.");
@@ -1924,6 +1956,7 @@ void DiscreteAdaptTC::SetParDiscreteTargetSpec(const ParGridFunction &tspec_)
 
 void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
 {
+   dbg();
    const int vdim = tspec_.FESpace()->GetVDim(),
              ndof = tspec_.FESpace()->GetNDofs();
    ncomp += vdim;
@@ -1947,6 +1980,7 @@ void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
 
 void DiscreteAdaptTC::SetTspecAtIndex(int idx, const GridFunction &tspec_)
 {
+   dbg();
    const int vdim = tspec_.FESpace()->GetVDim(),
              ndof = tspec_.FESpace()->GetNDofs();
    MFEM_VERIFY(ndof == tspec.Size()/ncomp, "Inconsistency in SetTspecAtIndex.");
@@ -2000,6 +2034,7 @@ void DiscreteAdaptTC::SetSerialDiscreteTargetOrientation(const GridFunction &o)
 
 void DiscreteAdaptTC::FinalizeSerialDiscreteTargetSpec(const GridFunction &t)
 {
+   dbg();
    MFEM_VERIFY(adapt_eval, "SetAdaptivityEvaluator() has not been called!")
    MFEM_VERIFY(ncomp > 0, "No target specifications have been set!");
 
@@ -2026,7 +2061,7 @@ void DiscreteAdaptTC::GetDiscreteTargetSpec(GridFunction &tspec_, int idx)
              vdim = tspec_.FESpace()->GetVDim();
    MFEM_VERIFY(ndof == tspec.Size()/ncomp,
                "Inconsistency in GetSerialDiscreteTargetSpec.");
-
+   assert(false);
    for (int i = 0; i < ndof*vdim; i++)
    {
       tspec_(i) = tspec(i + idx*ndof);
@@ -2035,6 +2070,7 @@ void DiscreteAdaptTC::GetDiscreteTargetSpec(GridFunction &tspec_, int idx)
 
 void DiscreteAdaptTC::UpdateAfterMeshTopologyChange()
 {
+   dbg();
    tspec_fesv->Update();
    tspec_gf->Update();
    tspec.SetDataAndSize(tspec_gf->GetData(), tspec_gf->Size());
@@ -2054,19 +2090,25 @@ void DiscreteAdaptTC::UpdateTargetSpecification(const Vector &new_x,
                                                 bool reuse_flag,
                                                 int new_x_ordering)
 {
+   dbg();
+   dbg("\x1B[33m tspec:{} {}",tspec*tspec, tspec.Norml2());
    if (reuse_flag && good_tspec) { return; }
 
+   dbg("new_x:{}",new_x*new_x);
+
    MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
    adapt_eval->ComputeAtNewPosition(new_x, tspec, new_x_ordering);
    tspec_sav = tspec;
 
    good_tspec = reuse_flag;
+   dbg("\x1B[33m tspec:{} {}",tspec*tspec, tspec.Norml2());
 }
 
 void DiscreteAdaptTC::UpdateTargetSpecification(Vector &new_x,
                                                 Vector &IntData,
                                                 int new_x_ordering)
 {
+   dbg();
    adapt_eval->ComputeAtNewPosition(new_x, IntData, new_x_ordering);
 }
 
@@ -2075,6 +2117,7 @@ void DiscreteAdaptTC::UpdateTargetSpecificationAtNode(const FiniteElement &el,
                                                       int dofidx, int dir,
                                                       const Vector &IntData)
 {
+   dbg();
    MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
 
    Array<int> dofs;
@@ -2090,11 +2133,13 @@ void DiscreteAdaptTC::UpdateTargetSpecificationAtNode(const FiniteElement &el,
 void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
                                                        int dofidx)
 {
+   dbg();
    MFEM_VERIFY(tspec.Size() > 0, "Target specification is not set!");
 
    Array<int> dofs;
    tspec_fesv->GetElementDofs(T.ElementNo, dofs);
    const int cnt = tspec.Size()/ncomp;
+   assert(false);
    for (int i = 0; i < ncomp; i++)
    {
       tspec(dofs[dofidx] + i*cnt) = tspec_sav(dofs[dofidx] + i*cnt);
@@ -2104,6 +2149,7 @@ void DiscreteAdaptTC::RestoreTargetSpecificationAtNode(ElementTransformation &T,
 void DiscreteAdaptTC::SetTspecFromIntRule(int e_id,
                                           const IntegrationRule &intrule)
 {
+   dbg();
    switch (target_type)
    {
       case IDEAL_SHAPE_GIVEN_SIZE:
@@ -2140,11 +2186,13 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
                                             const Vector &elfun,
                                             DenseTensor &Jtr) const
 {
+   // dbg("e_id:{} target_type:{}",e_id, static_cast<uint32_t>(target_type));
    MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
    const int dim = fe.GetDim(),
              nqp = ir.GetNPoints();
    Jtrcomp.SetSize(dim, dim, 4*nqp);
 
+   // dbg("tspec:{} {}",tspec*tspec, tspec.Norml2());
    FiniteElementSpace *src_fes = tspec_fesv;
 
    switch (target_type)
@@ -2163,8 +2211,10 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
          Array<int> dofs;
          DenseMatrix D_rho(dim), Q_phi(dim), R_theta(dim);
          tspec_fesv->GetElementVDofs(e_id, dofs);
-         tspec.UseDevice(true);
+         // tspec.UseDevice(true);
+         tspec.HostRead();
          tspec.GetSubVector(dofs, tspec_vals);
+         // dbg("tspec_vals:{} {}",tspec_vals*tspec_vals, tspec_vals.Norml2());
          if (tspec_refine.NumCols() > 0) // Refinement
          {
             MFEM_VERIFY(amr_el >= 0, " Target being constructed for an AMR element.");
@@ -2701,6 +2751,7 @@ void DiscreteAdaptTC::
 UpdateGradientTargetSpecification(const Vector &x, real_t dx,
                                   bool reuse_flag, int x_ordering)
 {
+   dbg();
    if (reuse_flag && good_tspec_grad) { return; }
 
    const int dim = tspec_fesv->GetFE(0)->GetDim(),
@@ -2738,6 +2789,7 @@ void DiscreteAdaptTC::
 UpdateHessianTargetSpecification(const Vector &x, real_t dx,
                                  bool reuse_flag, int x_ordering)
 {
+   dbg();
    if (reuse_flag && good_tspec_hess) { return; }
 
    const int dim    = tspec_fesv->GetFE(0)->GetDim(),

fem/tmop/tmop_pa.cpp

@@ -17,6 +17,10 @@
 #include "../../general/forall.hpp"
 #include "../../linalg/kernels.hpp"
 
+#define NVTX_COLOR ::gpu::nvtx::color_names::kYellow
+#include "../../general/nvtx.hpp"
+
+
 namespace mfem
 {
 
@@ -126,12 +130,14 @@ void TargetConstructor::ComputeAllElementTargets(const FiniteElementSpace &fes,
                                                  const Vector &xe,
                                                  DenseTensor &Jtr) const
 {
+   dbg("xe:{}", xe.Norml2());
    MFEM_VERIFY(Jtr.SizeI() == Jtr.SizeJ() && Jtr.SizeI() > 1, "");
    const int dim = Jtr.SizeI();
    bool done = false;
    if (dim == 2) { done = ComputeAllElementTargets<2>(fes, ir, xe, Jtr); }
    if (dim == 3) { done = ComputeAllElementTargets<3>(fes, ir, xe, Jtr); }
 
+   dbg("done:{}", done);
    if (!done) { ComputeAllElementTargets_Fallback(fes, ir, xe, Jtr); }
 }
 
@@ -140,6 +146,7 @@ void AnalyticAdaptTC::ComputeAllElementTargets(const FiniteElementSpace &fes,
                                                const Vector &xe,
                                                DenseTensor &Jtr) const
 {
+   dbg();
    ComputeAllElementTargets_Fallback(fes, ir, xe, Jtr);
 }
 
@@ -166,6 +173,7 @@ void AnalyticAdaptTC::ComputeAllElementTargets(const FiniteElementSpace &fes,
 //                                    Jtr(i) *= R_theta        (orientation)
 void TMOP_Integrator::ComputeAllElementTargets(const Vector &xe) const
 {
+   dbg("xe:{}",xe.Norml2());
    PA.Jtr_needs_update = false;
    PA.Jtr_debug_grad = false;
    const FiniteElementSpace *fes = PA.fes;
@@ -174,16 +182,23 @@ void TMOP_Integrator::ComputeAllElementTargets(const Vector &xe) const
 
    // Compute PA.Jtr for all elements
    targetC->ComputeAllElementTargets(*fes, ir, xe, PA.Jtr);
+
+   const int NE = PA.ne;
+   const int Q1D = PA.maps->nqpt;
+   DenseTensor &J = PA.Jtr;
+   Vector j(J.HostReadWrite(), 2* 2* Q1D* Q1D* NE);
+   dbg("j:{}",j*j);
 }
 
 void TMOP_Integrator::UpdateCoefficientsPA(const Vector &x_loc)
 {
+   x_loc.HostRead();
    // Both are constant or not specified.
    if (PA.MC.Size() == 1 && PA.C0.Size() == 1) { return; }
 
    // Coefficients are always evaluated on the CPU for now.
-   PA.MC.HostWrite();
-   PA.C0.HostWrite();
+   PA.MC.HostReadWrite();
+   PA.C0.HostReadWrite();
 
    const IntegrationRule &ir = *PA.ir;
    auto T = new IsoparametricTransformation;
@@ -361,6 +376,7 @@ void TMOP_Integrator::AddMultGradPA(const Vector &re, Vector &ce) const
 
 real_t TMOP_Integrator::GetLocalStateEnergyPA(const Vector &xe) const
 {
+   // dbg("xe:{} lim_coeff:{}", xe*xe, fmt::ptr(lim_coeff));
    // This method must be called after AssemblePA().
 
    real_t energy = 0.0;
@@ -373,7 +389,9 @@ real_t TMOP_Integrator::GetLocalStateEnergyPA(const Vector &xe) const
    if (PA.dim == 2)
    {
       energy = GetLocalStateEnergyPA_2D(xe);
+      // dbg("energy: {}", energy);
       if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_2D(xe); }
+      // dbg("energy: {}", energy);
    }
 
    if (PA.dim == 3)
@@ -382,6 +400,7 @@ real_t TMOP_Integrator::GetLocalStateEnergyPA(const Vector &xe) const
       if (lim_coeff) { energy += GetLocalStateEnergyPA_C0_3D(xe); }
    }
 
+   // dbg("energy: {}", energy);
    return energy;
 }
 

fem/tmop/tmop_pa.hpp

@@ -167,7 +167,7 @@ struct K##kernel##_T {\
      I==0 ? 0x22 : I==1 ? 0x23 : I==2 ? 0x24 : I==3 ? 0x25 : I==4 ? 0x26 :\
      I==5 ? 0x33 : I==6 ? 0x34 : I==7 ? 0x35 : I==8 ? 0x36  :\
      I==9 ? 0x44 : I==10 ? 0x45 : I==11 ? 0x46 :\
-     I==12 ? 0x55 : I==13 ? 0x56 : 0; }\
+     I==12 ? 0x55 : I==13 ? 0x56 :0; }\
    template<Key_t K> static constexpr Kernel_t GetKer() noexcept\
    { return &kernel<(K>>4)&0xF, K&0xF>; }\
 };\

fem/tmop/tmop_pa_tc2.cpp

@@ -16,6 +16,9 @@
 #include "../../general/forall.hpp"
 #include "../../linalg/kernels.hpp"
 
+#define NVTX_COLOR ::gpu::nvtx::color_names::kTomato
+#include "../../general/nvtx.hpp"
+
 using namespace mfem;
 
 namespace mfem
@@ -115,6 +118,7 @@ TargetConstructor::ComputeAllElementTargets<2>(const FiniteElementSpace &fes,
                                                const Vector &,
                                                DenseTensor &Jtr) const
 {
+   dbg();
    MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != nullptr, "");
    const Mesh *mesh = fes.GetMesh();
    const int NE = mesh->GetNE();

fem/tmop/tmop_pa_w2.cpp

@@ -16,6 +16,9 @@
 #include "../../linalg/kernels.hpp"
 #include "../../linalg/dinvariants.hpp"
 
+#define NVTX_COLOR ::gpu::nvtx::color_names::kLightSkyBlue
+#include "../../general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -160,7 +163,6 @@ MFEM_REGISTER_TMOP_KERNELS(real_t, EnergyPA_2D,
                mid == 77 ? EvalW_077(Jpt) :
                mid == 80 ? EvalW_080(Jpt, metric_data) :
                mid == 94 ? EvalW_094(Jpt, metric_data) : 0.0;
-
             E(qx,qy,e) = weight * EvalW;
          }
       }

fem/tmop_tools.cpp

@@ -14,6 +14,9 @@
 #include "pnonlinearform.hpp"
 #include "../general/osockstream.hpp"
 
+#define NVTX_COLOR ::gpu::nvtx::color_names::kCyan
+#include "../general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -30,35 +33,49 @@ void AdvectorCG::ComputeAtNewPosition(const Vector &new_nodes,
                                       Vector &new_field,
                                       int new_nodes_ordering)
 {
+   dbg("new_nodes: {}",new_nodes*new_nodes);
+   dbg("field0: {}",field0*field0);
    FiniteElementSpace *space = fes;
 #ifdef MFEM_USE_MPI
    if (pfes) { space = pfes; }
 #endif
    int fes_ordering = space->GetOrdering(),
        ncomp = space->GetVDim();
+   dbg("fes_ordering:{}",fes_ordering);
 
    // TODO: Implement for AMR meshes.
    const int pnt_cnt = field0.Size() / ncomp;
 
    new_field = field0;
+   dbg("new_field: {}",new_field*new_field);
    Vector new_field_temp;
+   new_field_temp =  new_field;
+
    for (int i = 0; i < ncomp; i++)
    {
       if (fes_ordering == Ordering::byNODES)
       {
-         new_field_temp.MakeRef(new_field, i*pnt_cnt, pnt_cnt);
+         new_field_temp.MakeRef(new_field, i*pnt_cnt, pnt_cnt); // 🔥🔥🔥🔥
       }
       else
       {
+         assert(false);
          new_field_temp.SetSize(pnt_cnt);
          for (int j = 0; j < pnt_cnt; j++)
          {
             new_field_temp(j) = new_field(i + j*ncomp);
          }
       }
+      dbg("new_nodes: {}",new_nodes*new_nodes);
+      dbg("new_field_temp: {}",new_field_temp*new_field_temp);
+      dbg("new_field: {}",new_field*new_field);
       ComputeAtNewPositionScalar(new_nodes, new_field_temp);
+      // dbg("new_field_temp: {}",new_field_temp*new_field_temp);
+      dbg("new_field: {}",new_field*new_field);
+
       if (fes_ordering == Ordering::byVDIM)
       {
+         assert(false);
          for (int j = 0; j < pnt_cnt; j++)
          {
             new_field(i + j*ncomp) = new_field_temp(j);
@@ -67,12 +84,15 @@ void AdvectorCG::ComputeAtNewPosition(const Vector &new_nodes,
    }
 
    field0 = new_field;
+   dbg("field0: {}",field0*field0);
    nodes0 = new_nodes;
+   dbg("nodes0: {}",nodes0*nodes0);
 }
 
 void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
                                             Vector &new_field)
 {
+   dbg();
    Mesh *m = mesh;
 #ifdef MFEM_USE_MPI
    if (pmesh) { m = pmesh; }
@@ -84,6 +104,8 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
    GridFunction *mesh_nodes = m->GetNodes();
    *mesh_nodes = nodes0;
    real_t minv = new_field.Min(), maxv = new_field.Max();
+   dbg("\x1B[32m minv:{}",minv);
+   dbg("\x1B[32m maxv:{}",minv);
 
    // Velocity of the positions.
    GridFunction u(mesh_nodes->FESpace());
@@ -117,9 +139,11 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
    {
       h_min = std::min(h_min, m->GetElementSize(i));
    }
+   dbg("\x1B[32m h_min:{}",h_min);
    real_t v_max = 0.0;
    const int s  = u.Size()/m->Dimension();
 
+   dbg("\x1B[32m u:{}",u*u);
    u.HostReadWrite();
    for (int i = 0; i < s; i++)
    {
@@ -130,6 +154,7 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
       }
       v_max = std::max(v_max, vel);
    }
+   dbg("\x1B[32m v_max:{}",v_max);
 
 #ifdef MFEM_USE_MPI
    if (pfes)
@@ -141,6 +166,8 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
                     pfes->GetComm());
    }
 #endif
+   dbg("\x1B[32m h_min:{}",h_min);
+   dbg("\x1B[32m v_max:{}",v_max);
 
    if (v_max == 0.0) // No need to change the field.
    {
@@ -154,7 +181,9 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
 
    v_max = std::sqrt(v_max);
    real_t dt = dt_scale * h_min / v_max;
+   dbg("\x1B[32m dt:{}",dt);
 
+   dbg("\x1B[32m new_field:{}",new_field*new_field);
    real_t t = 0.0;
    bool last_step = false;
    while (!last_step)
@@ -166,6 +195,7 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
       }
       ode_solver.Step(new_field, t, dt);
    }
+   dbg("\x1B[32m new_field:{}",new_field*new_field);
 
    real_t glob_minv = minv,
           glob_maxv = maxv;
@@ -178,6 +208,8 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
                     pfes->GetComm());
    }
 #endif
+   dbg("\x1B[32m glob_minv:{}",glob_minv);
+   dbg("\x1B[32m glob_maxv:{}",glob_maxv);
 
    // Trim the overshoots and undershoots.
    new_field.HostReadWrite();
@@ -186,6 +218,8 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
       if (new_field(i) < glob_minv) { new_field(i) = glob_minv; }
       if (new_field(i) > glob_maxv) { new_field(i) = glob_maxv; }
    }
+   new_field.Read(); // 🔥🔥🔥🔥
+   // dbg("\x1B[32m new_field:{}",new_field*new_field);
 
    delete oper;
    delete fess;
@@ -346,6 +380,7 @@ void ParAdvectorCGOper::Mult(const Vector &ind, Vector &di_dt) const
 void InterpolatorFP::SetInitialField(const Vector &init_nodes,
                                      const Vector &init_field)
 {
+   dbg();
    nodes0 = init_nodes;
    Mesh *m = mesh;
    FiniteElementSpace *f = fes;
@@ -353,6 +388,7 @@ void InterpolatorFP::SetInitialField(const Vector &init_nodes,
    if (pmesh) { m = pmesh; }
    if (pfes)  { f = pfes; }
 #endif
+   nodes0.HostReadWrite();
    m->SetNodes(nodes0);
 
    const real_t rel_bbox_el = 0.1;
@@ -390,11 +426,20 @@ void InterpolatorFP::ComputeAtNewPosition(const Vector &new_nodes,
                                           Vector &new_field,
                                           int new_nodes_ordering)
 {
+   dbg("\x1B[32m new_nodes:{}",new_nodes*new_nodes);
+   new_field.HostReadWrite();
+
    // Get physical node locations corresponding to field0_gf
    if (fes_field_nodes)
    {
       Vector mapped_nodes;
       GetFieldNodesPosition(new_nodes, mapped_nodes);
+#warning 🔥 gslib
+      mapped_nodes.HostReadWrite();
+      field0_gf.HostReadWrite();
+      new_field.HostReadWrite();
+      // finder->SetGPUtoCPUFallback(true);
+
       finder->Interpolate(mapped_nodes, field0_gf, new_field,
                           fes_field_nodes->GetOrdering());
    }
@@ -402,11 +447,13 @@ void InterpolatorFP::ComputeAtNewPosition(const Vector &new_nodes,
    {
       finder->Interpolate(new_nodes, field0_gf, new_field, new_nodes_ordering);
    }
+   dbg("\x1B[32m  new_field:{}",new_field*new_field);
 }
 
 void InterpolatorFP::GetFieldNodesPosition(const Vector &mesh_nodes,
                                            Vector &nodes_pos) const
 {
+   dbg();
    MFEM_VERIFY(fes_field_nodes, "InterpolatorFP: fes_field_nodes is not set.");
 
    Mesh *m = fes_field_nodes->GetMesh();
@@ -416,8 +463,15 @@ void InterpolatorFP::GetFieldNodesPosition(const Vector &mesh_nodes,
    if (nelem == 0) { return; }
    Array<int> dofs;
    Vector e_xyz;
+
+#warning 🔥
+   dbg("\x1B[37m mesh_nodes:{}",mesh_nodes*mesh_nodes);
+   mesh_nodes.HostRead();
+
    nodes_pos.SetSize(n_f_nodes*dim);
    nodes_pos.UseDevice(mesh_nodes.UseDevice());
+   nodes_pos.HostReadWrite();
+
    const FiniteElementSpace *mesh_fes = m->GetNodalFESpace();
 
    for (int e = 0; e < nelem; e++)
@@ -444,6 +498,8 @@ void InterpolatorFP::GetFieldNodesPosition(const Vector &mesh_nodes,
       fes_field_nodes->GetElementVDofs(e, dofs);
       nodes_pos.SetSubVector(dofs, gf_xyz);
    }
+   nodes_pos.Read();
+   dbg("\x1B[37m nodes_pos:{}",nodes_pos*nodes_pos);
 }
 
 #endif
@@ -451,6 +507,8 @@ void InterpolatorFP::GetFieldNodesPosition(const Vector &mesh_nodes,
 real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
                                               const Vector &b) const
 {
+   NVTX_MARK_FUNCTION;
+   dbg("x:{} b:{}", x*x, b*b);
    const FiniteElementSpace *fes = NULL;
    real_t energy_in = 0.0;
 #ifdef MFEM_USE_MPI
@@ -470,6 +528,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       fes = nlf->FESpace();
       energy_in = nlf->GetEnergy(x);
    }
+   dbg("\x1B[33m energy_in: {}", energy_in); // ✅
 
    // Get the local prolongation of the solution vector.
    Vector x_out_loc(fes->GetVSize(),
@@ -486,12 +545,14 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       fes->GetProlongationMatrix()->Mult(x, x_out_loc);
    }
 #endif
+   dbg("\x1B[31m x_out_loc: {}", x_out_loc*x_out_loc); //
 
    real_t scale = 1.0;
    bool fitting = IsSurfaceFittingEnabled();
    real_t init_fit_avg_err, init_fit_max_err = 0.0;
    if (fitting && surf_fit_converge_error)
    {
+      assert(false);
       GetSurfaceFittingError(x_out_loc, init_fit_avg_err, init_fit_max_err);
       // Check for convergence
       if (init_fit_max_err < surf_fit_max_err_limit)
@@ -521,6 +582,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
    // Check if the starting mesh (given by x) is inverted. Note that x hasn't
    // been modified by the Newton update yet.
    const real_t min_detT_in = ComputeMinDet(x_out_loc, *fes);
+   dbg("\x1B[33m min_detT_in: {}", min_detT_in); // ✅
    const bool untangling = (min_detT_in <= 0.0) ? true : false;
    const real_t untangle_factor = 1.5;
    if (untangling)
@@ -551,6 +613,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
    // Perform the line search.
    for (int i = 0; i < 12; i++)
    {
+      dbg("i:{}", i);
       avg_fit_err = 0.0;
       max_fit_err = 0.0;
 
@@ -565,9 +628,12 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
 #ifdef MFEM_USE_MPI
       else { fes->GetProlongationMatrix()->Mult(x_out, x_out_loc); }
 #endif
+      dbg("\x1B[34m x_out_loc: {}", x_out_loc*x_out_loc); //
 
       // Check the changes in detJ.
       min_detT_out = ComputeMinDet(x_out_loc, *fes);
+      dbg("\x1B[33m min_detT_out: {}", min_detT_out); //
+
       if (untangling == false && min_detT_out <= min_detJ_limit)
       {
          // No untangling, and detJ got negative (or small) -- no good.
@@ -594,11 +660,13 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
 
       // Check the changes in total energy.
       ProcessNewState(x_out);
+      dbg("\x1B[33m x_out: {}", x_out*x_out); //
 
       // Ensure sufficient decrease in fitting error if we are trying to
       // converge based on error.
       if (fitting && surf_fit_converge_error)
       {
+         dbg("fitting");
          GetSurfaceFittingError(x_out_loc, avg_fit_err, max_fit_err);
          if (max_fit_err >= 1.2*init_fit_max_err)
          {
@@ -610,6 +678,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
          }
       }
 
+      // dbg("Serial:{}", serial);
       if (serial)
       {
          energy_out = nlf->GetGridFunctionEnergy(x_out_loc);
@@ -620,6 +689,10 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
          energy_out = p_nlf->GetParGridFunctionEnergy(x_out_loc);
       }
 #endif
+      dbg("energy_out: {}", energy_out);
+      // assert(AlmostEq(energy_out, 0.7791684091777566) && "❌❌❌❌");
+      // assert(false && "✅✅✅✅");
+
       if (energy_out > energy_in + 0.2*fabs(energy_in) ||
           std::isnan(energy_out) != 0)
       {

fem/tmop_tools.hpp

@@ -29,7 +29,7 @@ private:
    Vector field0;
    const real_t dt_scale;
    const AssemblyLevel al;
-   MemoryType opt_mt = MemoryType::DEFAULT;
+   MemoryType opt_mt = MemoryType::PRESERVE;//DEFAULT;
 
    void ComputeAtNewPositionScalar(const Vector &new_nodes, Vector &new_field);
 public:

general/nvtx.hpp

@@ -0,0 +1,481 @@
+#pragma once
+
+#include <fmt/format.h>
+
+#include <array>
+#include <cassert>
+#include <cstddef>
+#include <cstdint>
+#include <iomanip>
+#include <iostream>
+#include <memory>
+#include <mutex>
+#include <stack>
+#include <string>
+#include <utility>
+
+#ifdef MFEM_USE_CALIPER
+#include <caliper/cali.h>
+#endif
+
+#ifdef MFEM_USE_CUDA
+#include <cudaProfiler.h>
+#include <cuda_runtime_api.h>
+#include <nvToolsExt.h>
+#else
+struct nvtxEventAttributes_t
+{
+   int version;
+   int size;
+   int category;
+   int colorType;
+   uint32_t color;
+   int payloadType;
+   uint64_t payload;
+   int messageType;
+   struct
+   {
+      std::string ascii;
+   } message;
+};
+#define NVTX_VERSION 1
+#define NVTX_EVENT_ATTRIB_STRUCT_SIZE 256
+#define NVTX_COLOR_ARGB 0
+#define NVTX_MESSAGE_TYPE_ASCII 0
+#define nvtxRangePushEx(...)
+#define nvtxRangePop(...)
+#define cudaStreamSynchronize(...)
+#endif
+
+//// ///////////////////////////////////////////////////////////////////////////
+template <class T>
+std::enable_if_t<!std::numeric_limits<T>::is_integer, bool>
+AlmostEq(T x, T y, T tolerance = 10.0*std::numeric_limits<T>::epsilon())
+{
+   const T neg = std::abs(x - y);
+   constexpr T min = std::numeric_limits<T>::min();
+   constexpr T eps = std::numeric_limits<T>::epsilon();
+   const T min_abs = std::min(std::abs(x), std::abs(y));
+   if (std::abs(min_abs) == 0.0) { return neg < eps; }
+   return (neg / (1.0 + std::max(min, min_abs))) < tolerance;
+}
+
+namespace gpu::nvtx
+{
+
+///////////////////////////////////////////////////////////////////////////////
+// https://en.wikipedia.org/wiki/Web_colors#Extended_colors
+// http://www.calmar.ws/vim/256-xterm-24bit-rgb-color-chart.html
+// clang-format off
+enum color_names
+{
+   kBlack = 0, kNavyBlue, kDarkBlue, kMediumBlue, kBlue, kDarkGreen, kWebGreen, kTeal,
+   kDarkCyan, kDeepSkyBlue, kDarkTurquoise, kMediumSpringGreen, kGreen, kLime,
+   kSpringGreen, kAqua, kCyan, kMidnightBlue, kDodgerBlue, kLightSeaGreen, kForestGreen,
+   kSeaGreen, kDarkSlateGray, kLimeGreen, kMediumSeaGreen, kTurquoise, kRoyalBlue,
+   kSteelBlue, kDarkSlateBlue, kMediumTurquoise, kIndigo, kDarkOliveGreen, kCadetBlue,
+   kCornflower, kRebeccaPurple, kMediumAquamarine, kDimGray, kSlateBlue, kOliveDrab,
+   kSlateGray, kLightSlateGray, kMediumSlateBlue, kLawnGreen, kWebMaroon, kWebPurple,
+   kChartreuse, kAquamarine, kOlive, kWebGray, kSkyBlue, kLightSkyBlue, kBlueViolet,
+   kDarkRed, kDarkMagenta, kSaddleBrown, kDarkSeaGreen, kLightGreen, kMediumPurple,
+   kDarkViolet, kPaleGreen, kDarkOrchid, kYellowGreen, kPurple, kSienna, kBrown,
+   kDarkGray, kLightBlue, kGreenYellow, kPaleTurquoise, kMaroon, kLightSteelBlue,
+   kPowderBlue, kFirebrick, kDarkGoldenrod, kMediumOrchid, kRosyBrown, kDarkKhaki,
+   kGray, kSilver, kMediumVioletRed, kIndianRed, kPeru, kChocolate, kTan, kLightGray,
+   kThistle, kOrchid, kGoldenrod, kPaleVioletRed, kCrimson, kGainsboro, kPlum, kBurlywood,
+   kLightCyan, kLavender, kDarkSalmon, kViolet, kPaleGoldenrod, kLightCoral, kKhaki,
+   kAliceBlue, kHoneydew, kAzure, kSandyBrown, kWheat, kBeige, kWhiteSmoke, kMintCream,
+   kGhostWhite, kSalmon, kAntiqueWhite, kLinen, kLightGoldenrod, kOldLace, kRed,
+   kFuchsia, kMagenta, kDeepPink, kOrangeRed, kTomato, kHotPink, kCoral, kDarkOrange,
+   kLightSalmon, kOrange, kLightPink, kPink, kGold, kPeachPuff, kNavajoWhite, kMoccasin,
+   kBisque, kMistyRose, kBlanchedAlmond, kPapayaWhip, kLavenderBlush, kSeashell,
+   kCornsilk, kLemonChiffon, kFloralWhite, kSnow, kYellow, kLightYellow, kIvory, kWhite,
+   kNvidia
+};
+// clang-format on
+
+static constexpr int kNumHexColors = 146;
+static constexpr std::array<uint32_t, kNumHexColors> kHexColors =
+{
+   {
+      0x000000, 0x000080, 0x00008B, 0x0000CD, 0x0000FF, 0x006400, 0x008000,
+      0x008080, 0x008B8B, 0x00BFFF, 0x00CED1, 0x00FA9A, 0x00FF00, 0x00FF00,
+      0x00FF7F, 0x00FFFF, 0x00FFFF, 0x191970, 0x1E90FF, 0x20B2AA, 0x228B22,
+      0x2E8B57, 0x2F4F4F, 0x32CD32, 0x3CB371, 0x40E0D0, 0x4169E1, 0x4682B4,
+      0x483D8B, 0x48D1CC, 0x4B0082, 0x556B2F, 0x5F9EA0, 0x6495ED, 0x663399,
+      0x66CDAA, 0x696969, 0x6A5ACD, 0x6B8E23, 0x708090, 0x778899, 0x7B68EE,
+      0x7CFC00, 0x7F0000, 0x7F007F, 0x7FFF00, 0x7FFFD4, 0x808000, 0x808080,
+      0x87CEEB, 0x87CEFA, 0x8A2BE2, 0x8B0000, 0x8B008B, 0x8B4513, 0x8FBC8F,
+      0x90EE90, 0x9370DB, 0x9400D3, 0x98FB98, 0x9932CC, 0x9ACD32, 0xA020F0,
+      0xA0522D, 0xA52A2A, 0xA9A9A9, 0xADD8E6, 0xADFF2F, 0xAFEEEE, 0xB03060,
+      0xB0C4DE, 0xB0E0E6, 0xB22222, 0xB8860B, 0xBA55D3, 0xBC8F8F, 0xBDB76B,
+      0xBEBEBE, 0xC0C0C0, 0xC71585, 0xCD5C5C, 0xCD853F, 0xD2691E, 0xD2B48C,
+      0xD3D3D3, 0xD8BFD8, 0xDA70D6, 0xDAA520, 0xDB7093, 0xDC143C, 0xDCDCDC,
+      0xDDA0DD, 0xDEB887, 0xE0FFFF, 0xE6E6FA, 0xE9967A, 0xEE82EE, 0xEEE8AA,
+      0xF08080, 0xF0E68C, 0xF0F8FF, 0xF0FFF0, 0xF0FFFF, 0xF4A460, 0xF5DEB3,
+      0xF5F5DC, 0xF5F5F5, 0xF5FFFA, 0xF8F8FF, 0xFA8072, 0xFAEBD7, 0xFAF0E6,
+      0xFAFAD2, 0xFDF5E6, 0xFF0000, 0xFF00FF, 0xFF00FF, 0xFF1493, 0xFF4500,
+      0xFF6347, 0xFF69B4, 0xFF7F50, 0xFF8C00, 0xFFA07A, 0xFFA500, 0xFFB6C1,
+      0xFFC0CB, 0xFFD700, 0xFFDAB9, 0xFFDEAD, 0xFFE4B5, 0xFFE4C4, 0xFFE4E1,
+      0xFFEBCD, 0xFFEFD5, 0xFFF0F5, 0xFFF5EE, 0xFFF8DC, 0xFFFACD, 0xFFFAF0,
+      0xFFFAFA, 0xFFFF00, 0xFFFFE0, 0xFFFFF0, 0xFFFFFF, 0x76B900
+   }
+};
+
+///////////////////////////////////////////////////////////////////////////////
+constexpr size_t static_strlen(const char *str)
+{
+   return *str == '\0' ? 0 : static_strlen(str + 1) + 1;
+}
+
+constexpr uint8_t static_checksum8(const char *bfr)
+{
+   unsigned int chk = 0;
+   size_t len = static_strlen(bfr);
+   for (; len; len--, bfr++) { chk += static_cast<unsigned int>(*bfr); }
+   return static_cast<uint8_t>(chk);
+}
+
+constexpr char *static_strrnchr(const char *str, const char c, int n)
+{
+   size_t len = static_strlen(str);
+   char *p = const_cast<char *>(str) + len - 1;
+   for (; n; n--, p--, len--)
+   {
+      for (; len; p--, len--)
+      {
+         if (*p == c) { break; }
+      }
+      if (!len) { return nullptr; }
+      if (n == 1) { return p; }
+   }
+   return nullptr;
+}
+
+inline uint32_t static_color(const uint8_t COLOR, const int RANK,
+                             const char *FILE)
+{
+   constexpr auto kMpiColorShift = 1;
+   const auto rank_shift = kMpiColorShift * RANK;
+   if (COLOR > 0) { return kHexColors[COLOR + rank_shift]; }
+   const auto file_color = static_checksum8(FILE);
+   return kHexColors[(file_color + rank_shift) % kNumHexColors];
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Helpers to generate unique variable names
+#define NVTX_FLF __FILE__, __LINE__, __FUNCTION__
+#define NVTX_PRIVATE_NAME(prefix) NVTX_PRIVATE_CONCAT(prefix, __LINE__)
+#define NVTX_PRIVATE_CONCAT(a, b) NVTX_PRIVATE_CONCAT2(a, b)
+#define NVTX_PRIVATE_CONCAT2(a, b) a##b
+
+#ifndef NVTX_COLOR
+#define NVTX_COLOR ::gpu::nvtx::kBlack
+#endif
+
+///////////////////////////////////////////////////////////////////////////////
+struct Debug
+{
+   const bool debug = false, end = true;
+
+   inline Debug() = default;
+
+   inline Debug(const int RANK, const char *FILE, const int LINE,
+                const char *FUNC, uint8_t COLOR, bool ini = true,
+                bool END = true): debug(true), end(END)
+   {
+      const char *base = static_strrnchr(FILE, '/', 2);
+      const char *file = base ? base + 1 : FILE;
+      const uint32_t rgb = static_color(COLOR, RANK, FILE);
+      const uint8_t r = (rgb >> 16) & 0xFF, g = (rgb >> 8) & 0xFF,
+                    b = rgb & 0xFF;
+      std::cout << "\033[38;2;";
+      std::cout << std::to_string(r) << ";";
+      std::cout << std::to_string(g) << ";";
+      std::cout << std::to_string(b) << "m";
+      if (ini)
+      {
+         std::cout << RANK << std::setw(64) << file << ":";
+         std::cout << "\033[2m" << std::setw(4) << std::left << LINE
+                   << "\033[22m: ";
+         if (FUNC) { std::cout << "[" << FUNC << "] "; }
+      }
+      std::cout << std::right << "\033[1m";
+   }
+
+   inline ~Debug()
+   {
+      if (debug) { std::cout << "\033[m" << (end ? "\n" : "") << std::flush; }
+   }
+
+   template <typename T>
+   inline void operator<<(const T &arg) const noexcept
+   {
+      if (debug) { std::cout << arg; }
+   }
+
+   template <typename T>
+   inline void operator()(const T &arg) const noexcept
+   {
+      if (debug) { this->operator<<(arg); }
+   }
+
+   template <typename... Args>
+   inline void operator()(const char *fmt, Args &&...args) const noexcept
+   {
+      if (debug) { std::cout << fmt::format(fmt, std::forward<Args>(args)...); }
+   }
+
+   inline void operator()() const noexcept {}
+
+   static Debug Set(const char *FILE, const int LINE, const char *FUNC,
+                    uint8_t COLOR, bool INI = true, bool END = true)
+   {
+      static int mpi_rank = 0, dbg_mpi_rank = 0;
+      static bool env_mpi = false, env_dbg = false;
+      if (static bool ini = false; !std::exchange(ini, true))
+      {
+         env_dbg = (::getenv("MFEM_DEBUG") != nullptr);
+         env_mpi = ::getenv("MFEM_DEBUG_MPI") != nullptr;
+         // int mpi_flag = 0;
+         // MPI_Initialized(&mpi_flag);
+         // if (mpi_flag) { MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); }
+         dbg_mpi_rank = atoi(env_mpi ? ::getenv("MFEM_DEBUG_MPI") : "0");
+      }
+      const bool debug = (env_dbg && (!env_mpi || (dbg_mpi_rank == mpi_rank)));
+      return debug ? Debug(mpi_rank, FILE, LINE, FUNC, COLOR, INI, END)
+             : Debug();
+   }
+};
+
+// Debug console traces, unnamed
+#define NVTX_DEBUG(...)                                                  \
+   ::gpu::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR).operator()(__VA_ARGS__)
+
+#define NVTX_DEBUG_NO_INI(...)                                \
+   ::gpu::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR, false, true) \
+      .operator()(__VA_ARGS__)
+
+#define NVTX_DEBUG_APPEND(...)                                 \
+   ::gpu::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR, false, false) \
+      .operator()(__VA_ARGS__)
+
+#define NVTX_DEBUG_NO_END(...)                                \
+   ::gpu::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR, true, false) \
+      .operator()(__VA_ARGS__)
+
+///////////////////////////////////////////////////////////////////////////////
+struct Nvtx
+{
+   const bool nvtx = false, enforce_kernel_sync = false;
+   const char *base, *file;
+   const uint32_t color = kBlack;
+   mutable std::string ascii;
+   mutable nvtxEventAttributes_t event;
+   mutable bool pushed = false;
+
+   inline Nvtx() = default;
+
+   Nvtx(bool enforce_kernel_sync, const char *FILE, const int LINE,
+        const char *FUNC, uint8_t COLOR):
+      nvtx(true), enforce_kernel_sync(enforce_kernel_sync),
+      base(static_strrnchr(FILE, '/', 2)), file(base ? base + 1 : FILE),
+      color(COLOR), ascii(file), event({})
+   {
+      event.version = NVTX_VERSION;
+      event.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+      event.colorType = NVTX_COLOR_ARGB;
+      event.color = static_color(COLOR, 0, FILE);
+      event.messageType = NVTX_MESSAGE_TYPE_ASCII;
+
+      ascii += ":";
+      ascii += std::to_string(LINE);
+      ascii += ":[";
+      ascii += FUNC;
+      ascii += "] ";
+
+      pushed = false;
+   }
+
+   explicit Nvtx(const char *title, uint8_t color = kWheat,
+                 bool enforce_kernel_sync = true):
+      nvtx(true), enforce_kernel_sync(enforce_kernel_sync), color(color),
+      ascii(title), event({})
+   {
+      event.version = NVTX_VERSION;
+      event.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+      event.colorType = NVTX_COLOR_ARGB;
+      event.color = static_color(color, 0, "");
+      event.messageType = NVTX_MESSAGE_TYPE_ASCII;
+      event.message.ascii = ascii.c_str();
+      nvtxRangePushEx(&event);
+      pushed = true;
+   }
+
+   inline void operator()() const
+   {
+      if (!nvtx) { return; }
+      event.message.ascii = ascii.c_str();
+      assert(!pushed);
+      nvtxRangePushEx(&event);
+      pushed = true;
+   }
+
+   template <typename T>
+   inline void operator()(const T &arg) const
+   {
+      if (!nvtx) { return; }
+      this->operator<<(arg);
+      event.message.ascii = ascii.c_str();
+      assert(!pushed);
+      nvtxRangePushEx(&event);
+      pushed = true;
+   }
+
+   template <typename... Args>
+   inline void operator()(fmt::format_string<Args...> fmt,
+                          Args &&...args) const
+   {
+      if (!nvtx) { return; }
+      ascii += fmt::format(fmt, std::forward<Args>(args)...);
+      event.message.ascii = ascii.c_str();
+      assert(!pushed);
+      nvtxRangePushEx(&event);
+      pushed = true;
+   }
+
+   template <typename T>
+   inline void operator<<(const T &arg) const
+   {
+      if (nvtx) { ascii += arg; }
+   }
+
+   inline ~Nvtx()
+   {
+      if (!nvtx) { return; }
+      if (enforce_kernel_sync)
+      {
+         nvtxEventAttributes_t eks = {};
+         eks.version = NVTX_VERSION;
+         eks.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+         eks.category = 0; // user value
+         eks.colorType = NVTX_COLOR_ARGB;
+         eks.messageType = NVTX_MESSAGE_TYPE_ASCII;
+         eks.message.ascii = "!"; // enforce kernel synchronization
+         eks.color = kHexColors[kYellow];
+         nvtxRangePushEx(&eks);
+         cudaStreamSynchronize(nullptr);
+         nvtxRangePop(/*eks*/);
+      }
+      assert(pushed);
+      nvtxRangePop(/*event*/);
+   }
+
+   using nvtx_ptr = std::unique_ptr<Nvtx>;
+   using nvtx_stack_t = std::stack<nvtx_ptr>;
+
+   static nvtx_ptr Set(const char *FILE, const int LINE, const char *FUNC,
+                       uint8_t COLOR)
+   {
+      static bool nvtx = false, eks = false;
+      if (static bool ini = false; !std::exchange(ini, true))
+      {
+         eks = ::getenv("MFEM_EKS") != nullptr;
+         nvtx = ::getenv("MFEM_NVTX") != nullptr;
+         Nvtx force_first_eks("Init EKS", kYellow, true);
+      }
+      return nvtx_ptr(nvtx ? new Nvtx(eks, FILE, LINE, FUNC, COLOR)
+                      : new Nvtx());
+   }
+
+   static nvtx_stack_t &Stack()
+   {
+      auto nvtx_events = []() -> nvtx_stack_t &
+      {
+         static nvtx_stack_t events;
+         return events;
+      };
+      static std::once_flag ready;
+      // one touch to guarantee the object is ready
+      std::call_once(ready, [&] { nvtx_events(); });
+      return nvtx_events();
+   }
+};
+
+// Temporary object only alive for the current statement
+#define NVTX_(COLOR, ...)                                       \
+   NVTX_DEBUG(__VA_ARGS__);                                     \
+   std::unique_ptr<::gpu::nvtx::Nvtx> NVTX_PRIVATE_NAME(nvtx) = \
+      ::gpu::nvtx::Nvtx::Set(NVTX_FLF, COLOR);                  \
+   NVTX_PRIVATE_NAME(nvtx)->operator()(__VA_ARGS__)
+
+// Temporary object only alive for the current statement
+#define NVTX(...) NVTX_(NVTX_COLOR, __VA_ARGS__)
+
+// Begin(with color)/End NVTX event traces
+#define NVTX_BEGIN_(COLOR, ...)                                              \
+   NVTX_DEBUG(__VA_ARGS__);                                                  \
+   ::gpu::nvtx::Nvtx::Stack().push(::gpu::nvtx::Nvtx::Set(NVTX_FLF, COLOR)); \
+   ::gpu::nvtx::Nvtx::Stack().top()->operator()(__VA_ARGS__)
+
+// Begin/End NVTX event traces
+#define NVTX_BEGIN(...) NVTX_BEGIN_(NVTX_COLOR, __VA_ARGS__);
+
+#define NVTX_END(...)                        \
+   ::gpu::nvtx::Nvtx::Stack().top().reset(); \
+   ::gpu::nvtx::Nvtx::Stack().pop()
+
+#ifdef USE_CALIPER
+// CALIPER & NVTX marks
+#define NVTX_MARK_FUNCTION                                            \
+      NVTX();                                                            \
+      std::unique_ptr<cali::Function> __cali_ann##__func__;              \
+      __cali_ann##__func__ = std::make_unique<cali::Function>(__func__);
+
+#define NVTX_MARK(...)                                                   \
+      NVTX(__VA_ARGS__);                                                    \
+      std::unique_ptr<cali::Function> __cali_ann##__func__;                 \
+      __cali_ann##__func__ = std::make_unique<cali::Function>(__VA_ARGS__);
+
+#define NVTX_MARK_FUNCTION_NAME(STR_NAME)                                \
+      NVTX(STR_NAME);                                                       \
+      std::unique_ptr<cali::Function> __cali_ann##__func__;                 \
+      if (g_caliper) {                                                      \
+         __cali_ann##__func__ = std::make_unique<cali::Function>(STR_NAME); \
+      }
+
+#define NVTX_MARK_BEGIN(...)     \
+      CALI_MARK_BEGIN(__VA_ARGS__); \
+      NVTX_BEGIN(__VA_ARGS__);
+
+#define NVTX_MARK_END(...)     \
+      NVTX_END(__VA_ARGS__);      \
+      CALI_MARK_END(__VA_ARGS__);
+#else
+#define NVTX_MARK_FUNCTION NVTX()
+#define NVTX_MARK(...) NVTX(__VA_ARGS__)
+#define NVTX_MARK_FUNCTION_NAME(...) NVTX(__VA_ARGS__)
+#define NVTX_MARK_BEGIN(...) NVTX_BEGIN(__VA_ARGS__)
+#define NVTX_MARK_END(...) NVTX_END(__VA_ARGS__)
+#endif
+
+} // namespace gpu::nvtx
+
+// Debug console traces, unnamed
+#if 1
+#define dbg(...) NVTX_DEBUG(__VA_ARGS__)
+#define dbl(...) NVTX_DEBUG_NO_END(__VA_ARGS__)
+#define dba(...) NVTX_DEBUG_APPEND(__VA_ARGS__)
+#define dbc(...) NVTX_DEBUG_NO_INI(__VA_ARGS__)
+#else
+#define dbg(...)
+#define dbl(...) (void)0
+#define dba(...)
+#define dbc(...)
+#endif
+
+inline bool ClearScreen()
+{
+   dbg("\x1B[2J\x1B[3J\x1B[H");
+   return true;
+}

makefile

@@ -473,6 +473,10 @@ $(MINIAPP_USE_COMMON): miniapps/common
 $(EM_DIRS) $(TEST_DIRS): lib
 	$(MAKE) -C $(BLD)$(@)
 
+.PHONY: mesh-optimizer
+mesh-optimizer: lib
+	$(MAKE) -C miniapps/meshing mesh-optimizer -j 10
+
 .PHONY: doc
 doc:
 	$(MAKE) -C $(BLD)$(@)

miniapps/meshing/pmesh-optimizer.cpp

@@ -114,11 +114,15 @@
 #include <fstream>
 #include "mesh-optimizer.hpp"
 
+#define NVTX_COLOR ::gpu::nvtx::color_names::kYellow
+#include "../../general/nvtx.hpp"
+
 using namespace mfem;
 using namespace std;
 
 int main (int argc, char *argv[])
 {
+   NVTX_MARK_FUNCTION;
    // 0. Initialize MPI and HYPRE.
    Mpi::Init(argc, argv);
    int myid = Mpi::WorldRank();