Skip to content

Commit

Permalink
separate routines into separate files so predictor and corrector aren…
Browse files Browse the repository at this point in the history 
…'t so long
  • Loading branch information
@asalmgren
asalmgren committed Jun 15, 2024
1 parent 844091a commit 0645f47
Show file tree
Hide file tree
Showing 7 changed files with 186 additions and 226 deletions.
4 changes: 3 additions & 1 deletion src/Make.package
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,10 @@ 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 += incflo_tagging.cpp
CEXE_sources += incflo_update_density.cpp
CEXE_sources += incflo_update_tracer.cpp
CEXE_sources += main.cpp

ifeq ($(USE_EB), TRUE)
Expand Down
12 changes: 10 additions & 2 deletions src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,10 @@
#include <DiffusionTensorOp.H>
#include <DiffusionScalarOp.H>

enum struct StepType {
Predictor, Corrector
};

class incflo : public amrex::AmrCore
{
public:
Expand Down Expand Up @@ -131,8 +135,8 @@ public:
//
///////////////////////////////////////////////////////////////////////////

void ApplyPredictor(bool incremental_projection = false);
void ApplyCorrector();
void ApplyPredictor (bool incremental_projection = false);
void ApplyCorrector ();
void compute_convective_term (amrex::Vector<amrex::MultiFab*> const& conv_u,
amrex::Vector<amrex::MultiFab*> const& conv_r,
amrex::Vector<amrex::MultiFab*> const& conv_t,
Expand All @@ -158,6 +162,9 @@ public:
void tracer_explicit_update(amrex::Vector<amrex::MultiFab> const& tra_forces);
void tracer_explicit_update_corrector(amrex::Vector<amrex::MultiFab> const& tra_forces);

void update_density (StepType step_type);
void update_tracer (StepType step_type, amrex::Vector<amrex::MultiFab>& tra_forces);

///////////////////////////////////////////////////////////////////////////
//
// derive
Expand Down Expand Up @@ -589,6 +596,7 @@ private:
amrex::MultiFab density;
amrex::MultiFab density_eb;
amrex::MultiFab density_o;
amrex::MultiFab density_nph;
amrex::MultiFab tracer;
amrex::MultiFab tracer_eb;
amrex::MultiFab tracer_o;
Expand Down
100 changes: 4 additions & 96 deletions src/incflo_apply_corrector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@ void incflo::ApplyCorrector()
// in constructing the advection term
// **********************************************************************************************
Vector<MultiFab> vel_forces, tra_forces;
Vector<MultiFab> vel_eta, tra_eta;
Vector<MultiFab> vel_eta;
for (int lev = 0; lev <= finest_level; ++lev) {
vel_forces.emplace_back(grids[lev], dmap[lev], AMREX_SPACEDIM, nghost_force(),
MFInfo(), Factory(lev));
Expand All @@ -116,9 +116,6 @@ void incflo::ApplyCorrector()
MFInfo(), Factory(lev));
}
vel_eta.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev));
if (m_advect_tracer) {
tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev));
}
}

// *************************************************************************************
Expand Down Expand Up @@ -146,7 +143,6 @@ void incflo::ApplyCorrector()
compute_viscosity(GetVecOfPtrs(vel_eta),
get_density_new(), get_velocity_new(),
new_time, 1);
compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1);

// Here we create divtau of the (n+1,*) state that was computed in the predictor;
// we use this laps only if DiffusionType::Explicit
Expand All @@ -156,11 +152,6 @@ void incflo::ApplyCorrector()
get_density_new_const(), GetVecOfConstPtrs(vel_eta));
}

if (m_advect_tracer && m_diff_type == DiffusionType::Explicit) {
compute_laps(get_laps_new(), get_tracer_new_const(),
get_density_new_const(), GetVecOfConstPtrs(tra_eta));
}

// *************************************************************************************
// Define local variables for lambda to capture.
// *************************************************************************************
Expand All @@ -169,95 +160,12 @@ void incflo::ApplyCorrector()
// *************************************************************************************
// Update density first
// *************************************************************************************
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);
} else {
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.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box const& bx = mfi.tilebox();
Array4<Real const> const& rho_o = ld.density_o.const_array(mfi);
Array4<Real> const& rho_n = ld.density.array(mfi);
Array4<Real> const& rho_nph = density_nph[lev].array(mfi);
Array4<Real const> const& drdt_o = ld.conv_density_o.const_array(mfi);
Array4<Real const> const& drdt = ld.conv_density.const_array(mfi);

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
const Real rho_old = rho_o(i,j,k);

Real rho_new = rho_old + l_dt * m_half*(drdt(i,j,k)+drdt_o(i,j,k));
rho_nph(i,j,k) = m_half * (rho_old + rho_new);

rho_n (i,j,k) = rho_new;
});
} // mfi
} // lev

// Average down solution
for (int lev = finest_level-1; lev >= 0; --lev) {
#ifdef AMREX_USE_EB
amrex::EB_average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density,
0, 1, refRatio(lev));
#else
amrex::average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density,
0, 1, refRatio(lev));
#endif
}
} // not constant density

// *************************************************************************************
// Compute the tracer forcing terms (forcing for (rho s), not for s)
// *************************************************************************************
if (m_advect_tracer)
compute_tra_forces(GetVecOfPtrs(tra_forces), GetVecOfConstPtrs(density_nph));
update_density(StepType::Corrector);

// *************************************************************************************
// Update the tracer next (note that dtdt already has rho in it)
// (rho trac)^new = (rho trac)^old + dt * (
// div(rho trac u) + div (mu grad trac) + rho * f_t
// Update tracer next
// *************************************************************************************
if (m_advect_tracer)
{
tracer_explicit_update_corrector(tra_forces);
}

// *************************************************************************************
// Solve diffusion equation for tracer
// *************************************************************************************
if ( m_advect_tracer )
{
if (m_diff_type == DiffusionType::Crank_Nicolson || m_diff_type == DiffusionType::Implicit)
{
const int ng_diffusion = 1;
for (int lev = 0; lev <= finest_level; ++lev)
fillphysbc_tracer(lev, new_time, m_leveldata[lev]->tracer, ng_diffusion);

Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : m_half*m_dt;
diffuse_scalar(get_tracer_new(), get_density_new(), GetVecOfConstPtrs(tra_eta), dt_diff);
}
else
{
// Need to average down tracer since the diffusion solver didn't do it for us.
for (int lev = finest_level-1; lev >= 0; --lev) {
#ifdef AMREX_USE_EB
amrex::EB_average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer,
0, m_ntrac, refRatio(lev));
#else
amrex::average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer,
0, m_ntrac, refRatio(lev));
#endif
}
}
} // if (m_advect_tracer)

update_tracer(StepType::Corrector, tra_forces);

// *************************************************************************************
// Define the forcing terms to use in the final update (using half-time density)
Expand Down
120 changes: 5 additions & 115 deletions src/incflo_apply_predictor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ void incflo::ApplyPredictor (bool incremental_projection)
// Forcing terms
Vector<MultiFab> vel_forces, tra_forces;

Vector<MultiFab> vel_eta, tra_eta;
Vector<MultiFab> vel_eta;

// *************************************************************************************
// Allocate space for the forcing terms
Expand All @@ -120,9 +120,6 @@ void incflo::ApplyPredictor (bool incremental_projection)
MFInfo(), Factory(lev));
}
vel_eta.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev));
if (m_advect_tracer) {
tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev));
}
}

// *************************************************************************************
Expand All @@ -135,27 +132,16 @@ void incflo::ApplyPredictor (bool incremental_projection)
compute_viscosity(GetVecOfPtrs(vel_eta),
get_density_old(), get_velocity_old(),
m_cur_time, 1);
compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1);

// *************************************************************************************
// Compute explicit viscous term
// Note that for !advect_momentum, this actually computes divtau / rho
// *************************************************************************************
if (need_divtau() || use_tensor_correction )
{
if (need_divtau() || use_tensor_correction ) {
compute_divtau(get_divtau_old(),get_velocity_old_const(),
get_density_old_const(),GetVecOfConstPtrs(vel_eta));
}

// *************************************************************************************
// Compute explicit diffusive terms
// *************************************************************************************
if (m_advect_tracer && need_divtau())
{
compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(),
GetVecOfConstPtrs(tra_eta));
}

// **********************************************************************************************
// Compute the MAC-projected velocities at all levels
// *************************************************************************************
Expand Down Expand Up @@ -188,107 +174,12 @@ void incflo::ApplyPredictor (bool incremental_projection)
// *************************************************************************************
// Update density first
// *************************************************************************************
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, 1);
} else {
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.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box const& bx = mfi.tilebox();
Array4<Real const> const& rho_o = ld.density_o.const_array(mfi);
Array4<Real> const& rho_new = ld.density.array(mfi);
Array4<Real const> const& drdt = ld.conv_density_o.const_array(mfi);

amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
rho_new(i,j,k) = rho_o(i,j,k) + l_dt * drdt(i,j,k);
});
} // mfi

// Fill ghost cells of the new density field so that we can define density_nph
// on the valid region grown by 1
int ng = 1;
fillpatch_density(lev, m_t_new[lev], ld.density, ng);

for (MFIter mfi(ld.velocity,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box const& gbx = mfi.growntilebox(1);
Array4<Real const> const& rho_old = ld.density_o.const_array(mfi);
Array4<Real const> const& rho_new = ld.density.const_array(mfi);
Array4<Real> const& rho_nph = density_nph[lev].array(mfi);

amrex::ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
rho_nph(i,j,k) = m_half * (rho_old(i,j,k) + rho_new(i,j,k));
});
} // mfi

} // lev

// Average down solution
for (int lev = finest_level-1; lev >= 0; --lev) {
#ifdef AMREX_USE_EB
amrex::EB_average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density,
0, 1, refRatio(lev));
#else
amrex::average_down(m_leveldata[lev+1]->density, m_leveldata[lev]->density,
0, 1, refRatio(lev));
#endif
}
} // not constant density


// *************************************************************************************
// Compute (or if Godunov, re-compute) the tracer forcing terms
// (forcing for (rho s) if conservative)
// *************************************************************************************
if (m_advect_tracer)
compute_tra_forces(GetVecOfPtrs(tra_forces), GetVecOfConstPtrs(density_nph));
update_density(StepType::Predictor);

// *************************************************************************************
// Update the tracer next
// Update tracer next
// *************************************************************************************
if (m_advect_tracer)
{
tracer_explicit_update(tra_forces);
}

// *************************************************************************************
// Solve diffusion equation for tracer
// *************************************************************************************
if ( m_advect_tracer )
{
if (m_diff_type == DiffusionType::Crank_Nicolson || m_diff_type == DiffusionType::Implicit)
{
const int ng_diffusion = 1;
for (int lev = 0; lev <= finest_level; ++lev)
fillphysbc_tracer(lev, new_time, m_leveldata[lev]->tracer, ng_diffusion);

Real dt_diff = (m_diff_type == DiffusionType::Implicit) ? m_dt : m_half*m_dt;
diffuse_scalar(get_tracer_new(), get_density_new(), GetVecOfConstPtrs(tra_eta), dt_diff);
}
else
{
// Need to average down tracer since the diffusion solver didn't do it for us.
for (int lev = finest_level-1; lev >= 0; --lev) {
#ifdef AMREX_USE_EB
amrex::EB_average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer,
0, m_ntrac, refRatio(lev));
#else
amrex::average_down(m_leveldata[lev+1]->tracer, m_leveldata[lev]->tracer,
0, m_ntrac, refRatio(lev));
#endif
}
}
} // if (m_advect_tracer)

update_tracer(StepType::Predictor, tra_forces);

// *************************************************************************************
// Define (or if advection_type != "MOL", re-define) the forcing terms, without the viscous terms
Expand All @@ -298,7 +189,6 @@ void incflo::ApplyPredictor (bool incremental_projection)
GetVecOfConstPtrs(density_nph),
get_tracer_old_const(), get_tracer_new_const());


// *************************************************************************************
// Update the velocity
// *************************************************************************************
Expand Down
Loading

0 comments on commit 0645f47

Please sign in to comment.