Generify hydro MOL for anisotropic EB

Exec/RegTests/EB-C7/eb-c7-unequal-dx.inp

@@ -0,0 +1,94 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 100000
+stop_time = 0.1
+
+# PROBLEM SIZE & GEOMETRY
+geometry.is_periodic =  0   0  1
+geometry.coord_sys   =  0       # 0 => cart
+geometry.prob_lo     = 0.0 -0.30000000000000004 0.0
+geometry.prob_hi     = 2.0 1.5 0.1
+amr.n_cell           = 320 144 8
+#amr.n_cell           = 160 144 8
+
+
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+# Interior, UserBC, Symmetry, SlipWall, NoSlipWall
+# >>>>>>>>>>>>>  BC KEYWORDS <<<<<<<<<<<<<<<<<<<<<<
+pelec.lo_bc       =  "FOExtrap"  "NoSlipWall"  "Interior"
+pelec.hi_bc       =  "FOExtrap"  "NoSlipWall"  "Interior"
+
+# Problem setup
+pelec.eb_boundary_T = 300.
+pelec.eb_isothermal = 1
+
+# WHICH PHYSICS
+pelec.do_hydro = 1
+pelec.do_mol = 1
+pelec.do_react = 0
+pelec.allow_negative_energy = 0
+pelec.diffuse_temp = 0
+pelec.diffuse_vel  = 0
+pelec.diffuse_spec = 0
+pelec.diffuse_enth = 0
+
+# TIME STEP CONTROL
+pelec.dt_cutoff      = 5.e-20  # level 0 timestep below which we halt
+pelec.cfl            = 0.3     # cfl number for hyperbolic system
+pelec.init_shrink    = 0.8    # scale back initial timestep
+pelec.change_max     = 1.05     # maximum increase in dt over successive steps
+
+# DIAGNOSTICS & VERBOSITY
+pelec.sum_interval   = 1       # timesteps between computing mass
+pelec.v              = 1       # verbosity in PeleC cpp files
+amr.v                = 1       # verbosity in Amr.cpp
+#amr.grid_log         = grdlog  # name of grid logging file
+
+# REFINEMENT / REGRIDDING
+amr.max_level       = 1       # maximum level number allowed
+amr.ref_ratio       = 2 2 2 2 # refinement ratio
+amr.regrid_int      = 2 2 2 2 # how often to regrid
+amr.blocking_factor = 8       # block factor in grid generation
+amr.max_grid_size   = 32
+
+# CHECKPOINT FILES
+amr.checkpoint_files_output = 0
+amr.check_file      = chk      # root name of checkpoint file
+amr.check_int       = -1       # number of timesteps between checkpoints
+
+# PLOTFILES
+amr.plot_files_output = 1
+amr.plot_file       = plt
+amr.plot_int        = 100
+amr.derive_plot_vars=ALL
+
+# PROBLEM DEFINITION
+prob.pl = 1.0
+prob.rhol = 1.0
+prob.pr = 0.1
+prob.rhor = 0.125
+prob.angle = 30.0
+
+# TAGGING
+tagging.denerr = 3
+tagging.dengrad = 0.01
+tagging.max_denerr_lev = 3
+tagging.max_dengrad_lev = 3
+
+tagging.presserr = 3
+tagging.pressgrad = 0.01
+tagging.max_presserr_lev = 3
+tagging.max_pressgrad_lev = 3
+
+#extruded triangles lets the user create a maximum of 5 triangles 
+#in 2D that will be extruded in the z direction
+#make sure the coordinates are in anti-clockwise direction 
+eb2.geom_type = "extruded_triangles"
+extruded_triangles.num_tri = 2
+extruded_triangles.tri_0_point_0 = -20.0 -11.662475437630437 0.0
+extruded_triangles.tri_0_point_1 = 22.0 -20.0 0.0
+extruded_triangles.tri_0_point_2 = 22.0 12.586235868333839 0.0
+
+extruded_triangles.tri_1_point_0 = -20.0 -11.431535329954588 0.0
+extruded_triangles.tri_1_point_1 = 22.0 12.817175976009688 0.0
+extruded_triangles.tri_1_point_2 = -20.0 20.0 0.0
+ebd.boundary_grad_stencil_type = 1

Source/EB.H

@@ -216,7 +216,7 @@ void pc_fill_sv_ebg(
 
 void pc_fill_bndry_grad_stencil_ls(
   const amrex::Box& /*bx*/,
-  const amrex::Real /*dx*/,
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& /*dx*/,
   const int /*unused*/,
   const EBBndryGeom* /*ebg*/,
   const int /*Nsten*/,
@@ -231,7 +231,7 @@ void pc_check_bndry_grad_stencil(
 
 void pc_fill_bndry_grad_stencil_quadratic(
   const amrex::Box& /*bx*/,
-  const amrex::Real /*dx*/,
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& /*dx*/,
   const int /*unused*/,
   const EBBndryGeom* /*ebg*/,
   const int /*Nsten*/,

Source/EB.cpp

@@ -66,14 +66,24 @@ pc_fill_sv_ebg(
 void
 pc_fill_bndry_grad_stencil_quadratic(
   const amrex::Box& bx,
-  const amrex::Real dx,
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx,
   const int /*Nebg*/,
   const EBBndryGeom* ebg,
   const int Nsten,
   EBBndrySten* grad_stencil)
 {
-  const amrex::Real area = std::pow(dx, AMREX_SPACEDIM - 1);
-  const amrex::Real fac = area / dx;
+  const amrex::Real small = 1.e-13;
+  const bool unequal_dx =
+    amrex::max<amrex::Real>(AMREX_D_DECL(
+      static_cast<amrex::Real>(0.0),
+      static_cast<amrex::Real>(std::abs(dx[0] - dx[1])),
+      static_cast<amrex::Real>(std::abs(dx[0] - dx[2])))) > small * dx[0];
+  if (unequal_dx) {
+    amrex::Abort(
+      "dx != dy != dz not supported with pc_fill_bndry_grad_stencil_quadratic");
+  }
+  const amrex::Real area = std::pow(dx[0], AMREX_SPACEDIM - 1);
+  const amrex::Real fac = area / dx[0];
 
   amrex::ParallelFor(Nsten, [=] AMREX_GPU_DEVICE(int L) {
     if (bx.contains(ebg[L].iv)) {
@@ -253,7 +263,7 @@ pc_fill_bndry_grad_stencil_quadratic(
 void
 pc_fill_bndry_grad_stencil_ls(
   const amrex::Box& bx,
-  const amrex::Real dx,
+  const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& dx,
   const int /*Nebg*/,
   const EBBndryGeom* ebg,
   const int Nsten,
@@ -263,8 +273,8 @@ pc_fill_bndry_grad_stencil_ls(
 
   AMREX_ASSERT(AMREX_SPACEDIM > 1);
 
-  const amrex::Real area = std::pow(dx, AMREX_SPACEDIM - 1);
-  const amrex::Real fac = area / dx;
+  const amrex::Real area = std::pow(dx[0], AMREX_SPACEDIM - 1);
+  const amrex::Real fac = area / dx[0];
 
   amrex::ParallelFor(Nsten, [=] AMREX_GPU_DEVICE(int L) {
     if (bx.contains(ebg[L].iv)) {

Source/InitEB.cpp

@@ -132,7 +132,7 @@ PeleC::initialize_eb2_structs()
 
       // Fill in boundary gradient for cut cells in this grown tile
       sv_eb_bndry_grad_stencil[iLocal].resize(ncutcells);
-      const amrex::Real dx = geom.CellSize()[0];
+      const auto& dx = geom.CellSizeArray();
       if (bgs == 0) {
         pc_fill_bndry_grad_stencil_quadratic(
           tbox, dx, ncutcells, sv_eb_bndry_geom[iLocal].data(), ncutcells,
@@ -562,12 +562,13 @@ PeleC::extend_signed_distance(
   // signed distance and propagates it manually up to the point where we need to
   // have it for derefining.
   BL_PROFILE("PeleC::extend_signed_distance()");
-  const auto geomdata = parent->Geom(0).data();
   amrex::Real maxSignedDist = signDist->max(0);
   const auto& ebfactory =
     dynamic_cast<amrex::EBFArrayBoxFactory const&>(signDist->Factory());
   const auto& flags = ebfactory.getMultiEBCellFlagFab();
   int nGrowFac = flags.nGrow() + 1;
+  const auto& dx = parent->Geom(0).CellSizeArray();
+  const amrex::Real dx_max = *std::max_element(dx.begin(), dx.end());
 
   // First set the region far away at the max value we need
   auto const& sd_ccs = signDist->arrays();
@@ -577,8 +578,7 @@ PeleC::extend_signed_distance(
     [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
       const auto& sd_cc = sd_ccs[nbx];
       if (sd_cc(i, j, k) >= maxSignedDist - 1e-12) {
-        const amrex::Real* dx = geomdata.CellSize();
-        sd_cc(i, j, k) = nGrowFac * dx[0] * extendFactor;
+        sd_cc(i, j, k) = nGrowFac * dx_max * extendFactor;
       }
     });
   amrex::Gpu::synchronize();
@@ -602,8 +602,7 @@ PeleC::extend_signed_distance(
       ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
         const auto glo = amrex::lbound(gbx);
         const auto ghi = amrex::ubound(gbx);
-        const amrex::Real* dx = geomdata.CellSize();
-        const amrex::Real extendedDist = dx[0] * extendFactor;
+        const amrex::Real extendedDist = dx_max * extendFactor;
         if (sd_cc(i, j, k) >= maxSignedDist - 1e-12) {
           amrex::Real closestEBDist = 1e12;
           for (int kk = glo.z; kk <= ghi.z; ++kk) {

Source/MOL.cpp

@@ -206,7 +206,6 @@ pc_compute_hyp_mol_flux_eb(
   const int nextra = 0;
   const int bc_test_val = 1;
 
-  const amrex::Real full_area = AMREX_D_PICK(1.0, dx[0], dx[0] * dx[1]);
   const amrex::Box bxg = amrex::grow(cbox, nextra - 1);
   const auto geomdata = geom.data();
   ProbParmDevice const* prob_parm = PeleC::d_prob_parm_device;
@@ -216,26 +215,30 @@ pc_compute_hyp_mol_flux_eb(
     if (bxg.contains(iv)) {
       amrex::Real ebnorm[AMREX_SPACEDIM] = {AMREX_D_DECL(
         ebg[L].eb_normal[0], ebg[L].eb_normal[1], ebg[L].eb_normal[2])};
-      const amrex::Real ebnorm_mag = std::sqrt(AMREX_D_TERM(
-        ebnorm[0] * ebnorm[0], +ebnorm[1] * ebnorm[1], +ebnorm[2] * ebnorm[2]));
-      for (amrex::Real& dir : ebnorm) {
-        dir /= ebnorm_mag;
+
+      amrex::Real ebnorm_local[AMREX_SPACEDIM] = {
+        AMREX_D_DECL(ebnorm[0] / dx[0], ebnorm[1] / dx[1], ebnorm[2] / dx[2])};
+      const amrex::Real ebnorm_mag_local = std::sqrt(AMREX_D_TERM(
+        ebnorm_local[0] * ebnorm_local[0], +ebnorm_local[1] * ebnorm_local[1],
+        +ebnorm_local[2] * ebnorm_local[2]));
+      for (amrex::Real& dir : ebnorm_local) {
+        dir /= ebnorm_mag_local;
       }
 
       amrex::Real flux_tmp[NVAR] = {0.0};
 
       if (!eb_problem_state) {
         AMREX_D_TERM(
-          flux_tmp[UMX] = -q(iv, QPRES) * ebnorm[0];
-          , flux_tmp[UMY] = -q(iv, QPRES) * ebnorm[1];
-          , flux_tmp[UMZ] = -q(iv, QPRES) * ebnorm[2];)
+          flux_tmp[UMX] = -q(iv, QPRES) * ebnorm_local[0];
+          , flux_tmp[UMY] = -q(iv, QPRES) * ebnorm_local[1];
+          , flux_tmp[UMZ] = -q(iv, QPRES) * ebnorm_local[2];)
       } else {
         auto eos = pele::physics::PhysicsType::eos();
         amrex::Real qtempl[R_NUM] = {0.0};
         amrex::Real* spl = &qtempl[R_Y];
         qtempl[R_UN] = -(AMREX_D_TERM(
-          q(iv, QU) * ebnorm[0], +q(iv, QV) * ebnorm[1],
-          +q(iv, QW) * ebnorm[2]));
+          q(iv, QU) * ebnorm_local[0], +q(iv, QV) * ebnorm_local[1],
+          +q(iv, QW) * ebnorm_local[2]));
         qtempl[R_UT1] = 0.0;
         qtempl[R_UT2] = 0.0;
         qtempl[R_P] = q(iv, QPRES);
@@ -276,8 +279,8 @@ pc_compute_hyp_mol_flux_eb(
 
         const bool do_ebfill = ProblemSpecificFunctions::problem_eb_state(
           geomdata, vfrac(iv), iv,
-          AMREX_D_DECL(ebnorm[0], ebnorm[1], ebnorm[2]), qtempl, spl, rhoe_l,
-          gamc_l, qtempr, spr, rhoe_r, gamc_r, prob_parm);
+          AMREX_D_DECL(ebnorm_local[0], ebnorm_local[1], ebnorm_local[2]),
+          qtempl, spl, rhoe_l, gamc_l, qtempr, spr, rhoe_r, gamc_r, prob_parm);
 
         if (do_ebfill) {
 
@@ -333,20 +336,32 @@ pc_compute_hyp_mol_flux_eb(
 
           const amrex::Real tmp_flx_umx = flux_tmp[UMX];
           AMREX_D_TERM(
-            flux_tmp[UMX] = -tmp_flx_umx * ebnorm[0];
-            , flux_tmp[UMY] = -tmp_flx_umx * ebnorm[1];
-            , flux_tmp[UMZ] = -tmp_flx_umx * ebnorm[2];)
+            flux_tmp[UMX] = -tmp_flx_umx * ebnorm_local[0];
+            , flux_tmp[UMY] = -tmp_flx_umx * ebnorm_local[1];
+            , flux_tmp[UMZ] = -tmp_flx_umx * ebnorm_local[2];)
 
         } else {
           AMREX_D_TERM(
-            flux_tmp[UMX] = -q(iv, QPRES) * ebnorm[0];
-            , flux_tmp[UMY] = -q(iv, QPRES) * ebnorm[1];
-            , flux_tmp[UMZ] = -q(iv, QPRES) * ebnorm[2];)
+            flux_tmp[UMX] = -q(iv, QPRES) * ebnorm_local[0];
+            , flux_tmp[UMY] = -q(iv, QPRES) * ebnorm_local[1];
+            , flux_tmp[UMZ] = -q(iv, QPRES) * ebnorm_local[2];)
         }
       }
 
       // Copy result into ebflux vector. Being a bit chicken here and only
       // copy values where ebg % iv is within box
+#if AMREX_SPACEDIM == 1
+      const amrex::Real full_area = 1.0;
+#elif AMREX_SPACEDIM == 2
+      const amrex::Real full_area = std::sqrt(
+        (ebnorm[0] * ebnorm[0]) * dx[1] * dx[1]
+        +(ebnorm[1] * ebnorm[1]) * dx[0] * dx[0]);
+#elif AMREX_SPACEDIM == 3
+      const amrex::Real full_area = std::sqrt(
+        (ebnorm[0] * ebnorm[0]) * (dx[1] * dx[2]) * (dx[1] * dx[2])
+        +(ebnorm[1] * ebnorm[1]) * (dx[0] * dx[2]) * (dx[0] * dx[2])
+        +(ebnorm[2] * ebnorm[2]) * (dx[0] * dx[1]) * (dx[0] * dx[1]));
+#endif
       for (int n = 0; n < NVAR; n++) {
         ebflux[n * nebflux + L] += flux_tmp[n] * ebg[L].eb_area * full_area;
       }

Source/PeleC.cpp

@@ -1798,17 +1798,18 @@ PeleC::errorEst(
 
     // Estimate how far I need to derefine
     const amrex::Real safetyFac = tagging_parm->detag_eb_factor;
-    amrex::Real clearTagDist =
-      parent->Geom(tagging_parm->max_eb_refine_lev).CellSize(0) *
-      static_cast<amrex::Real>(
-        parent->nErrorBuf(tagging_parm->max_eb_refine_lev)) *
-      safetyFac;
+    const auto& dx =
+      parent->Geom(tagging_parm->max_eb_refine_lev).CellSizeArray();
+    const amrex::Real dx_max = *std::max_element(dx.begin(), dx.end());
+    amrex::Real clearTagDist = dx_max *
+                               static_cast<amrex::Real>(parent->nErrorBuf(
+                                 tagging_parm->max_eb_refine_lev)) *
+                               safetyFac;
     const int finest_level = parent->finestLevel();
     for (int ilev = tagging_parm->max_eb_refine_lev + 1; ilev <= finest_level;
          ++ilev) {
       clearTagDist +=
-        static_cast<amrex::Real>(parent->nErrorBuf(ilev)) *
-        parent->Geom(tagging_parm->max_eb_refine_lev).CellSize(0) * safetyFac;
+        static_cast<amrex::Real>(parent->nErrorBuf(ilev)) * dx_max * safetyFac;
     }
 
     // Untag cells too close to EB