EulerOperator_impl.hpp

Go to the documentation of this file.

// Copyright (c) 2025-2026 Board of Trustees of the University of Illinois
//
// This file is part of Theseus.
//
// SPDX-License-Identifier: BSD-3-Clause
 
namespace Theseus
{
  // This is a Inf/NAN detector helper
  int CBE(const mfem::Vector &v)
  {
    mfem::Vector one_bad(v.Size());
    const mfem::real_t *vd = v.Read();
    mfem::real_t *bd_w = one_bad.Write();
 
    mfem::forall(v.Size(), [=] MFEM_HOST_DEVICE (int i) {
      bd_w[i] = Theseus::Kernels::is_bad_value(vd[i]) ? 1.0 : 0.0;
    });
    const mfem::real_t *bd_r = one_bad.Read();
    int nbad = 0;
    for(int i = 0;i < v.Size();i++){
      nbad += bd_r[i];
    }
    return nbad;
  }
 
  // Assemble the Cartesian and optional axisymmetric volume RHS for all elements.
  template<typename PhysicsT>
  void EulerOperator<PhysicsT>::MultEuler_Volume(const mfem::Vector &pu, mfem::Vector &pdudt) const
  {
    Theseus::ScopedTimer timer("MultEuler_Volume");
 
    // This block is executed by the host
    int nval_restr = operator_cache.restr_v->Height();
    if(operator_cache.uVol.Size() != nval_restr){
      operator_cache.uVol.SetSize(nval_restr);
      operator_cache.uVol.UseDevice();
    }
    mfem::Vector &Ue(operator_cache.uVol);
    if(!operator_cache.u_vol_restr_ready){
      Theseus::ScopedTimer timer("VolRestrict");
      operator_cache.restr_v->Mult(pu, Ue);
      operator_cache.u_vol_restr_ready = true;
    }
    if(operator_cache.rhsVol.Size() != nval_restr){
      operator_cache.rhsVol.SetSize(nval_restr);
      operator_cache.rhsVol.UseDevice();
    }
    mfem::Vector &dUe(operator_cache.rhsVol);
    // Zero the array on-device (*needed for now*)
    {
      Theseus::ScopedTimer zerotim("ZeroRHSStorage");
      mfem::real_t *d = dUe.Write();
      mfem::forall(dUe.Size(), [=] MFEM_HOST_DEVICE (int i) { d[i] = mfem::real_t(0); });
    }
 
    const mfem::real_t *Ue_d = Ue.Read();
    mfem::real_t *dUe_d = dUe.Write();
 
    // Copy the device cache so that it is not member data
    auto dc = device_cache;
 
    // Device cache parameters
    const int dim = dc.dim;
    const int ne = dc.num_elements;
    const int ndof = dc.ndof_scalar_el;
    const int neq = dc.num_equations;
#ifdef SUBCELL_FV_BLENDING
    const int Np_x = dc.Np_x;
    const int Np_y = dc.Np_y;
    const int Np_z = dc.Np_z;
    const int npe_metric_xi = (Np_x + 1)*Np_y*Np_z;
    const int npe_metric_eta = Np_x*(Np_y + 1)*Np_z;
    const int npe_metric_zeta = Np_x * Np_y * (Np_z + 1);
    const mfem::real_t *metric_xi_d = dc.subcell_metric_xi_d;
    const mfem::real_t *metric_eta_d = (dim > 1 ? dc.subcell_metric_eta_d : nullptr);
    const mfem::real_t *metric_zeta_d = (dim > 2 ? dc.subcell_metric_zeta_d : nullptr);
 
    mfem::Vector &dUfv(operator_cache.dUfv);
    if(dUfv.Size() == 0){
      dUfv.SetSize(nval_restr);
      dUfv.UseDevice();
      operator_cache.dUfv_d = dUfv.Write();
    }
    mfem::real_t *dUfv_d = operator_cache.dUfv_d;  // zeroing not needed, kernel does it
 
    const mfem::real_t *alpha_d = operator_cache.alpha->Read();
 
#endif
    // Derived parameters
    const int metric_stride = ndof * dim * dim;
    const int jac_stride    = ndof;
    const int estride = ndof*neq;
 
    // Device cache data/arrays
    const int *elem_attr_d = dc.elem_attr_d;
    const int *attr_marker_d = dc.attr_marker_d;
    const mfem::real_t *elJac_d = dc.elJac_d;
    const mfem::real_t *elMetric_d = dc.elMetric_d;
    const mfem::real_t *elRadius_d = dc.elRadius_d;
 
    // Inside the FORALL below, executed on device
#ifdef POINT_PARALLEL_VOLUME
    const int npoints = ne * ndof;
    mfem::forall(npoints, [=] MFEM_HOST_DEVICE (int p)
    {
      const int e = p / ndof;
      const int point = p % ndof;
      const int attr = elem_attr_d[e];
      if (attr_marker_d[attr-1] == 0) {
        return;
      }
 
      const int element_offset = e * estride;
      DGSEMIntegrator::AssembleVolumePointKernel(
        dc, Ue_d + element_offset, elJac_d + e*jac_stride,
        elMetric_d + e*metric_stride, point, dUe_d + element_offset);
    });
 
    mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
    {
      const int attr = elem_attr_d[e];
      if (attr_marker_d[attr-1] == 0) {
        return;
      }
 
      const int element_offset = e * estride;
      const mfem::real_t *u_el = Ue_d + element_offset;
      mfem::real_t *du_el = dUe_d + element_offset;
      const mfem::real_t *jac_el = elJac_d + e*jac_stride;
      const mfem::real_t *metric_el = elMetric_d + e*metric_stride;
      const mfem::real_t *radius_el =
        dc.axisymmetric ? elRadius_d + e*jac_stride : nullptr;
#ifdef SUBCELL_FV_BLENDING
      const mfem::real_t alpha_fv = alpha_d[e];
      if (alpha_fv > 1e-16) {
        const mfem::real_t alpha_dg = 1.0 - alpha_fv;
        mfem::real_t *du_fv = dUfv_d + element_offset;
        const mfem::real_t *el_metric_xi =
          metric_xi_d + e*npe_metric_xi*dim;
        const mfem::real_t *el_metric_eta = dim > 1 ?
          metric_eta_d + e*npe_metric_eta*dim : nullptr;
        const mfem::real_t *el_metric_zeta = dim > 2 ?
          metric_zeta_d + e*npe_metric_zeta*dim : nullptr;
        DGSEMIntegrator::ComputeFVFluxesKernel(
          dc, u_el, jac_el, el_metric_xi, el_metric_eta,
          el_metric_zeta, du_fv);
        for (int value = 0; value < estride; ++value) {
          du_el[value] =
            alpha_dg*du_el[value] + alpha_fv*du_fv[value];
        }
      }
#endif
 
      AddAxisymmetricEulerElementSource(
                                        dc, u_el, radius_el, jac_el, metric_el, du_el);
    });
#else
    mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
    {
 
      const mfem::real_t *jac_el    = elJac_d    + e * jac_stride;
      const mfem::real_t *metric_el = elMetric_d + e * metric_stride;
      const mfem::real_t *radius_el = dc.axisymmetric ?
        elRadius_d + e * jac_stride : nullptr;
 
      const int attr = elem_attr_d[e];
      if (attr_marker_d[attr-1] == 0) {
        return;
      }
 
      const int eoff = e * estride;
      const mfem::real_t *u_el = Ue_d + eoff;
      mfem::real_t *du_el = dUe_d + eoff;
 
      // Additive on du_el (*zero first if needed*)
      DGSEMIntegrator::AssembleElementVolumeKernel(dc, u_el,
                                                   jac_el, metric_el, du_el);
#ifdef SUBCELL_FV_BLENDING
      mfem::real_t alpha_fv = alpha_d[e];
      if(alpha_fv > 1e-16){
        mfem::real_t alpha_inv = (1.0 - alpha_fv);
        mfem::real_t *du_fv = dUfv_d + eoff;
        const mfem::real_t *el_metric_xi = metric_xi_d + e * npe_metric_xi * dim;
        const mfem::real_t *el_metric_eta = (dim > 1 ? metric_eta_d + e * npe_metric_eta * dim :
                                             nullptr);
        const mfem::real_t *el_metric_zeta = (dim > 2 ? metric_zeta_d + e * npe_metric_zeta * dim :
                                              nullptr);
        DGSEMIntegrator::ComputeFVFluxesKernel(
          dc, u_el, jac_el, el_metric_xi, el_metric_eta,
          el_metric_zeta, du_fv);
 
        for(int ipt = 0;ipt < estride;ipt++){
          du_el[ipt] = alpha_inv * du_el[ipt] + alpha_fv * du_fv[ipt];
        }
 
      }
#endif
 
      AddAxisymmetricEulerElementSource(
                                        dc, u_el, radius_el, jac_el, metric_el, du_el);
 
    });
#endif
 
    // Scatter RHS back to storage
    operator_cache.restr_v->AddMultTranspose(dUe, pdudt);
 
  }
 
  template<typename PhysicsT>
  void EulerOperator<PhysicsT>::MultEuler_InteriorFaces(const mfem::Vector &pu, mfem::Vector &pdudt) const
  {
    Theseus::ScopedTimer timer("MultEuler_InteriorFaces");
    auto dc = device_cache;
    const int dim = dc.dim;
    const int neq = dc.num_equations;
    const int nfp = dc.num_face_points;
    const int nval_restr = operator_cache.restr_f->Height();
    const int nfaces = nval_restr / (nfp * neq * 2); // (+/-)
    const int npoints = nfaces * nfp;
    const int face_size = 2*nfp*neq;
 
    if(operator_cache.uInt.Size() != nval_restr){
      operator_cache.uInt.SetSize(nval_restr);
      operator_cache.uInt.UseDevice(true);
    }
    mfem::Vector &u_faces(operator_cache.uInt);
    if(!operator_cache.u_int_restr_ready){
      operator_cache.restr_f->Mult(pu, u_faces);
      operator_cache.u_int_restr_ready = true;
    }
 
    if(operator_cache.rhsInt.Size() != nval_restr){
      operator_cache.rhsInt.SetSize(nval_restr);
      operator_cache.rhsInt.UseDevice(true);
    }
    mfem::Vector &rhs_faces(operator_cache.rhsInt);
 
    // For now, just keep this copy - i think it is device-friendly
    mfem::Vector faces_dudt(pdudt);
    // Zeroing unneeded: we clobber anything there on assignment
    faces_dudt.UseDevice(true);
 
    const mfem::real_t *u_d = u_faces.Read();
    mfem::real_t *rhs_d = rhs_faces.Write();
 
    const mfem::real_t *nor_d   = dc.nor_d;      // size nfaces*nfp*dim
    const mfem::real_t *inv1_d  = dc.fw_minus_d; // size nfaces*nfp
    const mfem::real_t *inv2_d  = dc.fw_plus_d;  // size nfaces*nfp
 
#ifdef POINT_PARALLEL_INTERIOR_FACES
    mfem::forall(npoints, [=] MFEM_HOST_DEVICE (int p)
    {
      const int f = p / nfp;
      const int fp = p % nfp;
      const int face_offset = f*face_size;
      const int point_offset = f*nfp + fp;
 
      const mfem::real_t *u_face_d = u_d + face_offset;
      mfem::real_t *rhs_face_d = rhs_d + face_offset;
 
      DGSEMIntegrator::AssembleFacePointKernel(
        dc, u_face_d, nor_d + point_offset*dim,
        inv1_d[point_offset], inv2_d[point_offset], fp, rhs_face_d);
    });
#else
    const int norm_size = nfp*dim;
    mfem::forall(nfaces, [=] MFEM_HOST_DEVICE (int f)
    {
      const int face_offset = f*face_size;
      const int n_offset = f*norm_size;
      const int w_offset = f*nfp;
 
      const mfem::real_t *u_face_d = u_d + face_offset;
      mfem::real_t *rhs_face_d = rhs_d + face_offset;
      const mfem::real_t *nor_face_d = nor_d + n_offset;
      const mfem::real_t *w_minus_d = inv1_d + w_offset;
      const mfem::real_t *w_plus_d = inv2_d + w_offset;
 
      DGSEMIntegrator::AssembleElementFaceKernel(
        dc, u_face_d, nor_face_d, w_minus_d, w_plus_d, rhs_face_d);
 
    });
#endif
 
    operator_cache.restr_f->MultTranspose(rhs_faces, faces_dudt);
    pdudt += faces_dudt; // on device?
 
  }
 
  template<typename PhysicsT>
  void EulerOperator<PhysicsT>::MultEuler_BoundaryFaces(const mfem::Vector &pu, mfem::Vector &pdudt) const
  {
    Theseus::ScopedTimer timer("MultEuler_BoundaryFaces");
 
    auto dc = device_cache;
    const int dim = dc.dim;
    const int neq = dc.num_equations;
    const int nfp = dc.num_face_points;
    const int face_size = nfp * neq;
    const int restr_size = operator_cache.restr_b->Height();
    const int nfaces_restr = restr_size / face_size;
    const int norm_size = nfp * dc.dim;
    const int npoints_bnd = nfaces_restr * nfp;
 
    if(restr_size == 0){
      return;
    }
 
    if(operator_cache.uBnd.Size() != restr_size)
      {
        operator_cache.uBnd.SetSize(restr_size);
        operator_cache.uBnd.UseDevice();
      }
    mfem::Vector &u_faces(operator_cache.uBnd);
    if(!operator_cache.u_bnd_restr_ready){
      operator_cache.restr_b->Mult(pu, u_faces);
      operator_cache.u_bnd_restr_ready = true;
    }
    if(operator_cache.rhsBnd.Size() != restr_size){
      operator_cache.rhsBnd.SetSize(restr_size);
      operator_cache.rhsBnd.UseDevice(true);
    }
    if(operator_cache.dudtBnd.Size() != restr_size){
      operator_cache.dudtBnd.SetSize(pdudt.Size());
      operator_cache.dudtBnd.UseDevice(true);
    }
 
    mfem::Vector &rhs_faces(operator_cache.rhsBnd);
    mfem::Vector faces_dudt(pdudt);
    faces_dudt.UseDevice(true);
 
    // Zero on device:
    {
      mfem::real_t *rd = rhs_faces.Write();
      mfem::forall(rhs_faces.Size(), [=] MFEM_HOST_DEVICE (int i)
      { rd[i] = mfem::real_t(0);});
    }
 
 
    const mfem::real_t *u_d = u_faces.Read();
    mfem::real_t *rhs_d = rhs_faces.Write();
 
    const mfem::real_t *nor_d   = dc.bnd_nor_d;      // size nfaces*nfp*dim
    const mfem::real_t *inv1_d  = dc.bnd_wt_d; // size nfaces*nfp
    const int *bnd_marker_index_d = dc.bnd_marker_index_d;
    mfem::forall(npoints_bnd, [=] MFEM_HOST_DEVICE (int p)
    {
      const int f = p / nfp;
      const int fp = p % nfp;
 
      int bnd_face_marker_index = bnd_marker_index_d[f];
      if(bnd_face_marker_index < 0){
        return;
      }
      //    int bc_index = bnd_marker_to_bc_descr_d[bnd_face_marker_index];
      int bc_index = bnd_face_marker_index; // no mapping atm
      if(bc_index < 0){
        return;
      }
      const Theseus::BCDescriptor &bc = dc.bc_descr_d[bc_index];
      if (bc.type == int(Theseus::BCType::Invalid))
        {
          return;
        }
 
      const int face_offset = f * face_size;
      const int n_offset = f * norm_size;
      const int w_offset = f * nfp;
 
      const mfem::real_t *u_face_d = u_d + face_offset;
      mfem::real_t *rhs_face_d = rhs_d + face_offset;
      const mfem::real_t *nor_face_d = nor_d + n_offset;
      const mfem::real_t *w_minus_d = inv1_d + w_offset;
      const mfem::real_t *nor_point = nor_face_d + fp*dim;
      mfem::real_t scale = -w_minus_d[fp];
      mfem::real_t state1[Theseus::MAXEQ];
      mfem::real_t fluxN[Theseus::MAXEQ];
 
      Theseus::Kernels::el_gather_state(u_face_d, nfp, neq, fp, state1);
      Theseus::BC::ApplyBoundaryConditionInviscid(dc, bc, state1,
                                                  nor_point, fluxN);
      Theseus::Kernels::el_scatter_add(fluxN, nfp, neq, fp, scale, rhs_face_d);
 
    });
 
    operator_cache.restr_b->MultTranspose(rhs_faces, faces_dudt);
 
    pdudt += faces_dudt; // on device? (likely yes)
 
  }
 
 
  // Top level MULT for inviscid cases, called from DGSEMOperator
  template<typename PhysicsT>
  void EulerOperator<PhysicsT>::FlowMult(const mfem::Vector &u, mfem::Vector &pdudt) const
  {
    Theseus::ScopedTimer timer("EulerMult");
 
    auto report_bad = [&](const char *name, const mfem::Vector &v)
    {
      int nbad = CBE(v);
      if (nbad)
        {
          mfem::out << "BAD VALUES IN: (" << name << "), count=" << nbad << std::endl;
        }
    };
 
    const mfem::Vector &pu(this->Prolongate(u));
 
    // This step overwrites contents of pdudt
    MultEuler_Volume(pu, pdudt);
 
    MultEuler_InteriorFaces(pu, pdudt);
    // report_bad("int rhs", pdudt);
 
    MultEuler_BoundaryFaces(pu, pdudt);
    // report_bad("bnd rhs", pdudt);
 
  }
 
}