Setup EB Scalar DevTest.

Exec/DevTests/EB_Test/ERF_Prob.H

@@ -23,6 +23,9 @@ struct ProbParm : ProbParmDefaults {
   amrex::Real yc_frac = 0.5; // Location of "center" of scalar (multiplies domain length)
   amrex::Real zc_frac = 0.5; // Location of "center" of scalar (multiplies domain length)
 
+  amrex::Real xradius = 10.0; // x-radius of scalar bubble
+  amrex::Real zradius = 10.0; // z-radius of scalar bubble
+
   int prob_type = -1;
 }; // namespace ProbParm
 

Exec/DevTests/EB_Test/ERF_Prob.cpp

@@ -32,6 +32,9 @@ Problem::Problem()
 
     pp.query("prob_type", parms.prob_type);
 
+    pp.query("xradius", parms.xradius);
+    pp.query("zradius", parms.zradius);
+
     init_base_parms(parms.rho_0, parms.T_0);
 }
 
@@ -81,16 +84,71 @@ Problem::init_custom_pert(
 
     AMREX_ALWAYS_ASSERT(bx.length()[2] == khi+1);
 
-    ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    // ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
+    // {
+    //     // Set scalar = 0 everywhere
+    //     state_pert(i, j, k, RhoScalar_comp) = 0.0;
+
+    //     if (use_moisture) {
+    //         state_pert(i, j, k, RhoQ1_comp) = 0.0;
+    //         state_pert(i, j, k, RhoQ2_comp) = 0.0;
+    //     }
+    //   });
+
+    // Set the state_pert
+    ParallelFor(bx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
     {
-        // Set scalar = 0 everywhere
-        state_pert(i, j, k, RhoScalar_comp) = 0.0;
+        // Geometry
+        const Real* prob_lo = geomdata.ProbLo();
+        const Real* prob_hi = geomdata.ProbHi();
+        const Real* dx = geomdata.CellSize();
+        const Real x = prob_lo[0] + (i + 0.5) * dx[0];
+        const Real y = prob_lo[1] + (j + 0.5) * dx[1];
+        const Real z = prob_lo[2] + (k + 0.5) * dx[2];
+
+        // Define a point (xc,yc,zc) at the center of the domain
+        const Real xc = parms_d.xc_frac * (prob_lo[0] + prob_hi[0]);
+        const Real yc = parms_d.yc_frac * (prob_lo[1] + prob_hi[1]);
+        const Real zc = parms_d.zc_frac * (prob_lo[2] + prob_hi[2]);
+
+        // Define ellipse parameters
+        const Real r0   = parms_d.rad_0 * (prob_hi[0] - prob_lo[0]);
+        const Real r3d    = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc) + (z-zc)*(z-zc));
+        const Real r2d_xy = std::sqrt((x-xc)*(x-xc) + (y-yc)*(y-yc));
+        const Real r2d_xz = std::sqrt((x-xc)*(x-xc) + (z-zc)*(z-zc));
+        const Real r2d_xz_nd   = std::sqrt((x-xc)*(x-xc)/parms_d.xradius/parms_d.xradius
+                                        + (z-zc)*(z-zc)/parms_d.zradius/parms_d.zradius);
+
+        if (parms_d.prob_type == 10)
+        {
+            // Set scalar = A_0*exp(-10r^2), where r is distance from center of domain,
+            //            + B_0*sin(x)
+            // state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0 * exp(-10.*r3d*r3d) + parms_d.B_0*sin(x);
+            state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0 * exp(-0.1*r2d_xz*r2d_xz) + parms_d.B_0*sin(x);
+
+        } else if (parms_d.prob_type == 11) {
+            if (r2d_xz_nd < 1.0)
+            {
+                state_pert(i, j, k, RhoScalar_comp) = 0.5 * parms_d.A_0 * (1.0 + std::cos(PI*r2d_xz_nd));
+            } else {
+                state_pert(i, j, k, RhoScalar_comp) = 0.0;
+            }
+        } else {
+            // Set scalar = A_0 in a ball of radius r0 and 0 elsewhere
+            if (r3d < r0) {
+            state_pert(i, j, k, RhoScalar_comp) = parms_d.A_0;
+            } else {
+            state_pert(i, j, k, RhoScalar_comp) = 0.0;
+            }
+        }
+
+        state_pert(i, j, k, RhoScalar_comp) *= parms_d.rho_0;
 
         if (use_moisture) {
             state_pert(i, j, k, RhoQ1_comp) = 0.0;
             state_pert(i, j, k, RhoQ2_comp) = 0.0;
         }
-      });
+    });
 
     // Set the x-velocity
     ParallelFor(xbx, [=, parms_d=parms] AMREX_GPU_DEVICE(int i, int j, int k) noexcept

Exec/DevTests/EB_Test/inputs

@@ -1,6 +1,5 @@
 # ------------------  INPUTS TO MAIN PROGRAM  -------------------
-max_step = 100
-max_step = 10
+max_step = 800
 
 amr.max_grid_size = 256 256 256
 
@@ -17,16 +16,16 @@ fabarray.mfiter_tile_size = 1024 1024 1024
 
 # PROBLEM SIZE & GEOMETRY
 geometry.prob_lo     = 0.   0.  0.
-geometry.prob_hi     = 16.  1.  16.
+geometry.prob_hi     = 16.  4.  8.
 
-amr.n_cell           = 64   4   64
+amr.n_cell           = 64   4   32
 
 geometry.is_periodic = 0 1 0
 
 xlo.type = Inflow
 xhi.type = Outflow
 
-xlo.velocity = 1. 0. 0.
+xlo.velocity = 10. 0. 0.
 xlo.density  = 1.0
 xlo.theta    = 1.0
 xlo.scalar   = 0.
@@ -36,7 +35,7 @@ zhi.type = SlipWall
 
 # TIME STEP CONTROL
 erf.substepping_type = None
-erf.fixed_dt         = 1.e-5
+erf.fixed_dt         = 1.e-3
 
 # DIAGNOSTICS & VERBOSITY
 erf.sum_interval   = 1       # timesteps between computing mass
@@ -72,4 +71,10 @@ erf.init_type = "uniform"
 # PROBLEM PARAMETERS
 prob.rho_0 = 1.0
 prob.T_0   = 1.0
-prob.u_0   = 1.0
+prob.u_0   = 10.0
+prob.A_0 = 1.0
+prob.prob_type = 11
+prob.xradius = 3.0
+prob.zradius = 1.5
+prob.xc_frac = 0.25
+prob.zc_frac = 0.15

Source/EB/ERF_redistribute.cpp → Source/EB/ERF_Redistribute.cpp

@@ -9,10 +9,11 @@ using namespace amrex;
 
 void
 ERF::redistribute_term ( int lev,
-                         MultiFab& result,
-                         MultiFab& result_tmp, // Saves doing a MF::copy. does this matter???
-                         MultiFab const& state,
-                         BCRec const* bc) // this is bc for the state (needed for SRD slopes)
+                    MultiFab& result,
+                    MultiFab& result_tmp, // Saves doing a MF::copy. does this matter???
+                    MultiFab const& state,
+                    BCRec const* bc, // this is bc for the state (needed for SRD slopes)
+                    Real const dt)
 {
     // ************************************************************************
     // Redistribute result_tmp and pass out result
@@ -26,16 +27,17 @@ ERF::redistribute_term ( int lev,
 #endif
     for (MFIter mfi(state,TilingIfNotGPU()); mfi.isValid(); ++mfi)
     {
-        redistribute_term(mfi, result, result_tmp, state, bc, lev);
+        redistribute_term(mfi, lev, result, result_tmp, state, bc, dt);
     }
 }
 
 void
 ERF::redistribute_term ( MFIter const& mfi, int lev,
-                         MultiFab& result,
-                         MultiFab& result_tmp,
-                         MultiFab const& state,
-                         BCRec const* bc) // this is bc for the state (needed for SRD slopes)
+                    MultiFab& result,
+                    MultiFab& result_tmp,
+                    MultiFab const& state,
+                    BCRec const* bc, // this is bc for the state (needed for SRD slopes)
+                    Real const dt)
 {
     AMREX_ASSERT(result.nComp() == state.nComp());
 
@@ -57,13 +59,13 @@ ERF::redistribute_term ( MFIter const& mfi, int lev,
         auto const& vfrac = ebfact.getVolFrac().const_array(mfi);
         auto const& ccc   = ebfact.getCentroid().const_array(mfi);
 
-        auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi);,
-        auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi);,
-        auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi););
+        auto const& apx = ebfact.getAreaFrac()[0]->const_array(mfi);
+        auto const& apy = ebfact.getAreaFrac()[1]->const_array(mfi);
+        auto const& apz = ebfact.getAreaFrac()[2]->const_array(mfi);
 
-        auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi);,
-        auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi);,
-        auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi););
+        auto const& fcx = ebfact.getFaceCent()[0]->const_array(mfi);
+        auto const& fcy = ebfact.getFaceCent()[1]->const_array(mfi);
+        auto const& fcz = ebfact.getFaceCent()[2]->const_array(mfi);
 
         Box gbx = bx; gbx.grow(3);
 
@@ -77,15 +79,15 @@ ERF::redistribute_term ( MFIter const& mfi, int lev,
         //     scratch(i,j,k) = 1.;
         // });
 
-        std::string redistribution_type = "StateRedistribution";
+        std::string redistribution_type = "StateRedist";
 
         // State redist acts on the state.
         Array4<Real const> state_arr = state.const_array(mfi);
         ApplyRedistribution( bx, ncomp, out, in, state_arr,
                              scratch, flag,
                              apx, apy, apz, vfrac,
                              fcx, fcy, fcz, ccc,
-                             bc, geom[lev], m_dt, edistribution_type);
+                             bc, geom[lev], dt, redistribution_type);
     }
     else
     {

Source/EB/Make.package

@@ -6,6 +6,5 @@ CEXE_sources += ERF_EBCylinder.cpp
 CEXE_sources += ERF_EBRegular.cpp
 CEXE_sources += ERF_Redistribute.cpp
 CEXE_sources += ERF_WriteEBSurface.cpp
-
 CEXE_headers += ERF_EBIF.H
 CEXE_headers += ERF_TerrainIF.H

Source/ERF.H

@@ -475,6 +475,18 @@ public:
     void make_eb_box ();
     void make_eb_cylinder ();
     void make_eb_regular ();
+    void redistribute_term ( int lev,
+                        amrex::MultiFab& result,
+                        amrex::MultiFab& result_tmp,
+                        amrex::MultiFab const& state,
+                        amrex::BCRec const* bc,
+                        amrex::Real const dt);
+    void redistribute_term ( amrex::MFIter const& mfi, int lev,
+                        amrex::MultiFab& result,
+                        amrex::MultiFab& result_tmp,
+                        amrex::MultiFab const& state,
+                        amrex::BCRec const* bc,
+                        amrex::Real const dt);
 #endif
 
     // more flexible version of AverageDown() that lets you average down across multiple levels

Source/TimeIntegration/ERF_TI_no_substep_fun.H

@@ -151,7 +151,13 @@
                 }
             } // mfi
 #ifdef ERF_USE_EB
-            // call to redistribute_term -- pass in ssum[IntVars::cons] which is MF
+            // // call to redistribute_term -- pass in ssum[IntVars::cons] which is MF
+            // MultiFab dUdt_tmp (ba, dm, n_data, 3, MFInfo(), Factory(level));
+            // dUdt_tmp.setVal(0.);
+            // dUdt_tmp.FillBoundary(fine_geom.periodicity());
+            // const BCRec* bc_ptr_h = domain_bcs_type.data(); // Should be host or device or both?
+            // redistribute_term ( level, F_slow[IntVars::cons], dUdt_tmp,
+            //                     S_sum[IntVars::cons], bc_ptr_h, slow_dt);
 #endif
         } // omp