Merge branch 'development' into add_derived

src/CMakeLists.txt

@@ -10,6 +10,7 @@ target_sources(incflo
  incflo_compute_forces.cpp
  incflo_correct_small_cells.cpp
  incflo_redistribute.cpp
+ incflo_explicit_update.cpp
  incflo_tagging.cpp
  incflo_regrid.cpp
  main.cpp

src/Make.package

@@ -4,6 +4,7 @@ CEXE_sources += incflo_apply_predictor.cpp
 CEXE_sources += incflo_apply_corrector.cpp
 CEXE_sources += incflo_compute_dt.cpp
 CEXE_sources += incflo_compute_forces.cpp
+CEXE_sources += incflo_explicit_update.cpp
 CEXE_sources += incflo_tagging.cpp
 CEXE_sources += incflo_regrid.cpp
 CEXE_sources += main.cpp

src/convection/incflo_compute_advection_term.cpp

@@ -28,9 +28,19 @@ void incflo::init_advection ()
  m_iconserv_density.resize(1, 1);
  m_iconserv_density_d.resize(1, 1);
 
- // We update (rho * tracer), not tracer itself, hence we update conservatively
+ // Advect scalars conservatively?
  m_iconserv_tracer.resize(m_ntrac, 1);
- m_iconserv_tracer_d.resize(m_ntrac, 1);
+ ParmParse pp("incflo");
+ pp.queryarr("trac_is_conservative", m_iconserv_tracer, 0, m_ntrac );
+ m_iconserv_tracer_d.resize(m_ntrac);
+ // copy
+#ifdef AMREX_USE_GPU
+ Gpu::htod_memcpy
+#else
+ std::memcpy
+#endif
+ (m_iconserv_tracer_d.data(), m_iconserv_tracer.data(), sizeof(int)*m_ntrac);
+
 }
 
 void
@@ -80,6 +90,21 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  Vector<Array<MultiFab*,AMREX_SPACEDIM> > fluxes(finest_level+1);
  Vector<Array<MultiFab*,AMREX_SPACEDIM> > faces(finest_level+1);
 
+ bool any_conserv_trac = false;
+ for (auto& i : m_iconserv_tracer){
+ if ( i == 1 ) {
+ any_conserv_trac = true;
+ break;
+ }
+ }
+ bool any_convective_trac = false;
+ for (auto& i : m_iconserv_tracer){
+ if ( i == 0 ) {
+ any_convective_trac = true;
+ break;
+ }
+ }
+
  for (int lev = 0; lev <= finest_level; ++lev) {
  AMREX_D_TERM(
  face_x[lev].define(u_mac[lev]->boxArray(),dmap[lev],n_flux_comp,0,MFInfo(),Factory(lev));,
@@ -94,7 +119,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  if (m_advect_momentum)
  rhovel[lev].define(vel[lev]->boxArray(),dmap[lev],AMREX_SPACEDIM,
  vel[lev]->nGrow(),MFInfo(),Factory(lev));
- if (m_advect_tracer && m_ntrac > 0)
+ if (m_advect_tracer && m_ntrac > 0 && any_conserv_trac)
  rhotrac[lev].define(vel[lev]->boxArray(),dmap[lev],tracer[lev]->nComp(),
  tracer[lev]->nGrow(),MFInfo(),Factory(lev));
 
@@ -151,7 +176,8 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  if (nghost_force() > 0)
  fillpatch_force(m_cur_time, vel_forces, nghost_force());
 
- // Note this is forcing for (rho s), not for s
+ // Note that for conservative tracers, this is forcing for (rho s)
+ // and for non-conservative, this is forcing for s
  if (m_advect_tracer)
  {
  compute_tra_forces(tra_forces, get_density_old_const());
@@ -409,7 +435,7 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  }
 
  // ************************************************************************
- // (Rho*Tracer)
+ // Tracer
  // ************************************************************************
  // Make a FAB holding (rho * tracer) that is the same size as the original tracer FAB
  if (m_advect_tracer && (m_ntrac>0)) {
@@ -421,13 +447,22 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
 
  Array4<Real const> tra = tracer[lev]->const_array(mfi);
  Array4<Real const> rho = density[lev]->const_array(mfi);
- Array4<Real > ro_trac = rhotrac[lev].array(mfi);
+ Array4<Real > trac_tmp;
 
- amrex::ParallelFor(bxg, m_ntrac,
- [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
- {
- ro_trac(i,j,k,n) = rho(i,j,k) * tra(i,j,k,n);
- });
+ auto const* iconserv = get_tracer_iconserv_device_ptr();
+ if ( any_conserv_trac ) {
+ trac_tmp = rhotrac[lev].array(mfi);
+
+ amrex::ParallelFor(bxg, m_ntrac,
+ [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
+ {
+ if ( iconserv[n] ){
+ trac_tmp(i,j,k,n) = rho(i,j,k) * tra(i,j,k,n);
+ } else {
+ trac_tmp(i,j,k,n) = tra(i,j,k,n);
+ }
+ });
+ }
 
  if (m_constant_density)
  face_comp = AMREX_SPACEDIM;
@@ -437,7 +472,8 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  is_velocity = false;
  Array4<int const> const& tracbc_arr = tracBC_MF ? (*tracBC_MF).const_array(mfi)
  : Array4<int const>{};
- HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi, ro_trac,
+ HydroUtils::ComputeFluxesOnBoxFromState( bx, ncomp, mfi,
+ any_conserv_trac ? trac_tmp : tracer[lev]->const_array(mfi),
  AMREX_D_DECL(flux_x[lev].array(mfi,face_comp),
  flux_y[lev].array(mfi,face_comp),
  flux_z[lev].array(mfi,face_comp)),
@@ -616,13 +652,10 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  } // mfi
  } // not constant density
 
- // Note: (rho*trac) is always updated conservatively -- we do not provide an option for
- // updating (rho*trac) convectively
  if (m_advect_tracer && m_ntrac > 0)
  {
  int flux_comp = (m_constant_density) ? AMREX_SPACEDIM : AMREX_SPACEDIM+1;
 
-//Was this OMP intentionally left off?
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
@@ -650,15 +683,31 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  (flagfab.getType(bx) != FabType::regular) ?
  ebfact->getBndryNormal().const_array(mfi) : Array4<Real const>{});
 #else
-  auto const& update_arr = conv_t[lev]->array(mfi);
-  HydroUtils::ComputeDivergence(bx, update_arr,
-  AMREX_D_DECL(flux_x[lev].const_array(mfi,flux_comp),
-  flux_y[lev].const_array(mfi,flux_comp),
-  flux_z[lev].const_array(mfi,flux_comp)),
-  m_ntrac, geom[lev], mult,
-  fluxes_are_area_weighted);
+ auto const& update_arr = conv_t[lev]->array(mfi);
+ HydroUtils::ComputeDivergence(bx, update_arr,
+ AMREX_D_DECL(flux_x[lev].const_array(mfi,flux_comp),
+ flux_y[lev].const_array(mfi,flux_comp),
+ flux_z[lev].const_array(mfi,flux_comp)),
+ m_ntrac, geom[lev], mult,
+ fluxes_are_area_weighted);
 #endif
 
+ if ( any_convective_trac )
+ {
+ // If convective, we define u dot grad trac = div (u trac) - trac div(u)
+ HydroUtils::ComputeConvectiveTerm(bx, m_ntrac, mfi,
+ tracer[lev]->array(mfi,0),
+ AMREX_D_DECL(face_x[lev].array(mfi),
+ face_y[lev].array(mfi),
+ face_z[lev].array(mfi)),
+ divu[lev].array(mfi),
+ update_arr,
+ get_tracer_iconserv_device_ptr(),
+#ifdef AMREX_USE_EB
+ *ebfact,
+#endif
+ m_advection_type);
+ }
  } // mfi
  } // advect tracer
 
@@ -699,7 +748,8 @@ incflo::compute_convective_term (Vector<MultiFab*> const& conv_u,
  if (m_advect_tracer) {
  auto const& bc_tra = get_tracer_bcrec_device_ptr();
  redistribute_term(mfi, *conv_t[lev], dtdt_tmp,
- rhotrac[lev], bc_tra, lev);
+ any_conserv_trac ? rhotrac[lev] : *tracer[lev],
+ bc_tra, lev);
  }
  } // mfi
 #endif

src/diffusion/DiffusionScalarOp.H

@@ -28,7 +28,7 @@ public:
 
  void compute_laps (amrex::Vector<amrex::MultiFab*> const& laps,
  amrex::Vector<amrex::MultiFab const*> const& a_scalar,
- amrex::Vector<amrex::MultiFab const*> const& density,
+ amrex::Vector<amrex::MultiFab const*> const& /*density*/,
  amrex::Vector<amrex::MultiFab const*> const& eta);
 
  void compute_divtau (amrex::Vector<amrex::MultiFab*> const& a_divtau,

src/diffusion/DiffusionScalarOp.cpp

@@ -434,7 +434,7 @@ DiffusionScalarOp::diffuse_vel_components (Vector<MultiFab*> const& vel,
 
 void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
  Vector<MultiFab const*> const& a_scalar,
- Vector<MultiFab const*> const& a_density,
+ Vector<MultiFab const*> const& /*a_density*/,
  Vector<MultiFab const*> const& a_eta)
 {
  BL_PROFILE("DiffusionScalarOp::compute_laps");
@@ -524,11 +524,6 @@ void DiffusionScalarOp::compute_laps (Vector<MultiFab*> const& a_laps,
  // We want to return div (mu grad)) phi
  m_reg_scal_apply_op->setScalars(0.0, -1.0);
 
- // This should have no effect since the first scalar is 0
- for (int lev = 0; lev <= finest_level; ++lev) {
- m_reg_scal_apply_op->setACoeffs(lev, *a_density[lev]);
- }
-
  for (int comp = 0; comp < m_incflo->m_ntrac; ++comp) {
 
  int eta_comp = comp;
@@ -651,11 +646,6 @@ void DiffusionScalarOp::compute_divtau (Vector<MultiFab*> const& a_divtau,
  // We want to return div (mu grad)) phi
  m_reg_vel_apply_op->setScalars(0.0, -1.0);
 
- // This should have no effect since the first scalar is 0
- for (int lev = 0; lev <= finest_level; ++lev) {
- m_reg_vel_apply_op->setACoeffs(lev, *a_density[lev]);
- }
-
  int eta_comp = 0;
  Vector<MultiFab> divtau_single;
  Vector<MultiFab> vel_single;

src/incflo.H

@@ -81,8 +81,6 @@ public:
  // Delete level data
  void ClearLevel (int lev) override;
 
-// for cuda
-
  void ComputeDt (int initialization, bool explicit_diffusion);
 
  amrex::Real vol_wgt_sum (amrex::Vector<amrex::MultiFab*> const& mf, int icomp);
@@ -115,10 +113,6 @@ public:
  std::string const& field);
  amrex::iMultiFab make_nodalBC_mask (int lev);
  amrex::Vector<amrex::MultiFab> make_robinBC_MFs(int lev, amrex::MultiFab* state = nullptr);
- // void make_ccBC_mask (int lev, const amrex::BoxArray& ba,
- // const amrex::DistributionMapping& dm);
- // void make_nodalBC_mask (int lev, const amrex::BoxArray& ba,
- // const amrex::DistributionMapping& dm);
 
 #ifdef AMREX_USE_EB
  void set_eb_velocity (int lev, amrex::Real time, amrex::MultiFab& eb_vel, int nghost);
@@ -156,6 +150,9 @@ public:
  amrex::Vector<amrex::MultiFab*> const& vel_forces,
  amrex::Real time);
 
+ void tracer_explicit_update(amrex::Vector<amrex::MultiFab> const& tra_forces);
+ void tracer_explicit_update_corrector(amrex::Vector<amrex::MultiFab> const& tra_forces);
+
  ///////////////////////////////////////////////////////////////////////////
  //
  // derive
@@ -204,6 +201,8 @@ public:
  [[nodiscard]] amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>
  average_scalar_eta_to_faces (int lev, int comp, amrex::MultiFab const& cc_eta) const;
 
+ // FIXME? compute_laps never used density in any meaningful way, but
+ // then removing density would break anything downstream....
  void compute_laps (amrex::Vector<amrex::MultiFab *> const& laps,
  amrex::Vector<amrex::MultiFab const*> const& scalar,
  amrex::Vector<amrex::MultiFab const*> const& density,

src/incflo.cpp

@@ -59,7 +59,7 @@ void incflo::InitData ()
  InitialRedistribution();
 #endif
 
- if (m_do_initial_proj) {
+  if (m_do_initial_proj) {
  InitialProjection();
  }
 

src/incflo_apply_corrector.cpp

@@ -166,13 +166,11 @@ void incflo::ApplyCorrector()
  // Define local variables for lambda to capture.
  // *************************************************************************************
  Real l_dt = m_dt;
- bool l_constant_density = m_constant_density;
- int l_ntrac = (m_advect_tracer) ? m_ntrac : 0;
 
  // *************************************************************************************
  // Update density first
  // *************************************************************************************
- if (l_constant_density)
+ if (m_constant_density)
  {
  for (int lev = 0; lev <= finest_level; lev++)
  MultiFab::Copy(density_nph[lev], m_leveldata[lev]->density_o, 0, 0, 1, 0);
@@ -229,80 +227,8 @@ void incflo::ApplyCorrector()
  // *************************************************************************************
  if (m_advect_tracer)
  {
- for (int lev = 0; lev <= finest_level; lev++)
- {
- auto& ld = *m_leveldata[lev];
-#ifdef _OPENMP
-#pragma omp parallel if (Gpu::notInLaunchRegion())
-#endif
- for (MFIter mfi(ld.tracer,TilingIfNotGPU()); mfi.isValid(); ++mfi)
- {
- Box const& bx = mfi.tilebox();
- Array4<Real const> const& tra_o = ld.tracer_o.const_array(mfi);
- Array4<Real const> const& rho_o = ld.density_o.const_array(mfi);
- Array4<Real > const& tra = ld.tracer.array(mfi);
- Array4<Real const> const& rho = ld.density.const_array(mfi);
- Array4<Real const> const& dtdt_o = ld.conv_tracer_o.const_array(mfi);
- Array4<Real const> const& dtdt = ld.conv_tracer.const_array(mfi);
- Array4<Real const> const& tra_f = (l_ntrac > 0) ? tra_forces[lev].const_array(mfi)
- : Array4<Real const>{};
-
- if (m_diff_type == DiffusionType::Explicit)
- {
- Array4<Real const> const& laps_o = (l_ntrac > 0) ? ld.laps_o.const_array(mfi)
- : Array4<Real const>{};
- Array4<Real const> const& laps = (l_ntrac > 0) ? ld.laps.const_array(mfi)
- : Array4<Real const>{};
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
- {
- for (int n = 0; n < l_ntrac; ++n)
- {
- Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * (
- m_half*( dtdt(i,j,k,n) + dtdt_o(i,j,k,n))
- +m_half*(laps_o(i,j,k,n) + laps(i,j,k,n))
- + tra_f(i,j,k,n) );
-
- tra_new /= rho(i,j,k);
- tra(i,j,k,n) = tra_new;
- }
- });
- }
- else if (m_diff_type == DiffusionType::Crank_Nicolson)
- {
- Array4<Real const> const& laps_o = (l_ntrac > 0) ? ld.laps_o.const_array(mfi)
- : Array4<Real const>{};
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
- {
- for (int n = 0; n < l_ntrac; ++n)
- {
- Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * (
- m_half*( dtdt(i,j,k,n) + dtdt_o(i,j,k,n))
- +m_half*(laps_o(i,j,k,n) )
- + tra_f(i,j,k,n) );
-
- tra_new /= rho(i,j,k);
- tra(i,j,k,n) = tra_new;
- }
- });
- }
- else if (m_diff_type == DiffusionType::Implicit)
- {
- amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
- {
- for (int n = 0; n < l_ntrac; ++n)
- {
- Real tra_new = rho_o(i,j,k)*tra_o(i,j,k,n) + l_dt * (
- m_half*( dtdt(i,j,k,n)+dtdt_o(i,j,k,n))
- + tra_f(i,j,k,n) );
-
- tra_new /= rho(i,j,k);
- tra(i,j,k,n) = tra_new;
- }
- });
- }
- } // mfi
- } // lev
- } // if (m_advect_tracer)
+ tracer_explicit_update_corrector(tra_forces);
+ }
 
  // *************************************************************************************
  // Solve diffusion equation for tracer