1) initial pressure projection now plays correctly with gp0, 2) the p… (

#110)

…lotfile variables gpx, gpy, gpz now contain (gp0+gp) instead of just
the gp variables

src/boundary_conditions/incflo_set_bcs.cpp

@@ -89,7 +89,7 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> const& bcs,
  if (m_bc_type[olo] == BC::mixed) {
  prob_set_BC_MF(olo, blo, bc_arr, lev, inflow, outflow, field);
  } else {
- amrex::ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(blo, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  for (int n = 0; n < ncomp; n++) {
  bc_arr(i,j,k,n) = bc_ptr[n].lo(dir);
@@ -99,7 +99,7 @@ incflo::make_BC_MF(int lev, amrex::Gpu::DeviceVector<amrex::BCRec> const& bcs,
  if (m_bc_type[ohi] == BC::mixed) {
  prob_set_BC_MF(ohi, bhi, bc_arr, lev, inflow, outflow, field);
  } else {
- amrex::ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bhi, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  for (int n = 0; n < ncomp; n++) {
  bc_arr(i,j,k,n) = bc_ptr[n].hi(dir);
@@ -189,7 +189,7 @@ incflo::make_robinBC_MFs(int lev, MultiFab* state)
 
 #ifdef AMREX_USE_EB
 void
-incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int nghost)
+incflo::set_eb_velocity (int lev, Real /*time*/, MultiFab& eb_vel, int nghost)
 {
  Geometry const& gm = Geom(lev);
  eb_vel.setVal(0.);
@@ -224,7 +224,7 @@ incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int ng
 
  bool has_comps = false;
  Real eb_vel_mag(0.0);
- GpuArray<amrex::Real,3> eb_vel_comps{0.};
+ GpuArray<Real,3> eb_vel_comps{0.};
  // Flow is specified as a velocity magnitude
  if ( m_eb_flow.is_mag ) {
  eb_vel_mag = m_eb_flow.vel_mag;
@@ -281,7 +281,7 @@ incflo::set_eb_velocity (int lev, amrex::Real /*time*/, MultiFab& eb_vel, int ng
 }
 
 void
-incflo::set_eb_density (int lev, amrex::Real /*time*/, MultiFab& eb_density, int nghost)
+incflo::set_eb_density (int lev, Real /*time*/, MultiFab& eb_density, int nghost)
 {
  Geometry const& gm = Geom(lev);
  eb_density.setVal(0.);
@@ -350,7 +350,7 @@ incflo::set_eb_density (int lev, amrex::Real /*time*/, MultiFab& eb_density, int
 }
 
 void
-incflo::set_eb_tracer (int lev, amrex::Real /*time*/, MultiFab& eb_tracer, int nghost)
+incflo::set_eb_tracer (int lev, Real /*time*/, MultiFab& eb_tracer, int nghost)
 {
  Geometry const& gm = Geom(lev);
  eb_tracer.setVal(0.);

src/convection/incflo_compute_MAC_projected_velocities.cpp

@@ -36,15 +36,15 @@ incflo::compute_MAC_projected_velocities (
  Array4<Real const> const& rho = density[lev]->array(mfi);
  Array4<Real const> const& divtau = ld.divtau_o.const_array(mfi);
  if (m_advect_momentum) {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  AMREX_D_TERM(vel_f(i,j,k,0) += divtau(i,j,k,0)/rho(i,j,k);,
  vel_f(i,j,k,1) += divtau(i,j,k,1)/rho(i,j,k);,
  vel_f(i,j,k,2) += divtau(i,j,k,2)/rho(i,j,k););
  });
  }
  else {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  AMREX_D_TERM(vel_f(i,j,k,0) += divtau(i,j,k,0);,
  vel_f(i,j,k,1) += divtau(i,j,k,1);,
@@ -80,7 +80,7 @@ incflo::compute_MAC_projected_velocities (
  EB_interp_CellCentroid_to_FaceCentroid (*density[lev], GetArrOfPtrs(inv_rho[lev]),
  0, 0, 1, geom[lev], get_density_bcrec());
 #else
- amrex::average_cellcenter_to_face(GetArrOfPtrs(inv_rho[lev]), *density[lev], geom[lev]);
+ average_cellcenter_to_face(GetArrOfPtrs(inv_rho[lev]), *density[lev], geom[lev]);
 #endif
 
  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {

src/convection/incflo_compute_advection_term.cpp

@@ -154,15 +154,15 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  Array4<Real const> const& rho = density[lev]->array(mfi);
  Array4<Real const> const& divtau = ld.divtau_o.const_array(mfi);
  if (m_advect_momentum) {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  AMREX_D_TERM(vel_f(i,j,k,0) += divtau(i,j,k,0)/rho(i,j,k);,
  vel_f(i,j,k,1) += divtau(i,j,k,1)/rho(i,j,k);,
  vel_f(i,j,k,2) += divtau(i,j,k,2)/rho(i,j,k););
  });
  }
  else {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  AMREX_D_TERM(vel_f(i,j,k,0) += divtau(i,j,k,0);,
  vel_f(i,j,k,1) += divtau(i,j,k,1);,
@@ -345,7 +345,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  Array4<Real const> rho = density[lev]->const_array(mfi);
  Array4<Real > rho_vel = rhovel[lev].array(mfi);
 
- amrex::ParallelFor(bxg, AMREX_SPACEDIM,
+ ParallelFor(bxg, AMREX_SPACEDIM,
  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
  {
  rho_vel(i,j,k,n) = rho(i,j,k) * U(i,j,k,n);
@@ -358,7 +358,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  Array4<Real const> vf = vel_forces[lev]->const_array(mfi);
  Array4<Real > rho_vel_f = rhovel_f.array();
 
- amrex::ParallelFor(fbx, AMREX_SPACEDIM,
+ ParallelFor(fbx, AMREX_SPACEDIM,
  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
  {
  rho_vel_f(i,j,k,n) = rho(i,j,k) * vf(i,j,k,n);
@@ -456,7 +456,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  if ( any_conserv_trac ) {
  trac_tmp = rhotrac[lev].array(mfi);
 
- amrex::ParallelFor(bxg, m_ntrac,
+ ParallelFor(bxg, m_ntrac,
  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
  {
  if ( iconserv[n] ){
@@ -560,7 +560,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  Array4<Real const> rho = get_density_eb()[lev]->const_array(mfi);
  Array4<Real > rho_vel = rhovel_fab.array();
 
- amrex::ParallelFor(bxg, AMREX_SPACEDIM,
+ ParallelFor(bxg, AMREX_SPACEDIM,
  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
  {
  rho_vel(i,j,k,n) = rho(i,j,k) * U(i,j,k,n);
@@ -741,7 +741,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  *density[lev], bc_den, lev);
  } else {
  auto const& drdt = conv_r[lev]->array(mfi);
- amrex::ParallelFor(bx,
+ ParallelFor(bx,
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  drdt(i,j,k) = 0.;

src/derive/incflo_derive.cpp

@@ -73,23 +73,23 @@ void incflo::compute_strainrate_at_level (int /*lev*/,
  auto typ = flag_fab.getType(bx);
  if (typ == FabType::covered)
  {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  sr_arr(i,j,k) = Real(0.0);
  });
  }
  else if (typ == FabType::singlevalued)
  {
  auto const& flag_arr = flag_fab.const_array();
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  sr_arr(i,j,k) = incflo_strainrate_eb(i,j,k,AMREX_D_DECL(idx,idy,idz),vel_arr,flag_arr(i,j,k));
  });
  }
  else
 #endif
  {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  sr_arr(i,j,k) = incflo_strainrate(i,j,k,AMREX_D_DECL(idx,idy,idz),vel_arr);
  });
@@ -171,7 +171,7 @@ void incflo::ComputeMagVel (int /*lev*/,
  auto typ = flags.getType(bx);
  if (typ == FabType::covered)
  {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  magvel_fab(i,j,k) = Real(0.0);
  });
@@ -240,15 +240,15 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
  auto typ = flags.getType(bx);
  if (typ == FabType::covered)
  {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  vort_fab(i,j,k) = Real(0.0);
  });
  }
  else if (typ == FabType::singlevalued)
  {
  const auto& flag_fab = flags.const_array();
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  constexpr Real c0 = Real(-1.5);
  constexpr Real c1 = Real( 2.0);
@@ -306,7 +306,7 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
  else
 #endif
  {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  Real vx = Real(0.5) * (ccvel_fab(i+1,j,k,1) - ccvel_fab(i-1,j,k,1)) * idx;
  Real uy = Real(0.5) * (ccvel_fab(i,j+1,k,0) - ccvel_fab(i,j-1,k,0)) * idy;
@@ -343,15 +343,15 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
  auto typ = flags.getType(bx);
  if (typ == FabType::covered)
  {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  vort_fab(i,j,k) = Real(0.0);
  });
  }
  else if (typ == FabType::singlevalued)
  {
  const auto& flag_fab = flags.const_array();
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  constexpr Real c0 = Real(-1.5);
  constexpr Real c1 = Real( 2.0);
@@ -447,7 +447,7 @@ void incflo::ComputeVorticity (int lev, Real /*time*/, MultiFab& vort, MultiFab
  else
 #endif
  {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  Real vx = Real(0.5) * (ccvel_fab(i+1,j,k,1) - ccvel_fab(i-1,j,k,1)) * idx;
  Real wx = Real(0.5) * (ccvel_fab(i+1,j,k,2) - ccvel_fab(i-1,j,k,2)) * idx;

src/derive/incflo_error.cpp

@@ -20,7 +20,7 @@ void incflo::DiffFromExact (int /*lev*/, Geometry& lev_geom, Real time, Real dt,
  // When we enter this routine, this holds the computed solution
  Array4<Real> const& err = error.array(mfi);
 
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  constexpr Real twopi = Real(2.0)*Real(3.1415926535897932);
  constexpr Real fourpi = Real(4.0)*Real(3.1415926535897932);
@@ -75,7 +75,7 @@ void incflo::DiffFromExact (int /*lev*/, Geometry& lev_geom, Real time, Real dt,
 
  Real omega = pi * pi * visc_coef;
 
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
 
  Real x = Real(i+0.5)*dx[0];
@@ -117,7 +117,7 @@ void incflo::DiffFromExact (int /*lev*/, Geometry& lev_geom, Real time, Real dt,
  // When we enter this routine, this holds the computed solution
  Array4<Real> const& err = error.array(mfi);
 
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  Real x = Real(i+0.5)*dx[0];
  Real y = Real(j+0.5)*dx[1];

src/diffusion/DiffusionScalarOp.cpp

@@ -260,8 +260,7 @@ DiffusionScalarOp::diffuse_scalar (Vector<MultiFab*> const& tracer,
  Array4<Real> const& rhs_a = rhs[lev].array(mfi);
  Array4<Real const> const& tra_a = tracer[lev]->const_array(mfi,comp);
  Array4<Real const> const& rho_a = density[lev]->const_array(mfi);
- amrex::ParallelFor(bx,
- [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  rhs_a(i,j,k) = rho_a(i,j,k) * tra_a(i,j,k);
  });
@@ -412,8 +411,7 @@ DiffusionScalarOp::diffuse_vel_components (Vector<MultiFab*> const& vel,
  Array4<Real> const& rhs_a = rhs[lev].array(mfi);
  Array4<Real const> const& vel_a = vel[lev]->const_array(mfi,comp);
  Array4<Real const> const& rho_a = density[lev]->const_array(mfi);
- amrex::ParallelFor(bx,
- [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  rhs_a(i,j,k) = rho_a(i,j,k) * vel_a(i,j,k);
  });
@@ -714,8 +712,7 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
  Box const& bx = mfi.tilebox();
  Array4<Real> const& divtau_arr = a_divtau[lev]->array(mfi);
  Array4<Real const> const& rho_arr = a_density[lev]->const_array(mfi);
- amrex::ParallelFor(bx,
- [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  Real rhoinv = Real(1.0)/rho_arr(i,j,k);
  AMREX_D_TERM(divtau_arr(i,j,k,0) *= rhoinv;,

src/diffusion/DiffusionTensorOp.cpp

@@ -164,7 +164,7 @@ DiffusionTensorOp::diffuse_velocity (Vector<MultiFab*> const& velocity,
  Array4<Real> const& rhs_a = rhs[lev].array(mfi);
  Array4<Real const> const& vel_a = velocity[lev]->const_array(mfi);
  Array4<Real const> const& rho_a = density[lev]->const_array(mfi);
- amrex::ParallelFor(bx,AMREX_SPACEDIM,
+ ParallelFor(bx,AMREX_SPACEDIM,
  [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
  {
  rhs_a(i,j,k,n) = rho_a(i,j,k) * vel_a(i,j,k,n);
@@ -301,8 +301,7 @@ void DiffusionTensorOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
  Box const& bx = mfi.tilebox();
  Array4<Real> const& divtau_arr = a_divtau[lev]->array(mfi);
  Array4<Real const> const& rho_arr = a_density[lev]->const_array(mfi);
- amrex::ParallelFor(bx,
- [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+ ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  Real rhoinv = Real(1.0)/rho_arr(i,j,k);
  AMREX_D_TERM(divtau_arr(i,j,k,0) *= rhoinv;,

src/diffusion/incflo_diffusion.cpp

@@ -325,14 +325,14 @@ incflo::fixup_eta_on_domain_faces (int lev, Array<MultiFab,AMREX_SPACEDIM>& fc,
  if (!gm.isPeriodic(idim)) {
  Array4<Real> const& fca = fc[idim].array(mfi);
  if (bx.smallEnd(idim) == domain.smallEnd(idim)) {
- amrex::ParallelFor(amrex::bdryLo(bx, idim),
+ ParallelFor(amrex::bdryLo(bx, idim),
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  fca(i,j,k) = cca(i,j,k);
  });
  }
  if (bx.bigEnd(idim) == domain.bigEnd(idim)) {
- amrex::ParallelFor(amrex::bdryHi(bx, idim),
+ ParallelFor(amrex::bdryHi(bx, idim),
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  fca(i,j,k) = cca(i-1,j,k);
@@ -344,14 +344,14 @@ incflo::fixup_eta_on_domain_faces (int lev, Array<MultiFab,AMREX_SPACEDIM>& fc,
  if (!gm.isPeriodic(idim)) {
  Array4<Real> const& fca = fc[idim].array(mfi);
  if (bx.smallEnd(idim) == domain.smallEnd(idim)) {
- amrex::ParallelFor(amrex::bdryLo(bx, idim),
+ ParallelFor(amrex::bdryLo(bx, idim),
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  fca(i,j,k) = cca(i,j,k);
  });
  }
  if (bx.bigEnd(idim) == domain.bigEnd(idim)) {
- amrex::ParallelFor(amrex::bdryHi(bx, idim),
+ ParallelFor(amrex::bdryHi(bx, idim),
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  fca(i,j,k) = cca(i,j-1,k);
@@ -364,14 +364,14 @@ incflo::fixup_eta_on_domain_faces (int lev, Array<MultiFab,AMREX_SPACEDIM>& fc,
  if (!gm.isPeriodic(idim)) {
  Array4<Real> const& fca = fc[idim].array(mfi);
  if (bx.smallEnd(idim) == domain.smallEnd(idim)) {
- amrex::ParallelFor(amrex::bdryLo(bx, idim),
+ ParallelFor(amrex::bdryLo(bx, idim),
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  fca(i,j,k) = cca(i,j,k);
  });
  }
  if (bx.bigEnd(idim) == domain.bigEnd(idim)) {
- amrex::ParallelFor(amrex::bdryHi(bx, idim),
+ ParallelFor(amrex::bdryHi(bx, idim),
  [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
  {
  fca(i,j,k) = cca(i,j,k-1);

src/incflo.H

@@ -316,6 +316,7 @@ public:
 
  static bool SteadyStateReached ();
 
+ void InitialPressureProjection ();
 
 private:
 
@@ -920,7 +921,6 @@ private:
  void ReadIOParameters ();
  void ResizeArrays (); // Resize arrays to fit (up to) max_level + 1 AMR levels
  void InitialProjection ();
- void InitialPressureProjection ();
  void InitialIterations ();
 #ifdef AMREX_USE_EB
  void InitialRedistribution ();