Merge SRD changes from amrex PR#3491

Fix a bug in StateRedistribution which occurred when the EB was
near a fine-fine boundary. In this case we weren't calculating
certain arrays in a wide enough band of ghost cells. This has been
fixed by enlarging two of the arrays and calculating them on the
correct number of ghost cells.

src/redistribution/hydro_redistribution.cpp

@@ -201,19 +201,18 @@ void Redistribution::Apply ( Box const& bx, int ncomp,
  } else if (redistribution_type == "StateRedist") {
 
  Box const& bxg1 = grow(bx,1);
- Box const& bxg2 = grow(bx,2);
  Box const& bxg3 = grow(bx,3);
  Box const& bxg4 = grow(bx,4);
 
  // We use the first component of this for the number of neighbors, and later
  // components identify the neighbors (utilizing the CellMap)
  IArrayBox itracker(bxg4,itracker_comp,The_Async_Arena());
 
- FArrayBox nrs_fab(bxg3,1,The_Async_Arena());
+ FArrayBox nrs_fab(bxg4,1,The_Async_Arena());
  FArrayBox alpha_fab(bxg3,2,The_Async_Arena());
 
  // Total volume of all cells in my nbhd
- FArrayBox nbhd_vol_fab(bxg2,1,The_Async_Arena());
+ FArrayBox nbhd_vol_fab(bxg3,1,The_Async_Arena());
 
  // Centroid of my nbhd
  FArrayBox cent_hat_fab(bxg3,AMREX_SPACEDIM,The_Async_Arena());
@@ -589,7 +588,6 @@ Redistribution::ApplyToInitialData ( Box const& bx, int ncomp,
  amrex::Error(msg);
  }
 
- Box const& bxg2 = grow(bx,2);
  Box const& bxg3 = grow(bx,3);
  Box const& bxg4 = grow(bx,4);
 
@@ -603,11 +601,11 @@ Redistribution::ApplyToInitialData ( Box const& bx, int ncomp,
  // use the first component of this for the number of neighbors
  IArrayBox itracker(bxg4,itracker_comp,The_Async_Arena());
 #endif
- FArrayBox nrs_fab(bxg3,1,The_Async_Arena());
+ FArrayBox nrs_fab(bxg4,1,The_Async_Arena());
  FArrayBox alpha_fab(bxg3,2,The_Async_Arena());
 
  // Total volume of all cells in my nbhd
- FArrayBox nbhd_vol_fab(bxg2,1,The_Async_Arena());
+ FArrayBox nbhd_vol_fab(bxg3,1,The_Async_Arena());
 
  // Centroid of my nbhd
  FArrayBox cent_hat_fab(bxg3,AMREX_SPACEDIM,The_Async_Arena());

src/redistribution/hydro_state_utils.cpp

@@ -71,11 +71,16 @@ Redistribution::MakeStateRedistUtils ( Box const& bx,
  if (is_periodic_z) domain_per_grown.grow(2,2);
 #endif
 
- amrex::ParallelFor(bxg3,
+ amrex::ParallelFor(bxg4,
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  // Everyone is in their own neighborhood at least
  nrs(i,j,k) = 1.;
+ });
+
+ amrex::ParallelFor(bxg3,
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ {
  alpha(i,j,k,0) = 1.;
  alpha(i,j,k,1) = 0.;
  });
@@ -91,14 +96,14 @@ Redistribution::MakeStateRedistUtils ( Box const& bx,
  int s = j+jmap[itracker(i,j,k,i_nbor)];
  int t = k+kmap[itracker(i,j,k,i_nbor)];
  if ( domain_per_grown.contains(IntVect(AMREX_D_DECL(r,s,t))) &&
- bxg3.contains(IntVect(AMREX_D_DECL(r,s,t))) )
+ bxg4.contains(IntVect(AMREX_D_DECL(r,s,t))) )
  {
  amrex::Gpu::Atomic::Add(&nrs(r,s,t),1.0_rt);
  }
  }
  });
 
- amrex::ParallelFor(bxg2,
+ amrex::ParallelFor(bxg3,
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  if (vfrac_old(i,j,k) > 0.0 || vfrac_new(i,j,k) > 0.0)
@@ -136,7 +141,7 @@ Redistribution::MakeStateRedistUtils ( Box const& bx,
  });
 
  // Define how much each cell keeps
- amrex::ParallelFor(bxg2,
+ amrex::ParallelFor(bxg3,
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  if (vfrac_new(i,j,k) > 0.0 || vfrac_old(i,j,k) > 0.0)
@@ -148,13 +153,15 @@ Redistribution::MakeStateRedistUtils ( Box const& bx,
  int r = i+imap[itracker(i,j,k,i_nbor)];
  int s = j+jmap[itracker(i,j,k,i_nbor)];
  int t = k+kmap[itracker(i,j,k,i_nbor)];
- amrex::Gpu::Atomic::Add(&alpha(r,s,t,0),-(alpha(i,j,k,1)/nrs(r,s,t)));
+ if ( bxg3.contains(IntVect(AMREX_D_DECL(r,s,t))) ) {
+ amrex::Gpu::Atomic::Add(&alpha(r,s,t,0),-(alpha(i,j,k,1)/nrs(r,s,t)));
+ }
  }
  }
  });
 
  // Redefine nbhd_vol
- amrex::ParallelFor(bxg2,
+ amrex::ParallelFor(bxg3,
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  if (vfrac_new(i,j,k) > 0.0 || vfrac_old(i,j,k) > 0.0)
@@ -188,7 +195,7 @@ Redistribution::MakeStateRedistUtils ( Box const& bx,
  });
 
  // Define xhat,yhat,zhat (from Berger and Guliani)
- amrex::ParallelFor(bxg3,
+ amrex::ParallelFor(bxg2,
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  if (vfrac_new(i,j,k) > 0.0)