From 0f5d519ec1fb0e3be4d4c12d2fc10e52a0b39b92 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Tue, 24 Sep 2024 12:57:43 -0600 Subject: [PATCH] Removed the replicate functions advance_nonlinear() and ApplyPredictorNonLinear(), and moved their functionality into advance() and ApplyPredictor() --- amr-wind/incflo.H | 6 +- amr-wind/incflo.cpp | 11 +- amr-wind/incflo_advance.cpp | 320 +++++------------------------------- amr-wind/setup/init.cpp | 2 +- 4 files changed, 44 insertions(+), 295 deletions(-) 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;