diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H
index fc3dd4060f..6b09ac1261 100644
--- a/amr-wind/incflo.H
+++ b/amr-wind/incflo.H
@@ -98,8 +98,7 @@ public:
void pre_advance_stage2();
void do_advance(int inonlin);
void post_overset_work();
- void advance();
- void advance_nonlinear();
+ void advance(int inonlin);
void prescribe_advance();
void post_advance_work();
@@ -135,8 +134,7 @@ public:
void ComputeDt(bool explicit_diffusion);
void ComputePrescribeDt();
- void ApplyPredictor(bool incremental_projection = false);
- void ApplyPredictorNonLinear(bool incremental_projection = false);
+ void ApplyPredictor(bool incremental_projection = false, int inonlin = 0);
void ApplyCorrector();
void ApplyPrescribeStep();
diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp
index 315bb1240e..3794f23cc1 100644
--- a/amr-wind/incflo.cpp
+++ b/amr-wind/incflo.cpp
@@ -318,14 +318,9 @@ void incflo::do_advance(int inonlin)
if (m_prescribe_vel && inonlin == 0) {
prescribe_advance();
} else {
- amrex::Print() << "INONLIN " << inonlin << std::endl;
- if (inonlin == 0) {
- advance();
- } else {
- amrex::Print() << "Iteration " << inonlin
- << " in advection iteration loop" << std::endl;
- advance_nonlinear();
- }
+ amrex::Print() << "Iteration " << inonlin
+ << " in advection iteration loop" << std::endl;
+ advance(inonlin);
}
if (m_sim.has_overset()) {
m_ovst_ops.post_advance_work();
diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp
index a9afdd0117..fb838f63ea 100644
--- a/amr-wind/incflo_advance.cpp
+++ b/amr-wind/incflo_advance.cpp
@@ -47,30 +47,18 @@ void incflo::pre_advance_stage2()
*
* \callgraph
*/
-void incflo::advance()
+void incflo::advance(int inonlin)
{
BL_PROFILE("amr-wind::incflo::Advance");
+ if (inonlin == 0) m_sim.pde_manager().advance_states();
- m_sim.pde_manager().advance_states();
-
- ApplyPredictor();
+ ApplyPredictor(false, inonlin);
if (!m_use_godunov) {
ApplyCorrector();
}
}
-void incflo::advance_nonlinear()
-{
- BL_PROFILE("amr_wind::incflo::Advance_Nonlinear");
-
- ApplyPredictorNonLinear();
-
- if (!m_use_godunov) {
- amrex::Abort("Non-linear iterations are not supported for MOL");
- }
-}
-
// Apply predictor step
//
// For Godunov, this completes the timestep. For MOL, this is the first part of
@@ -176,10 +164,9 @@ void incflo::advance_nonlinear()
*
\ref incflo::ApplyProjection "Apply projection"
*
*/
-void incflo::ApplyPredictor(bool incremental_projection)
+void incflo::ApplyPredictor(bool incremental_projection, int inonlin)
{
BL_PROFILE("amr-wind::incflo::ApplyPredictor");
-
// We use the new time value for things computed on the "*" state
Real new_time = m_time.new_time();
@@ -201,6 +188,16 @@ void incflo::ApplyPredictor(bool incremental_projection)
const auto& density_old = density_new.state(amr_wind::FieldState::Old);
auto& density_nph = density_new.state(amr_wind::FieldState::NPH);
+ if (inonlin > 0) {
+ icns().fields().field.fillpatch(m_time.current_time());
+ // Get n + 1/2 velocity
+ amr_wind::field_ops::lincomb(
+ icns().fields().field.state(amr_wind::FieldState::NPH), 0.5,
+ icns().fields().field.state(amr_wind::FieldState::Old), 0, 0.5,
+ icns().fields().field, 0, 0, icns().fields().field.num_comp(),
+ icns().fields().field.num_grow());
+ }
+
// *************************************************************************************
// Compute viscosity / diffusive coefficients
// *************************************************************************************
@@ -279,7 +276,11 @@ void incflo::ApplyPredictor(bool incremental_projection)
}
// Extrapolate and apply MAC projection for advection velocities
- icns().pre_advection_actions(amr_wind::FieldState::Old);
+ if (inonlin == 0) {
+ icns().pre_advection_actions(amr_wind::FieldState::Old);
+ } else {
+ icns().pre_advection_actions(amr_wind::FieldState::NPH);
+ }
// For scalars only first
// *************************************************************************************
@@ -347,7 +348,11 @@ void incflo::ApplyPredictor(bool incremental_projection)
}
// With scalars computed, compute advection of momentum
- icns().compute_advection_term(amr_wind::FieldState::Old);
+ if (inonlin == 0) {
+ icns().compute_advection_term(amr_wind::FieldState::Old);
+ } else {
+ icns().compute_advection_term(amr_wind::FieldState::NPH);
+ }
// *************************************************************************************
// Define (or if use_godunov, re-define) the forcing terms and viscous terms
@@ -408,273 +413,24 @@ void incflo::ApplyPredictor(bool incremental_projection)
if (m_verbose > 2) {
PrintMaxVelLocations("after nodal projection");
}
-}
-
-// Predictor step for subsequent non-linear iterations within each timestep
-void incflo::ApplyPredictorNonLinear(bool incremental_projection)
-{
- BL_PROFILE("amr-wind::incflo::ApplyPredictorNonLinear");
-
- // We use the new time value for things computed on the "*" state
- Real new_time = m_time.new_time();
-
- if (m_verbose > 2) {
- PrintMaxValues("before non-linear predictor step");
- }
-
- if (m_use_godunov) {
- amr_wind::io::print_mlmg_header("Godunov:");
- } else {
- amr_wind::io::print_mlmg_header("Predictor:");
- }
-
- auto& icns_fields = icns().fields();
- auto& velocity_new = icns_fields.field;
- auto& velocity_old = velocity_new.state(amr_wind::FieldState::Old);
-
- auto& density_new = density();
- const auto& density_old = density_new.state(amr_wind::FieldState::Old);
- auto& density_nph = density_new.state(amr_wind::FieldState::NPH);
-
- icns().fields().field.fillpatch(m_time.current_time());
- // Get n + 1/2 velocity
- amr_wind::field_ops::lincomb(
- icns().fields().field.state(amr_wind::FieldState::NPH), 0.5,
- icns().fields().field.state(amr_wind::FieldState::Old), 0, 0.5,
- icns().fields().field, 0, 0, icns().fields().field.num_comp(),
- icns().fields().field.num_grow());
-
- // *************************************************************************************
- // Compute viscosity / diffusive coefficients
- // *************************************************************************************
- // TODO: This sub-section has not been adjusted for mesh mapping - adjust in
- // corrector too
- m_sim.turbulence_model().update_turbulent_viscosity(
- amr_wind::FieldState::Old, m_diff_type);
- icns().compute_mueff(amr_wind::FieldState::Old);
- for (auto& eqns : scalar_eqns()) {
- eqns->compute_mueff(amr_wind::FieldState::Old);
- }
-
- // *************************************************************************************
- // Define the forcing terms to use in the Godunov prediction
- // *************************************************************************************
- // TODO: Godunov has not been adjusted for mesh mapping - adjust in
- // corrector too
- if (m_use_godunov) {
- icns().compute_source_term(amr_wind::FieldState::Old);
- for (auto& seqn : scalar_eqns()) {
- seqn->compute_source_term(amr_wind::FieldState::Old);
- }
- }
-
- // TODO: This sub-section has not been adjusted for mesh mapping - adjust in
- // corrector too
- if (need_divtau()) {
- // *************************************************************************************
- // Compute explicit viscous term using old density (1/rho)
- // *************************************************************************************
- // Reuse existing buffer to avoid creating new multifabs
- amr_wind::field_ops::copy(
- velocity_new, velocity_old, 0, 0, velocity_new.num_comp(), 1);
- icns().compute_diffusion_term(amr_wind::FieldState::Old);
- if (m_use_godunov) {
- auto& velocity_forces = icns_fields.src_term;
- // only the old states are used in predictor
- const auto& divtau = icns_fields.diff_term;
-
- amr_wind::field_ops::add(
- velocity_forces, divtau, 0, 0, AMREX_SPACEDIM, 0);
- }
- // *************************************************************************************
- // Compute explicit diffusive terms
- // *************************************************************************************
- for (auto& eqn : scalar_eqns()) {
- auto& field = eqn->fields().field;
- // Reuse existing buffer to avoid creating new multifabs
- amr_wind::field_ops::copy(
- field, field.state(amr_wind::FieldState::Old), 0, 0,
- field.num_comp(), 1);
-
- eqn->compute_diffusion_term(amr_wind::FieldState::Old);
-
- if (m_use_godunov) {
- amr_wind::field_ops::add(
- eqn->fields().src_term, eqn->fields().diff_term, 0, 0,
- field.num_comp(), 0);
- }
- }
- }
-
- if (m_use_godunov) {
- // ICNS source term (and viscous term) are in terms of momentum;
- // convert to velocity for MAC velocities by dividing by density
- amr_wind::field_ops::divide(
- icns().fields().src_term, density_old, 0, 0, 1, AMREX_SPACEDIM, 0);
-
- const int nghost_force = 1;
- IntVect ng(nghost_force);
- icns().fields().src_term.fillpatch(m_time.current_time(), ng);
-
- for (auto& eqn : scalar_eqns()) {
- eqn->fields().src_term.fillpatch(m_time.current_time(), ng);
- }
- }
-
- // Extrapolate and apply MAC projection for advection velocities
- icns().pre_advection_actions(amr_wind::FieldState::NPH);
-
- // For scalars only first
- // *************************************************************************************
- // if ( m_use_godunov) Compute the explicit advective terms
- // R_u^(n+1/2), R_s^(n+1/2) and R_t^(n+1/2)
- // if (!m_use_godunov) Compute the explicit advective terms
- // R_u^n , R_s^n and R_t^n
- // *************************************************************************************
- // TODO: Advection computation for scalar equations have not been adjusted
- // for mesh mapping
- for (auto& seqn : scalar_eqns()) {
- seqn->compute_advection_term(amr_wind::FieldState::Old);
- }
-
- // *************************************************************************************
- // Update density first
- // *************************************************************************************
- if (m_sim.pde_manager().constant_density()) {
- amr_wind::field_ops::copy(density_nph, density_old, 0, 0, 1, 1);
- }
-
- // TODO: This sub-section has not been adjusted for mesh mapping - adjust in
- // corrector too.
- // Perform scalar update one at a time. This is to allow an
- // updated density at `n+1/2` to be computed before other scalars use it
- // when computing their source terms.
- for (auto& eqn : scalar_eqns()) {
- // Compute (recompute for Godunov) the scalar forcing terms
- eqn->compute_source_term(amr_wind::FieldState::NPH);
-
- // Update the scalar (if explicit), or the RHS for implicit/CN
- eqn->compute_predictor_rhs(m_diff_type);
-
- auto& field = eqn->fields().field;
- if (m_diff_type != DiffusionType::Explicit) {
- amrex::Real dt_diff = (m_diff_type == DiffusionType::Implicit)
- ? m_time.delta_t()
- : 0.5 * m_time.delta_t();
-
- // Solve diffusion eqn. and update of the scalar field
- eqn->solve(dt_diff);
-
- // Post-processing actions after a PDE solve
- } else if (m_diff_type == DiffusionType::Explicit && m_use_godunov) {
- // explicit RK2
- std::unique_ptr diff_old =
- m_repo.create_scratch_field(1, 0, amr_wind::FieldLoc::CELL);
- auto& diff_new =
- eqn->fields().diff_term.state(amr_wind::FieldState::New);
- amr_wind::field_ops::copy(*diff_old, diff_new, 0, 0, 1, 0);
- eqn->compute_diffusion_term(amr_wind::FieldState::New);
- amrex::Real dto2 = 0.5 * m_time.delta_t();
- amr_wind::field_ops::saxpy(
- eqn->fields().field, -dto2, *diff_old, 0, 0, 1, 0);
- amr_wind::field_ops::saxpy(
- eqn->fields().field, +dto2, diff_new, 0, 0, 1, 0);
- }
- eqn->post_solve_actions();
+ // ScratchField to store the old np1
+ if (inonlin > 0) {
+ auto vel_np1_old = m_repo.create_scratch_field(
+ "vel_np1_old", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL);
- // Update scalar at n+1/2
+ auto vel_diff = m_repo.create_scratch_field(
+ "vel_diff", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL);
+ // lincomb to get old np1
amr_wind::field_ops::lincomb(
- field.state(amr_wind::FieldState::NPH), 0.5,
- field.state(amr_wind::FieldState::Old), 0, 0.5, field, 0, 0,
- field.num_comp(), 1);
- }
-
- // Compute advection of momentum at NPH
- icns().compute_advection_term(amr_wind::FieldState::NPH);
-
- // *************************************************************************************
- // Define (or if use_godunov, re-define) the forcing terms and viscous terms
- // independently for the right hand side, without 1/rho term
- // *************************************************************************************
- icns().compute_source_term(amr_wind::FieldState::NPH);
- icns().compute_diffusion_term(amr_wind::FieldState::Old);
-
- // *************************************************************************************
- // Evaluate right hand side and store in velocity
- // *************************************************************************************
- icns().compute_predictor_rhs(m_diff_type);
-
- // Recompute source and diffusion terms (if requested?)
- /* icns().compute_source_term(amr_wind::FieldState::NPH)
- icns().compute_diffusion_term(amr_wind::FieldState::NPH)*/
-
- // RHS
- icns().compute_predictor_rhs(m_diff_type);
- //
- if (m_verbose > 2) {
- PrintMaxVelLocations("after predictor rhs");
- }
-
- // *************************************************************************************
- // Solve diffusion equation for u* but using eta_old at old time
- // *************************************************************************************
- if (m_diff_type == DiffusionType::Crank_Nicolson ||
- m_diff_type == DiffusionType::Implicit) {
- Real dt_diff = (m_diff_type == DiffusionType::Implicit)
- ? m_time.delta_t()
- : 0.5 * m_time.delta_t();
- icns().solve(dt_diff);
- } else if (m_diff_type == DiffusionType::Explicit && m_use_godunov) {
- // explicit RK2
- std::unique_ptr diff_old =
- m_repo.create_scratch_field(
- AMREX_SPACEDIM, 0, amr_wind::FieldLoc::CELL);
-
- auto& diff_new =
- icns().fields().diff_term.state(amr_wind::FieldState::New);
- amr_wind::field_ops::copy(*diff_old, diff_new, 0, 0, AMREX_SPACEDIM, 0);
- icns().compute_diffusion_term(amr_wind::FieldState::New);
- amrex::Real dto2 = 0.5 * m_time.delta_t();
- amr_wind::field_ops::saxpy(
- icns().fields().field, -dto2, *diff_old, 0, 0, AMREX_SPACEDIM, 0);
- amr_wind::field_ops::saxpy(
- icns().fields().field, +dto2, diff_new, 0, 0, AMREX_SPACEDIM, 0);
- }
- icns().post_solve_actions();
+ (*vel_np1_old), -1.0,
+ icns().fields().field.state(amr_wind::FieldState::Old), 0, 2,
+ icns().fields().field.state(amr_wind::FieldState::NPH), 0, 0,
+ icns().fields().field.num_comp(), 1);
- if (m_verbose > 2) {
- PrintMaxVelLocations("after diffusion solve");
+ amr_wind::io::print_nonlinear_residual(m_sim, *vel_diff, *vel_np1_old);
+ vel_np1_old.reset();
+ vel_diff.reset();
}
-
- // ************************************************************************************
- //
- // Project velocity field, update pressure
- //
- // ************************************************************************************
- ApplyProjection(
- (density_new).vec_const_ptrs(), new_time, m_time.delta_t(),
- incremental_projection);
-
- if (m_verbose > 2) {
- PrintMaxVelLocations("after nodal projection");
- }
-
- // ScratchField to store the old np1
- auto vel_np1_old = m_repo.create_scratch_field(
- "vel_np1_old", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL);
-
- auto vel_diff = m_repo.create_scratch_field(
- "vel_diff", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL);
- // lincomb to get old np1
- amr_wind::field_ops::lincomb(
- (*vel_np1_old), -1.0,
- icns().fields().field.state(amr_wind::FieldState::Old), 0, 2,
- icns().fields().field.state(amr_wind::FieldState::NPH), 0, 0,
- icns().fields().field.num_comp(), 1);
-
- amr_wind::io::print_nonlinear_residual(m_sim, *vel_diff, *vel_np1_old);
- vel_np1_old.reset();
- vel_diff.reset();
}
//
// Apply corrector:
diff --git a/amr-wind/setup/init.cpp b/amr-wind/setup/init.cpp
index b669a8b2af..a5ed713735 100644
--- a/amr-wind/setup/init.cpp
+++ b/amr-wind/setup/init.cpp
@@ -103,7 +103,7 @@ void incflo::InitialIterations()
amrex::Print() << "In initial_iterations: iter = " << iter << "\n";
}
- ApplyPredictor(true);
+ ApplyPredictor(true,0);
{
auto& vel = icns().fields().field;