Streamline FPr BA (#1199)

* Fixes for use with amrex patch for face linear interp.

* not sure why this is needed with substepping.

* This requires amrex patch 3483. Also, do relaxation on all RK stages; from tests with 2d cylinder.

---------

Co-authored-by: Ann Almgren <[email protected]>

Exec/RegTests/DensityCurrent/inputs_periodic_twolevel

@@ -0,0 +1,68 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 999999
+stop_time = 900.0
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_lo     = -25600.    0.      0.
+geometry.prob_hi     =  25600.  400.  6400.
+
+amr.n_cell           =   512   12     64   # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000
+
+# periodic in x to match WRF setup
+# - as an alternative, could use symmetry at x=0 and outflow at x=25600
+geometry.is_periodic = 1 1 0
+zlo.type = "SlipWall"
+zhi.type = "SlipWall"
+
+# TIME STEP CONTROL
+erf.fixed_dt       = 0.1       # fixed time step [s] -- Straka et al 1993
+erf.fixed_fast_dt  = 0.025     # fixed time step [s] -- Straka et al 1993
+
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1       # timesteps between computing mass
+erf.v              = 1       # verbosity in ERF.cpp
+amr.v              = 1       # verbosity in Amr.cpp
+
+# MULTILEVEL FLAGS
+amr.max_level = 1
+amr.ref_ratio_vect = 2 2 1
+erf.refinement_indicators = box1
+erf.box1.max_level = 1
+erf.box1.in_box_lo =   -6400.  0.
+erf.box1.in_box_hi =    6400. 400.
+erf.coupling_type  = "OneWay"
+erf.cf_width = 5
+erf.cf_set_width = 1
+
+# CHECKPOINT FILES
+erf.check_file      = chk        # root name of checkpoint file
+erf.check_int       = -57600      # number of timesteps between checkpoints
+
+# PLOTFILES
+erf.plot_file_1     = plt        # prefix of plotfile name
+erf.plot_int_1      = 1       # number of timesteps between plotfiles
+erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta pres_hse dens_hse pert_pres pert_dens
+
+# SOLVER CHOICE
+erf.use_gravity = true
+erf.use_coriolis = false
+erf.use_rayleigh_damping = false
+
+#
+# diffusion coefficient from Straka, K = 75 m^2/s
+#
+erf.les_type         = "None"
+erf.molec_diff_type = "ConstantAlpha"
+erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities
+erf.dynamicViscosity = 75.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s
+erf.alpha_T = 75.0 # [m^2/s]
+
+erf.c_p = 1004.5
+
+# PROBLEM PARAMETERS (optional)
+prob.T_0 = 300.0
+prob.U_0 = 10.0

Exec/RegTests/Terrain3d_Hemisphere/inputs.amr

@@ -0,0 +1,94 @@
+# ------------------  INPUTS TO MAIN PROGRAM  -------------------
+max_step = 4000
+
+amrex.fpe_trap_invalid = 1
+
+fabarray.mfiter_tile_size = 1024 1024 1024
+
+# PROBLEM SIZE & GEOMETRY
+geometry.prob_lo     = 0.   0.  0. 
+geometry.prob_hi     = 10  10  10
+
+amr.n_cell           = 100  100   100   # dx=dy=dz=100 m, Straka et al 1993 / Xue et al 2000
+
+geometry.is_periodic = 0 1 0
+
+xlo.type = "Outflow"
+xhi.type = "Outflow"
+xlo.velocity = 1. 0. 0.
+xlo.density  = 1.16
+xlo.theta    = 300.
+xlo.scalar   = 0.
+
+zlo.type = "SlipWall"
+zhi.type = "Outflow"
+
+erf.sponge_strength = 100.0
+erf.use_xlo_sponge_damping = true
+erf.xlo_sponge_end = 2.0
+erf.use_xhi_sponge_damping = true
+erf.xhi_sponge_start = 8.0
+erf.use_ylo_sponge_damping = true
+erf.ylo_sponge_end = 2.0
+erf.use_yhi_sponge_damping = true
+erf.yhi_sponge_start = 8.0
+erf.use_zhi_sponge_damping = true
+erf.zhi_sponge_start = 8.0
+
+erf.sponge_density = 1.2
+erf.sponge_x_velocity = 10.0
+erf.sponge_y_velocity = 0.0
+erf.sponge_z_velocity = 0.0
+
+
+# TIME STEP CONTROL
+erf.no_substepping     = 1
+#erf.fixed_dt           = 1E-5
+
+erf.cfl = 0.3
+
+# DIAGNOSTICS & VERBOSITY
+erf.sum_interval   = 1        # timesteps between computing mass
+erf.v              = 1        # verbosity in ERF.cpp
+amr.v              = 1        # verbosity in Amr.cpp
+
+# REFINEMENT / REGRIDDING
+amr.max_level = 1
+amr.ref_ratio_vect = 2 2 1
+erf.refinement_indicators = box1
+erf.box1.max_level = 1
+erf.box1.in_box_lo =   3.   3.
+erf.box1.in_box_hi =   7.   7.
+erf.coupling_type  = "OneWay"
+erf.cf_width = 0 #5
+erf.cf_set_width = 0 #1
+
+# CHECKPOINT FILES
+erf.check_file      = chk     # root name of checkpoint file
+erf.check_int       = -57600  # number of timesteps between checkpoints
+
+# PLOTFILES
+erf.plot_file_1     = plt     # prefix of plotfile name
+erf.plot_int_1      = 100       # number of timesteps between plotfiles
+erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta pres_hse dens_hse pert_pres pert_dens z_phys detJ dpdx dpdy pres_hse_x pres_hse_y
+
+# SOLVER CHOICE
+erf.use_gravity = true
+erf.use_coriolis = false
+erf.use_rayleigh_damping = false
+erf.les_type = "None"
+
+# TERRRAIN GRID TYPE
+erf.use_terrain = true
+erf.terrain_smoothing = 0
+
+# Diffusion coefficient from Straka, K = 75 m^2/s
+erf.molec_diff_type = "ConstantAlpha"
+erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities
+erf.dynamicViscosity = 0.1 # [kg/(m-s)] ==> nu = 75.0 m^2/s
+erf.alpha_T = 0.0 # [m^2/s]
+
+# PROBLEM PARAMETERS (optional)
+prob.T_0   = 300.0
+prob.U_0   = 0.0
+prob.rho_0 = 1.16

Source/BoundaryConditions/ERF_FillPatcher.H

@@ -67,7 +67,7 @@ void
 ERFFillPatcher::fill (amrex::MultiFab& mf, amrex::Real time,
                       BC& cbc, amrex::Vector<amrex::BCRec> const& bcs, bool fill_subset)
 {
-    constexpr amrex::Real eps = std::numeric_limits<amrex::Real>::epsilon();
+    constexpr amrex::Real eps = std::numeric_limits<float>::epsilon();
     AMREX_ALWAYS_ASSERT((time >= m_crse_times[0]-eps) && (time <= m_crse_times[1]+eps));
 
     // Time interpolation factors

Source/BoundaryConditions/ERF_FillPatcher.cpp

@@ -116,49 +116,31 @@ void ERFFillPatcher::Define (BoxArray const& fba, DistributionMapping  fdm,
         }
     }
 
-    // Grown fine box list (expand interior halo to touch crse face)
-    BoxList gbl;
-    gbl.set(m_ixt);
-    for (auto bl_iter = bl.begin(); bl_iter < bl.end(); bl_iter++) {
-      const int* lo = (*bl_iter).loVect();
-      const int* hi = (*bl_iter).hiVect();
-      IntVect glo(0); IntVect ghi(0);
-      for (int idim(0); idim < AMREX_SPACEDIM; ++idim) {
-        int rlo = lo[idim]%m_ratio[idim];
-        int rhi = hi[idim]%m_ratio[idim];
-        glo[idim] = lo[idim]; ghi[idim] = hi[idim];
-        if (rlo > 0) glo[idim] = lo[idim] - rlo;
-        if (rhi > 0) ghi[idim] = hi[idim] - rhi + m_ratio[idim];
-      }
-      Box gbx(glo,ghi,m_ixt);
-      gbl.push_back(gbx);
-    }
-
     // Coarse box list (coincides with grown fine box list)
     BoxList cbl;
     cbl.set(m_ixt);
-    cbl.reserve(gbl.size());
-    for (auto const& b : gbl) {
+    cbl.reserve(bl.size());
+    for (auto const& b : bl) {
         cbl.push_back(interp->CoarseBox(b, m_ratio));
     }
 
     // Box arrays for fine and coarse
-    BoxArray gcf_fba(std::move(gbl));
-    BoxArray gcf_cba(std::move(cbl));
+    BoxArray cf_fba(std::move(bl));
+    BoxArray cf_cba(std::move(cbl));
     BoxArray cf_fba_s(std::move(bl_subset));
-    DistributionMapping gcf_dm(gcf_fba);
+    DistributionMapping cf_dm(cf_fba);
     DistributionMapping cf_dm_s(cf_fba_s);
 
 
     // Fine patch to hold the time-interpolated state
-    m_cf_fine_data = new MultiFab (gcf_fba, gcf_dm, m_ncomp, 0);
+    m_cf_fine_data = new MultiFab (cf_fba, cf_dm, m_ncomp, 0);
 
     // Fine subset patch to hold the time-interpolated state
-    m_cf_fine_subset_data = new MultiFab (cf_fba_s, gcf_dm, m_ncomp, 0);
+    m_cf_fine_subset_data = new MultiFab (cf_fba_s, cf_dm, m_ncomp, 0);
 
     // Two coarse patches to hold the data to be interpolated
-    m_cf_crse_data[0] = new MultiFab (gcf_cba, gcf_dm, m_ncomp, 0);
-    m_cf_crse_data[1] = new MultiFab (gcf_cba, gcf_dm, m_ncomp, 0);
+    m_cf_crse_data[0] = new MultiFab (cf_cba, cf_dm, m_ncomp, 0);
+    m_cf_crse_data[1] = new MultiFab (cf_cba, cf_dm, m_ncomp, 0);
 }
 
 

Source/TimeIntegration/TI_slow_rhs_fun.H

@@ -194,7 +194,7 @@
 
 #ifdef ERF_USE_NETCDF
         // Populate RHS for relaxation zones
-        if (init_type=="real" && level==0 && nrk==0) {
+        if (init_type=="real" && level==0) {
             wrfbdy_compute_interior_ghost_RHS(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt,
                                               wrfbdy_width-1, wrfbdy_set_width, fine_geom,
                                               S_rhs, S_data, bdy_data_xlo, bdy_data_xhi,
@@ -203,7 +203,7 @@
 #endif
 
         // Compute RHS for fine interior ghost
-        if (level>0 && coupling_type=="OneWay" && cf_width>0 && nrk==0) {
+        if (level>0 && coupling_type=="OneWay" && cf_width>0) {
             fine_compute_interior_ghost_RHS(new_stage_time, slow_dt,
                                             cf_width-1, cf_set_width, fine_geom,
                                             &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1],
@@ -342,7 +342,7 @@
 
 #ifdef ERF_USE_NETCDF
         // Populate RHS for relaxation zones
-        if (init_type=="real" && level==0 && nrk==0) {
+        if (init_type=="real" && level==0) {
             wrfbdy_compute_interior_ghost_RHS(bdy_time_interval, start_bdy_time, new_stage_time, slow_dt,
                                               wrfbdy_width-1, wrfbdy_set_width, fine_geom,
                                               S_rhs, S_data, bdy_data_xlo, bdy_data_xhi,
@@ -351,7 +351,7 @@
 #endif
 
         // Compute RHS for fine interior ghost
-        if (level>0 && coupling_type=="OneWay" && cf_width>0 && nrk==0) {
+        if (level>0 && coupling_type=="OneWay" && cf_width>0) {
             fine_compute_interior_ghost_RHS(new_stage_time, slow_dt,
                                             cf_width-1, cf_set_width,
                                             &FPr_c[level-1], &FPr_u[level-1], &FPr_v[level-1], &FPr_w[level-1],

Source/Utils/InteriorGhostCells.cpp

@@ -577,6 +577,7 @@ fine_compute_interior_ghost_RHS(const Real& time,
                     Box domain  = boxes_at_level[g_ind];
                     domain.convert(fmf.boxArray().ixType());
 
+                    // NOTE:: width+1 to get halo interior cell from fill
                     Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
                     compute_interior_ghost_bxs_xy(tbx, domain, width+1, 0,
                                                   bx_xlo, bx_xhi,
@@ -625,6 +626,7 @@ fine_compute_interior_ghost_RHS(const Real& time,
                     Box domain  = boxes_at_level[g_ind];
                     domain.convert(fmf.boxArray().ixType());
 
+                    // NOTE:: width+1 to get halo interior cell from fill
                     Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
                     compute_interior_ghost_bxs_xy(tbx, domain, width+1, 0,
                                                   bx_xlo, bx_xhi,
@@ -672,6 +674,7 @@ fine_compute_interior_ghost_RHS(const Real& time,
                     Box domain = boxes_at_level[g_ind];
                     domain.convert(fmf.boxArray().ixType());
 
+                    // NOTE:: width+1 to get halo interior cell from fill
                     Box bx_xlo, bx_xhi, bx_ylo, bx_yhi;
                     compute_interior_ghost_bxs_xy(tbx, domain, width+1, 0,
                                                   bx_xlo, bx_xhi,