fix initialization with wrfinput (#1166)

Co-authored-by: Aaron M. Lattanzi <[email protected]>

Source/ERF_make_new_level.cpp

@@ -171,9 +171,6 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
  z_t_rk[lev] = nullptr;
  }
 
- // Build the data structures for terrain-related quantities
- update_terrain_arrays(lev, time);
-
  // ********************************************************************************************
  // Define Theta_prim storage if using MOST BC
  // ********************************************************************************************
@@ -191,9 +188,24 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
 
  // ********************************************************************************************
  // Initialize the data itself
- // ********************************************************************************************
- if (restart_chkfile.empty()) {
- init_only(lev, start_time);
+ // If (init_type == "real") then we are initializing terrain and the initial data in
+ // the same call so we must call init_only before update_terrain_arrays
+ // If (init_type != "real") then we want to initialize the terrain before the initial data
+ // since we may need to use the grid information before constructing
+ // initial idealized data
+ // ********************************************************************************************
+ if (init_type == "real") {
+ if (restart_chkfile.empty()) {
+ init_only(lev, start_time);
+ }
+ // Build the data structures for terrain-related quantities
+ update_terrain_arrays(lev, time);
+ } else {
+ // Build the data structures for terrain-related quantities
+ update_terrain_arrays(lev, time);
+ if (restart_chkfile.empty()) {
+ init_only(lev, start_time);
+ }
  }
 }
 

Source/Initialization/ERF_init_from_wrfinput.cpp

@@ -79,7 +79,7 @@ init_msfs_from_wrfinput(int lev, FArrayBox& msfu_fab,
  const Vector<FArrayBox>& NC_MSFV_fab,
  const Vector<FArrayBox>& NC_MSFM_fab);
 void
-init_terrain_from_wrfinput(int lev, FArrayBox& z_phys,
+init_terrain_from_wrfinput(int lev, const Box& domain, FArrayBox& z_phys,
  const Vector<FArrayBox>& NC_PH_fab,
  const Vector<FArrayBox>& NC_PHB_fab);
 
@@ -180,13 +180,15 @@ ERF::init_from_wrfinput(int lev)
  NC_MSFU_fab, NC_MSFV_fab, NC_MSFM_fab);
  } // mf
 
+ const Box& domain = geom[lev].Domain();
+
  if (solverChoice.use_terrain)
  {
  std::unique_ptr<MultiFab>& z_phys = z_phys_nd[lev];
  for ( MFIter mfi(lev_new[Vars::cons], TilingIfNotGPU()); mfi.isValid(); ++mfi )
  {
  FArrayBox& z_phys_nd_fab = (*z_phys)[mfi];
- init_terrain_from_wrfinput(lev, z_phys_nd_fab, NC_PH_fab, NC_PHB_fab);
+ init_terrain_from_wrfinput(lev, domain, z_phys_nd_fab, NC_PH_fab, NC_PHB_fab);
  } // mf
 
  make_J (geom[lev],*z_phys_nd[lev],* detJ_cc[lev]);
@@ -230,8 +232,6 @@ ERF::init_from_wrfinput(int lev)
  // Without relaxation zones, we must augment this value by 1.
  if (wrfbdy_width == wrfbdy_set_width) wrfbdy_width += 1;
 
- const Box& domain = geom[lev].Domain();
-
  convert_wrfbdy_data(0,domain,bdy_data_xlo,
  NC_MUB_fab[0], NC_MSFU_fab[0], NC_MSFV_fab[0], NC_MSFM_fab[0],
  NC_PH_fab[0] , NC_PHB_fab[0],
@@ -408,7 +408,7 @@ init_base_state_from_wrfinput(int lev, const Box& valid_bx, const Real l_rdOcp,
  * @param NC_PHB_fab Vector of FArrayBox objects storing WRF terrain coordinate data (PHB)
  */
 void
-init_terrain_from_wrfinput(int lev, FArrayBox& z_phys,
+init_terrain_from_wrfinput(int lev, const Box& domain, FArrayBox& z_phys,
  const Vector<FArrayBox>& NC_PH_fab,
  const Vector<FArrayBox>& NC_PHB_fab)
 {
@@ -431,6 +431,9 @@ init_terrain_from_wrfinput(int lev, FArrayBox& z_phys,
  int klo = nodal_box.smallEnd()[2];
  int khi = nodal_box.bigEnd()[2]-1;
 
+ const auto domlo = amrex::lbound(domain);
+ const auto domhi = amrex::ubound(domain);
+
  //
  // We must be careful not to read out of bounds of the WPS data
  //
@@ -462,6 +465,20 @@ init_terrain_from_wrfinput(int lev, FArrayBox& z_phys,
  nc_ph_arr (ii,jj-1,k) + nc_ph_arr (ii-1,jj-1,k) +
  nc_phb_arr(ii,jj ,k) + nc_phb_arr(ii-1,jj ,k) +
  nc_phb_arr(ii,jj-1,k) + nc_phb_arr(ii-1,jj-1,k) ) / CONST_GRAV;
+
+ //
+ // Fill values outside the domain on the fly (we will need these to make detJ in ghost cells)
+ //
+ if (i == domlo.x) z_arr(i-1,j,k) = z_arr(i,j,k);
+ if (i == domhi.x) z_arr(i+1,j,k) = z_arr(i,j,k);
+ if (j == domlo.y) z_arr(i,j-1,k) = z_arr(i,j,k);
+ if (j == domhi.y) z_arr(i,j+1,k) = z_arr(i,j,k);
+
+ if (i == domlo.x && j == domlo.y) z_arr(i-1,j-1,k) = z_arr(i,j,k);
+ if (i == domlo.x && j == domhi.y) z_arr(i-1,j+1,k) = z_arr(i,j,k);
+ if (i == domhi.x && j == domlo.y) z_arr(i+1,j-1,k) = z_arr(i,j,k);
+ if (i == domhi.x && j == domhi.y) z_arr(i+1,j+1,k) = z_arr(i,j,k);
+
  } // k
  });
  } // idx