Skip to content

Commit

Permalink
Add nonconservative tracer (#94)
Browse files Browse the repository at this point in the history
  • Loading branch information
@cgilet
cgilet authored Apr 2, 2024
1 parent a732e84 commit 3c02cad
Show file tree
Hide file tree
Showing 12 changed files with 290 additions and 195 deletions.
1 change: 1 addition & 0 deletions src/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions src/Make.package
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
96 changes: 73 additions & 23 deletions src/convection/incflo_compute_advection_term.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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));,
Expand All @@ -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));

Expand Down Expand Up @@ -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());
Expand Down Expand Up @@ -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)) {
Expand All @@ -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;
Expand All @@ -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)),
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/diffusion/DiffusionScalarOp.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
12 changes: 1 addition & 11 deletions src/diffusion/DiffusionScalarOp.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down
11 changes: 5 additions & 6 deletions src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -199,6 +196,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,
Expand Down
2 changes: 1 addition & 1 deletion src/incflo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,7 @@ void incflo::InitData ()
InitialRedistribution();
#endif

if (m_do_initial_proj) {
if (m_do_initial_proj) {
InitialProjection();
}

Expand Down
80 changes: 3 additions & 77 deletions src/incflo_apply_corrector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 3c02cad

Please sign in to comment.