Skip to content

Commit

Permalink
fix to previous commit
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Jun 15, 2024
1 parent ff6374b commit 4844441
Show file tree
Hide file tree
Showing 4 changed files with 50 additions and 35 deletions.
3 changes: 2 additions & 1 deletion src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -163,7 +163,8 @@ public:
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);
void update_tracer (StepType step_type, amrex::Vector<amrex::MultiFab>& tra_eta,
amrex::Vector<amrex::MultiFab>& tra_forces);
void update_velocity (StepType step_type, amrex::Vector<amrex::MultiFab>& vel_eta,
amrex::Vector<amrex::MultiFab>& vel_forces);

Expand Down
22 changes: 14 additions & 8 deletions src/incflo_apply_corrector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,17 +97,22 @@ void incflo::ApplyCorrector()
// in constructing the advection term
// **********************************************************************************************
Vector<MultiFab> vel_forces, tra_forces;
Vector<MultiFab> vel_eta;
Vector<MultiFab> vel_eta , tra_eta;

for (int lev = 0; lev <= finest_level; ++lev) {
vel_forces.emplace_back(grids[lev], dmap[lev], AMREX_SPACEDIM, nghost_force(),
MFInfo(), Factory(lev));
if (m_advect_tracer) {
tra_forces.emplace_back(grids[lev], dmap[lev], m_ntrac, nghost_force(),
MFInfo(), Factory(lev));
}
vel_eta.emplace_back(grids[lev], dmap[lev], 1, 1, MFInfo(), Factory(lev));
}

if (m_advect_tracer) {
for (int lev = 0; lev <= finest_level; ++lev) {
tra_forces.emplace_back(grids[lev], dmap[lev], m_ntrac, nghost_force(),
MFInfo(), Factory(lev));
tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev));
} //lev
} // if

// *************************************************************************************
// Compute the MAC-projected velocities at all levels
// *************************************************************************************
Expand All @@ -133,6 +138,7 @@ 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 @@ -143,14 +149,14 @@ void incflo::ApplyCorrector()
}

// *************************************************************************************
// Update density first
// Update density
// *************************************************************************************
update_density(StepType::Corrector);

// *************************************************************************************
// Update tracer next
// Update tracer
// *************************************************************************************
update_tracer(StepType::Corrector, tra_forces);
update_tracer(StepType::Corrector, tra_eta, tra_forces);

// *************************************************************************************
// Update velocity
Expand Down
41 changes: 29 additions & 12 deletions src/incflo_apply_predictor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -94,10 +94,13 @@ void incflo::ApplyPredictor (bool incremental_projection)
}
}

// *************************************************************************************
// Allocate space for half-time density
// *************************************************************************************
// Forcing terms
Vector<MultiFab> vel_forces, tra_forces;

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

// *************************************************************************************
// Allocate space for the forcing terms
Expand All @@ -111,6 +114,9 @@ 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 @@ -123,23 +129,38 @@ 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
// Compute the forcing terms
// *************************************************************************************
bool include_pressure_gradient = !(m_use_mac_phi_in_godunov);
compute_vel_forces(GetVecOfPtrs(vel_forces), get_velocity_old_const(),
get_density_old_const(), get_tracer_old_const(), get_tracer_old_const(),
include_pressure_gradient);

// **********************************************************************************************
// Compute the MAC-projected velocities at all levels
// *************************************************************************************
compute_MAC_projected_velocities(get_velocity_old_const(), get_density_old_const(),
AMREX_D_DECL(GetVecOfPtrs(u_mac), GetVecOfPtrs(v_mac),
GetVecOfPtrs(w_mac)), GetVecOfPtrs(vel_forces), m_cur_time);
Expand All @@ -162,20 +183,18 @@ void incflo::ApplyPredictor (bool incremental_projection)
// *************************************************************************************
update_density(StepType::Predictor);

// *************************************************************************************
// **********************************************************************************************
// Update tracer
// *************************************************************************************
update_tracer(StepType::Predictor, tra_forces);
// **********************************************************************************************
update_tracer(StepType::Predictor, tra_eta, tra_forces);

// *************************************************************************************
// **********************************************************************************************
// Update velocity
// *************************************************************************************
// **********************************************************************************************
update_velocity(StepType::Predictor, vel_eta, vel_forces);

// **********************************************************************************************
//
// Project velocity field, update pressure
//
// **********************************************************************************************
ApplyProjection(get_density_nph_const(),new_time,m_dt,incremental_projection);

Expand All @@ -191,9 +210,7 @@ void incflo::ApplyPredictor (bool incremental_projection)

#ifdef AMREX_USE_EB
// **********************************************************************************************
//
// Over-write velocity in cells with vfrac < 1e-4
//
// **********************************************************************************************
if (m_advection_type == "MOL")
incflo_correct_small_cells(get_velocity_new(),
Expand Down
19 changes: 5 additions & 14 deletions src/incflo_update_tracer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,32 +2,23 @@

using namespace amrex;

void incflo::update_tracer (StepType step_type, Vector<MultiFab>& tra_forces)
void incflo::update_tracer (StepType step_type, Vector<MultiFab>& tra_eta, Vector<MultiFab>& tra_forces)
{
BL_PROFILE("incflo::update_tracer");

Vector<MultiFab> tra_eta;

Real new_time = m_cur_time + m_dt;

if (m_advect_tracer)
{
// *************************************************************************************
// Compute the tracer forcing terms ( forcing for (rho s) if conservative )
// Compute the tracer forcing terms (forcing for (rho s), not for s)
// *************************************************************************************
for (int lev = 0; lev <= finest_level; ++lev) {
tra_eta.emplace_back(grids[lev], dmap[lev], m_ntrac, 1, MFInfo(), Factory(lev));
}

compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1);
compute_tra_forces(GetVecOfPtrs(tra_forces), get_density_nph_const());

// *************************************************************************************
// Compute explicit diffusive terms
// Compute explicit diffusive term (if corrector)
// *************************************************************************************
if (step_type == StepType::Predictor && need_divtau()) {
compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(),
GetVecOfConstPtrs(tra_eta));
} else if (step_type == StepType::Corrector && m_diff_type == DiffusionType::Explicit) {
if (step_type == StepType::Corrector && m_diff_type == DiffusionType::Explicit) {
compute_laps(get_laps_new(), get_tracer_new_const(), get_density_new_const(),
GetVecOfConstPtrs(tra_eta));
}
Expand Down

0 comments on commit 4844441

Please sign in to comment.