Skip to content

Commit

Permalink
fix corrector
Browse files Browse the repository at this point in the history
  • Loading branch information
@asalmgren
asalmgren committed Jun 16, 2024
1 parent 4844441 commit 5d5ac55
Show file tree
Hide file tree
Showing 4 changed files with 160 additions and 103 deletions.
25 changes: 11 additions & 14 deletions src/incflo_apply_corrector.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,8 @@ void incflo::ApplyCorrector()
{
BL_PROFILE("incflo::ApplyCorrector");

constexpr Real m_half = Real(0.5);

// We use the new time value for things computed on the "*" state
Real new_time = m_cur_time + m_dt;

Expand Down Expand Up @@ -97,21 +99,19 @@ void incflo::ApplyCorrector()
// in constructing the advection term
// **********************************************************************************************
Vector<MultiFab> vel_forces, tra_forces;
Vector<MultiFab> vel_eta , tra_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));
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) {
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) {
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 @@ -135,10 +135,7 @@ void incflo::ApplyCorrector()
// *************************************************************************************
// Compute viscosity / diffusive coefficients
// *************************************************************************************
compute_viscosity(GetVecOfPtrs(vel_eta),
get_density_new(), get_velocity_new(),
new_time, 1);
compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1);
compute_viscosity(GetVecOfPtrs(vel_eta), get_density_new(), get_velocity_new(), new_time, 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 @@ -164,8 +161,8 @@ void incflo::ApplyCorrector()
update_velocity(StepType::Corrector, vel_eta, vel_forces);

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

Expand Down
17 changes: 8 additions & 9 deletions src/incflo_apply_predictor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,17 +119,12 @@ void incflo::ApplyPredictor (bool incremental_projection)
}
}

// *************************************************************************************
// We now define the forcing terms to use in the Godunov prediction inside the predictor
// *************************************************************************************

// *************************************************************************************
// Compute viscosity / diffusive coefficients
// *************************************************************************************
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
Expand All @@ -142,12 +137,15 @@ void incflo::ApplyPredictor (bool incremental_projection)
}

// *************************************************************************************
// Compute explicit diffusive terms
// Compute explicit diffusive term -- note this is used inside compute_convective_term
// *************************************************************************************
if (m_advect_tracer && need_divtau())
if (m_advect_tracer)
{
compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(),
GetVecOfConstPtrs(tra_eta));
compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1);
if (need_divtau()) {
compute_laps(get_laps_old(), get_tracer_old_const(), get_density_old_const(),
GetVecOfConstPtrs(tra_eta));
}
}

// **********************************************************************************************
Expand All @@ -170,6 +168,7 @@ void incflo::ApplyPredictor (bool incremental_projection)
// Compute the explicit advective terms R_u^(n+1/2), R_s^(n+1/2) and R_t^(n+1/2)
// if (advection_type == "MOL" )
// Compute the explicit advective terms R_u^n , R_s^n and R_t^n
// Note that if advection_type != "MOL" then we call compute_tra_forces inside this routine
// *************************************************************************************
compute_convective_term(get_conv_velocity_old(), get_conv_density_old(), get_conv_tracer_old(),
get_velocity_old_const(), get_density_old_const(), get_tracer_old_const(),
Expand Down
19 changes: 16 additions & 3 deletions src/incflo_update_tracer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,10 @@ void incflo::update_tracer (StepType step_type, Vector<MultiFab>& tra_eta, Vecto

if (m_advect_tracer)
{
// *************************************************************************************
// Compute diffusion coefficient for tracer
// *************************************************************************************

// *************************************************************************************
// Compute the tracer forcing terms (forcing for (rho s), not for s)
// *************************************************************************************
Expand All @@ -18,11 +22,20 @@ void incflo::update_tracer (StepType step_type, Vector<MultiFab>& tra_eta, Vecto
// *************************************************************************************
// Compute explicit diffusive term (if corrector)
// *************************************************************************************
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));
if (step_type == StepType::Corrector)
{
compute_tracer_diff_coeff(GetVecOfPtrs(tra_eta),1);
if (m_diff_type == DiffusionType::Explicit) {
compute_laps(get_laps_new(), get_tracer_new_const(), get_density_new_const(),
GetVecOfConstPtrs(tra_eta));
}
}

// *************************************************************************************
// 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
// *************************************************************************************
if (step_type == StepType::Predictor) {
tracer_explicit_update(tra_forces);
} else if (step_type == StepType::Corrector) {
Expand Down
Loading

0 comments on commit 5d5ac55

Please sign in to comment.