From e21c549e7fc78eec5a4ec2956efa892dc4feaf55 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 18 Aug 2023 16:45:26 -0600 Subject: [PATCH 01/47] ability to iterate over advection step --- amr-wind/convection/Godunov.H | 1 + amr-wind/equation_systems/icns/icns.H | 2 +- .../equation_systems/icns/icns_advection.H | 2 ++ amr-wind/incflo_advance.cpp | 21 +++++++++++++++++++ 4 files changed, 25 insertions(+), 1 deletion(-) diff --git a/amr-wind/convection/Godunov.H b/amr-wind/convection/Godunov.H index dcbf32bdbb..7c1b84c22b 100644 --- a/amr-wind/convection/Godunov.H +++ b/amr-wind/convection/Godunov.H @@ -8,6 +8,7 @@ namespace godunov { +<<<<<<< HEAD enum class scheme { PLM, PPM, PPM_NOLIM, BDS, WENO_JS, WENOZ, MINMOD, UPWIND }; } // namespace godunov diff --git a/amr-wind/equation_systems/icns/icns.H b/amr-wind/equation_systems/icns/icns.H index 9337fa78c4..703cdbf5c4 100644 --- a/amr-wind/equation_systems/icns/icns.H +++ b/amr-wind/equation_systems/icns/icns.H @@ -43,7 +43,7 @@ struct ICNS : VectorTransport static constexpr bool has_diffusion = true; // No n+1/2 state for velocity for now - static constexpr bool need_nph_state = false; + static constexpr bool need_nph_state = true; }; } // namespace amr_wind::pde diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 835e46e293..eab59b1783 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -239,6 +239,8 @@ struct AdvectionOp const auto& src_term = fields.src_term; // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term; + // if fstate is NPH, then n and n+1 are known, need only spatial interp + bool only_spatial = fstate == FieldState::NPH; const auto& dof_field = fields.field.state(fstate); auto flux_x = diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 608851f5fe..dc2f6ba9f6 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -350,6 +350,27 @@ void incflo::ApplyPredictor(bool incremental_projection) // ************************************************************************************* icns().compute_predictor_rhs(m_diff_type); + const int n_adv_iters = 2; // <-- Remove later + for (int ia = 1; ia < n_adv_iters; ++ia) { + // Fillpatch the velocity + icns().fields().field.fillpatch(0.0); + // 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()); + + // Recompute advection of momentum + icns().compute_advection_term(amr_wind::FieldState::NPH); + + // 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"); } From 2405405e8feeffb3cbeba65d0873f470cbebf815 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 18 Aug 2023 21:14:41 -0600 Subject: [PATCH 02/47] central scheme --- amr-wind/convection/Godunov.H | 1 - amr-wind/convection/incflo_godunov_central.H | 139 ++++++++++++++++++ .../equation_systems/icns/icns_advection.H | 6 +- 3 files changed, 144 insertions(+), 2 deletions(-) create mode 100644 amr-wind/convection/incflo_godunov_central.H diff --git a/amr-wind/convection/Godunov.H b/amr-wind/convection/Godunov.H index 7c1b84c22b..dcbf32bdbb 100644 --- a/amr-wind/convection/Godunov.H +++ b/amr-wind/convection/Godunov.H @@ -8,7 +8,6 @@ namespace godunov { -<<<<<<< HEAD enum class scheme { PLM, PPM, PPM_NOLIM, BDS, WENO_JS, WENOZ, MINMOD, UPWIND }; } // namespace godunov diff --git a/amr-wind/convection/incflo_godunov_central.H b/amr-wind/convection/incflo_godunov_central.H new file mode 100644 index 0000000000..58d1e96699 --- /dev/null +++ b/amr-wind/convection/incflo_godunov_central.H @@ -0,0 +1,139 @@ +#ifndef GODUNOV_CENTRAL_H +#define GODUNOV_CENTRAL_H + +#include +#include + +/* This header file contains the inlined __host__ __device__ functions required + for the scalar advection routines for 3D Godunov. It also contains function + declarations for controlling host functions. */ + +namespace { + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void central( + const amrex::Real sm1, + const amrex::Real s0, + const amrex::Real sp1, + amrex::Real& dsm, + amrex::Real& dsp) +{ + // Calculate gradients on both sides + dsp = sp1 - s0; + dsm = s0 - sm1; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_central_bc( + const int n, + const amrex::Real sm1, + const amrex::Real s0, + const amrex::Real sp1, + amrex::Real& dsm, + amrex::Real& dsp, + const int bclo, + const int bchi, + const int domlo, + const int domhi) +{ + using namespace amrex; + + if (bclo == BCType::ext_dir || bclo == BCType::hoextrap) { + if (n == domlo) { + // Ensure that left-side slope is used unchanged + dsm = s0 - sm1; + } + } + + if (bchi == BCType::ext_dir || bchi == BCType::hoextrap) { + if (n == domhi) { + // Ensure that the right-side slope is used unchanged + dsp = sp1 - s0; + } + } +} + +// !!! No upwinding for central scheme !!! +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_central_x( + const int i, + const int j, + const int k, + const int n, + const amrex::Real dx, + amrex::Real& Im, + amrex::Real& Ip, + const amrex::Array4& S, + const amrex::BCRec& bc, + const int domlo, + const int domhi) +{ + using namespace amrex; + amrex::Real sm1 = S(i - 1, j, k, n); + amrex::Real s0 = S(i, j, k, n); + amrex::Real sp1 = S(i + 1, j, k, n); + amrex::Real dsp = 0.0; + amrex::Real dsm = 0.0; + central(sm1, s0, sp1, dsm, dsp); + Godunov_central_bc( + i, sm1, s0, sp1, dsm, dsp, bc.lo(0), bc.hi(0), domlo, domhi); + + // Interpolate to edges + Ip = s0 + 0.5 * dsp; + Im = s0 - 0.5 * dsm; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_central_y( + const int i, + const int j, + const int k, + const int n, + const amrex::Real dx, + amrex::Real& Im, + amrex::Real& Ip, + const amrex::Array4& S, + const amrex::BCRec& bc, + const int domlo, + const int domhi) +{ + using namespace amrex; + amrex::Real sm1 = S(i, j - 1, k, n); + amrex::Real s0 = S(i, j, k, n); + amrex::Real sp1 = S(i, j + 1, k, n); + amrex::Real dsp = 0.0; + amrex::Real dsm = 0.0; + central(sm1, s0, sp1, dsm, dsp); + Godunov_central_bc( + j, sm1, s0, sp1, dsm, dsp, bc.lo(1), bc.hi(1), domlo, domhi); + + Ip = s0 + 0.5 * dsp; + Im = s0 - 0.5 * dsm; +} + +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_central_z( + const int i, + const int j, + const int k, + const int n, + const amrex::Real dx, + amrex::Real& Im, + amrex::Real& Ip, + const amrex::Array4& S, + const amrex::BCRec& bc, + const int domlo, + const int domhi) +{ + using namespace amrex; + amrex::Real sm1 = S(i, j, k - 1, n); + amrex::Real s0 = S(i, j, k, n); + amrex::Real sp1 = S(i, j, k + 1, n); + amrex::Real dsp = 0.0; + amrex::Real dsm = 0.0; + central(sm1, s0, sp1, dsm, dsp); + Godunov_central_bc( + k, sm1, s0, sp1, dsm, dsp, bc.lo(2), bc.hi(2), domlo, domhi); + + Ip = s0 + 0.5 * dsp; + Im = s0 - 0.5 * dsm; +} + +} // namespace + +#endif diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index eab59b1783..9bea07ff96 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -120,6 +120,8 @@ struct AdvectionOp godunov_scheme = godunov::scheme::WENO_JS; } else if (amrex::toLower(godunov_type) == "weno_z") { godunov_scheme = godunov::scheme::WENOZ; + } else if (amrex::toLower(godunov_type) == "central") { + godunov_scheme = godunov::scheme::CENTRAL; } else { amrex::Abort( "Invalid godunov_type specified. For godunov_type select " @@ -240,7 +242,9 @@ struct AdvectionOp // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term; // if fstate is NPH, then n and n+1 are known, need only spatial interp - bool only_spatial = fstate == FieldState::NPH; + bool only_spatial = + (fstate == FieldState::NPH || + godunov_scheme == godunov::scheme::CENTRAL); const auto& dof_field = fields.field.state(fstate); auto flux_x = From 9243e25fa46bc7ebeb1f3230c76e9d4b26c9973b Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 18 Aug 2023 21:37:17 -0600 Subject: [PATCH 03/47] input argument for advection iterations --- amr-wind/incflo.H | 3 +++ amr-wind/incflo.cpp | 4 ++++ amr-wind/incflo_advance.cpp | 3 +-- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 6ef1fc65d2..1faaef5158 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -172,6 +172,9 @@ private: // Prescribe advection velocity bool m_prescribe_vel = false; + // Advection iterations every timestep + int m_adv_iters{1}; + //! number of cells on all levels including covered cells amrex::Long m_cell_count{-1}; diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index f537155c7a..f27b6d5e1d 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -372,6 +372,10 @@ void incflo::init_physics_and_pde() if (activate_overset) { m_sim.activate_overset(); } + + // Get number of advection iterations + pp.query("advection_iterations", m_adv_iters); + m_adv_iters = std::max(1, m_adv_iters); } auto& pde_mgr = m_sim.pde_manager(); diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index dc2f6ba9f6..577d5316a8 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -350,8 +350,7 @@ void incflo::ApplyPredictor(bool incremental_projection) // ************************************************************************************* icns().compute_predictor_rhs(m_diff_type); - const int n_adv_iters = 2; // <-- Remove later - for (int ia = 1; ia < n_adv_iters; ++ia) { + for (int ia = 1; ia < m_adv_iters; ++ia) { // Fillpatch the velocity icns().fields().field.fillpatch(0.0); // Get n + 1/2 velocity From 086649a49e9ebfe6ad94db3b5812d374b7e0ee47 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Thu, 24 Aug 2023 09:11:41 -0600 Subject: [PATCH 04/47] something to try: weno-z used for first iteration, then central --- amr-wind/equation_systems/icns/icns_advection.H | 2 ++ 1 file changed, 2 insertions(+) diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 9bea07ff96..ba7f95605e 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -122,6 +122,8 @@ struct AdvectionOp godunov_scheme = godunov::scheme::WENOZ; } else if (amrex::toLower(godunov_type) == "central") { godunov_scheme = godunov::scheme::CENTRAL; + } else if (amrex::toLower(godunov_type) == "weno_z_central") { + godunov_scheme = godunov::scheme::WENOZ_CENTRAL; } else { amrex::Abort( "Invalid godunov_type specified. For godunov_type select " From cb423143cf6e82f69b4539e1cbbb293fe4a29114 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Tue, 6 Aug 2024 11:37:59 -0600 Subject: [PATCH 05/47] Adding advance_nonlinear() for non-linear iterations --- amr-wind/incflo_advance.cpp | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 577d5316a8..34dabd522d 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -60,6 +60,19 @@ void incflo::advance() } } +void incflo::advance_nonlinear(int it_nonlinear) +{ + BL_PROFILE("amr_wind::incflo::Advance_Nonlinear"); + + if (it_nonlinear == 0) m_sim.pde_manager().advance_states(); + + 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 From 7a4f4321fd7037012e8ce293e1c33298a04c152a Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Mon, 12 Aug 2024 13:18:30 -0600 Subject: [PATCH 06/47] Adding compute_fluxes_spatial() to AdvOp_Godunov.H for scalar equations --- amr-wind/equation_systems/AdvOp_Godunov.H | 87 ++++++++++++----------- 1 file changed, 44 insertions(+), 43 deletions(-) diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index e0871b33fb..6daa8a301a 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -168,6 +168,8 @@ struct AdvectionOp< (godunov_scheme == godunov::scheme::WENO_JS)) { m_allow_inflow_on_outflow = true; } + amrex::Real dt_arg = (only_spatial) ? 0.0 : dt; + HydroUtils::ComputeFluxesOnBoxFromState( bx, PDE::ndim, mfi, (PDE::multiply_rho ? rhotrac : tra_arr), @@ -177,7 +179,7 @@ struct AdvectionOp< known_edge_state, u_mac(lev).const_array(mfi), v_mac(lev).const_array(mfi), w_mac(lev).const_array(mfi), divu, - src_term(lev).const_array(mfi), geom[lev], dt, + src_term(lev).const_array(mfi), geom[lev], dt_arg, dof_field.bcrec(), dof_field.bcrec_device().data(), iconserv.data(), godunov_use_ppm, godunov_use_forces_in_trans, is_velocity, @@ -186,58 +188,57 @@ struct AdvectionOp< } else { amrex::Abort("Invalid godunov scheme"); } - amrex::Gpu::streamSynchronize(); } + amrex::Gpu::streamSynchronize(); } + } - amrex::Vector> fluxes( - repo.num_active_levels()); - for (int lev = 0; lev < repo.num_active_levels(); ++lev) { - fluxes[lev][0] = &(*flux_x)(lev); - fluxes[lev][1] = &(*flux_y)(lev); - fluxes[lev][2] = &(*flux_z)(lev); - } + amrex::Vector> + fluxes(repo.num_active_levels()); + for (int lev = 0; lev < repo.num_active_levels(); ++lev) { + fluxes[lev][0] = &(*flux_x)(lev); + fluxes[lev][1] = &(*flux_y)(lev); + fluxes[lev][2] = &(*flux_z)(lev); + } - // In order to enforce conservation across coarse-fine boundaries we - // must be sure to average down the fluxes before we use them - for (int lev = repo.num_active_levels() - 1; lev > 0; --lev) { - amrex::IntVect rr = - geom[lev].Domain().size() / geom[lev - 1].Domain().size(); - amrex::average_down_faces( - GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1], rr, - geom[lev - 1]); - } + // In order to enforce conservation across coarse-fine boundaries we + // must be sure to average down the fluxes before we use them + for (int lev = repo.num_active_levels() - 1; lev > 0; --lev) { + amrex::IntVect rr = + geom[lev].Domain().size() / geom[lev - 1].Domain().size(); + amrex::average_down_faces( + GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1], rr, geom[lev - 1]); + } - for (int lev = 0; lev < repo.num_active_levels(); ++lev) { + for (int lev = 0; lev < repo.num_active_levels(); ++lev) { #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for (amrex::MFIter mfi(dof_field(lev), amrex::TilingIfNotGPU()); - mfi.isValid(); ++mfi) { - const auto& bx = mfi.tilebox(); - - HydroUtils::ComputeDivergence( - bx, conv_term(lev).array(mfi), (*flux_x)(lev).array(mfi), - (*flux_y)(lev).array(mfi), (*flux_z)(lev).array(mfi), - PDE::ndim, geom[lev], amrex::Real(-1.0), - fluxes_are_area_weighted); - } + for (amrex::MFIter mfi(dof_field(lev), amrex::TilingIfNotGPU()); + mfi.isValid(); ++mfi) { + const auto& bx = mfi.tilebox(); + + HydroUtils::ComputeDivergence( + bx, conv_term(lev).array(mfi), (*flux_x)(lev).array(mfi), + (*flux_y)(lev).array(mfi), (*flux_z)(lev).array(mfi), PDE::ndim, + geom[lev], amrex::Real(-1.0), fluxes_are_area_weighted); } } - - PDEFields& fields; - Field& density; - Field& u_mac; - Field& v_mac; - Field& w_mac; - amrex::Gpu::DeviceVector iconserv; - - godunov::scheme godunov_scheme = godunov::scheme::WENOZ; - std::string godunov_type{"weno_z"}; - const bool fluxes_are_area_weighted{false}; - bool godunov_use_forces_in_trans{false}; - bool m_allow_inflow_on_outflow{false}; - std::string advection_type{"Godunov"}; +} + +PDEFields& fields; +Field & density; +Field & u_mac; +Field & v_mac; +Field & w_mac; +amrex::Gpu::DeviceVector iconserv; + +godunov::scheme godunov_scheme = godunov::scheme::WENOZ; +std::string godunov_type{"weno_z"}; +const bool fluxes_are_area_weighted{false}; +bool godunov_use_forces_in_trans{false}; +bool m_allow_inflow_on_outflow{false}; +std::string advection_type{"Godunov"}; }; } // namespace amr_wind::pde From 5812dc040bcd2410d97a69225ee40140cdeac3fa Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Wed, 14 Aug 2024 16:41:43 -0600 Subject: [PATCH 07/47] advance_nonlinear() and ApplyPredictorNonLinear() for performing non-linear iterations using the Godunov scheme --- amr-wind/equation_systems/AdvOp_Godunov.H | 2 + amr-wind/incflo.H | 2 + amr-wind/incflo.cpp | 10 +- amr-wind/incflo_advance.cpp | 257 +++++++++++++++++++++- 4 files changed, 267 insertions(+), 4 deletions(-) diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 6daa8a301a..6a6fcc2fd8 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -139,6 +139,7 @@ struct AdvectionOp< rho_arr(i, j, k) * tra_arr(i, j, k, n); }); } +<<<<<<< HEAD if ((godunov_scheme == godunov::scheme::PPM) || (godunov_scheme == godunov::scheme::PPM_NOLIM) || @@ -161,6 +162,7 @@ struct AdvectionOp< limiter_type = PPM::WENOZ; } else if (godunov_scheme == godunov::scheme::WENO_JS) { limiter_type = PPM::WENO_JS; +======= } else { limiter_type = PPM::default_limiter; } diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 1faaef5158..493ee8dd2e 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -98,6 +98,7 @@ public: void pre_advance_stage2(); void do_advance(); void advance(); + void advance_nonlinear(); void prescribe_advance(); void post_advance_work(); @@ -134,6 +135,7 @@ public: void ComputePrescribeDt(); void ApplyPredictor(bool incremental_projection = false); + void ApplyPredictorNonLinear(bool incremental_projection = false); void ApplyCorrector(); void ApplyPrescribeStep(); diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index f27b6d5e1d..3b1fac51d9 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -317,7 +317,15 @@ void incflo::do_advance() if (m_prescribe_vel) { prescribe_advance(); } else { - advance(); + + for (int it_nl = 1; it_nl < m_adv_iters; ++it_nl) { + + if (it_nl == 1) { + advance(); + } else { + advance_nonlinear(); + } + } } 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 34dabd522d..a6a1560e7b 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -60,12 +60,10 @@ void incflo::advance() } } -void incflo::advance_nonlinear(int it_nonlinear) +void incflo::advance_nonlinear() { BL_PROFILE("amr_wind::incflo::Advance_Nonlinear"); - if (it_nonlinear == 0) m_sim.pde_manager().advance_states(); - ApplyPredictorNonLinear(); if (!m_use_godunov) { @@ -432,6 +430,259 @@ void incflo::ApplyPredictor(bool incremental_projection) } } +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 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); + + // ************************************************************************************* + // 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::Old); + + // 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(); + + // Update scalar at n+1/2 + 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); + } + + // With scalars computed, compute advection of momentum + icns().compute_advection_term(amr_wind::FieldState::Old); + + // ************************************************************************************* + // 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); + + for (int ia = 1; ia < m_adv_iters; ++ia) { + // Fillpatch the velocity + icns().fields().field.fillpatch(0.0); + // 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()); + + // Recompute advection of momentum + icns().compute_advection_term(amr_wind::FieldState::NPH); + + // 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(); + + if (m_verbose > 2) { + PrintMaxVelLocations("after diffusion solve"); + } + + // ************************************************************************************ + // + // 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"); + } +} // // Apply corrector: // From 49832b0ff55f8ec770d75c130e5186e16a521fda Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Thu, 15 Aug 2024 11:45:17 -0600 Subject: [PATCH 08/47] m_adv_iters loop removed from ApplyPredictorNonLinear() --- amr-wind/incflo_advance.cpp | 57 +++++++++++++------------------------ 1 file changed, 19 insertions(+), 38 deletions(-) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index a6a1560e7b..843d91f4f8 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -361,26 +361,6 @@ void incflo::ApplyPredictor(bool incremental_projection) // ************************************************************************************* icns().compute_predictor_rhs(m_diff_type); - for (int ia = 1; ia < m_adv_iters; ++ia) { - // Fillpatch the velocity - icns().fields().field.fillpatch(0.0); - // 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()); - - // Recompute advection of momentum - icns().compute_advection_term(amr_wind::FieldState::NPH); - - // 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"); } @@ -430,6 +410,7 @@ void incflo::ApplyPredictor(bool incremental_projection) } } +// Predictor step for subsequent non-linear iterations within each timestep void incflo::ApplyPredictorNonLinear(bool incremental_projection) { BL_PROFILE("amr-wind::incflo::ApplyPredictorNonLinear"); @@ -438,7 +419,7 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) Real new_time = m_time.new_time(); if (m_verbose > 2) { - PrintMaxValues("before predictor step"); + PrintMaxValues("before non-linear predictor step"); } if (m_use_godunov) { @@ -615,26 +596,26 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) // ************************************************************************************* icns().compute_predictor_rhs(m_diff_type); - for (int ia = 1; ia < m_adv_iters; ++ia) { - // Fillpatch the velocity - icns().fields().field.fillpatch(0.0); - // 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()); + // Interpolate the NPH velocity from Old and New states + // Fillpatch the velocity + icns().fields().field.fillpatch(0.0); + // 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()); - // Recompute advection of momentum - icns().compute_advection_term(amr_wind::FieldState::NPH); + // Recompute advection of momentum + icns().compute_advection_term(amr_wind::FieldState::NPH); - // Recompute source and diffusion terms (if requested?) - /* icns().compute_source_term(amr_wind::FieldState::NPH) - icns().compute_diffusion_term(amr_wind::FieldState::NPH)*/ + // 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); - } + // RHS + icns().compute_predictor_rhs(m_diff_type); + // if (m_verbose > 2) { PrintMaxVelLocations("after predictor rhs"); } From b38a38cab393549ec39ad914d5c7b543ad921589 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Tue, 20 Aug 2024 17:44:59 -0600 Subject: [PATCH 09/47] - Flux prediction and preadvection routines without time extrapolation - Calculation of NPH velocity within the non-linear iteration loop --- amr-wind/convection/Godunov.H | 1 + amr-wind/equation_systems/AdvOp_Godunov.H | 6 + amr-wind/equation_systems/AdvOp_MOL.H | 6 + amr-wind/equation_systems/PDE.H | 8 + amr-wind/equation_systems/PDEBase.H | 3 + .../equation_systems/icns/icns_advection.H | 205 ++++++++++++++++++ amr-wind/equation_systems/vof/vof_advection.H | 6 + amr-wind/incflo_advance.cpp | 31 ++- 8 files changed, 250 insertions(+), 16 deletions(-) diff --git a/amr-wind/convection/Godunov.H b/amr-wind/convection/Godunov.H index dcbf32bdbb..a93df095c4 100644 --- a/amr-wind/convection/Godunov.H +++ b/amr-wind/convection/Godunov.H @@ -10,6 +10,7 @@ namespace godunov { enum class scheme { PLM, PPM, PPM_NOLIM, BDS, WENO_JS, WENOZ, MINMOD, UPWIND }; + } // namespace godunov #endif /* Godunov_H */ diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 6a6fcc2fd8..64b940aa27 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -80,6 +80,12 @@ struct AdvectionOp< const amrex::Real /*unused*/) {} + void preadvect_spatial( + const FieldState /*unused*/, + const amrex::Real /*unused*/, + const amrex::Real /*unused*/) + {} + void operator()(const FieldState fstate, const amrex::Real dt) { static_assert( diff --git a/amr-wind/equation_systems/AdvOp_MOL.H b/amr-wind/equation_systems/AdvOp_MOL.H index b276e39d78..97fc278ab1 100644 --- a/amr-wind/equation_systems/AdvOp_MOL.H +++ b/amr-wind/equation_systems/AdvOp_MOL.H @@ -43,6 +43,12 @@ struct AdvectionOp< const amrex::Real /*unused*/) {} + void preadvect_spatial( + const FieldState /*unused*/, + const amrex::Real /*unused*/, + const amrex::Real /*unused*/) + {} + void operator()(const FieldState fstate, const amrex::Real /*unused*/) { static_assert( diff --git a/amr-wind/equation_systems/PDE.H b/amr-wind/equation_systems/PDE.H index 485616d680..f74526c3ac 100644 --- a/amr-wind/equation_systems/PDE.H +++ b/amr-wind/equation_systems/PDE.H @@ -140,6 +140,14 @@ public: m_adv_op->preadvect(fstate, m_time.delta_t(), m_time.current_time()); } + void nonlinear_advection_actions(const FieldState fstate) override + { + BL_PROFILE( + "amr-wind::" + this->identifier() + "::pre_advection_actions"); + m_adv_op->preadvect_spatial( + fstate, m_time.delta_t(), m_time.current_time()); + } + void compute_predictor_rhs(const DiffusionType difftype) override { BL_PROFILE( diff --git a/amr-wind/equation_systems/PDEBase.H b/amr-wind/equation_systems/PDEBase.H index 5abbd19278..19206026c9 100644 --- a/amr-wind/equation_systems/PDEBase.H +++ b/amr-wind/equation_systems/PDEBase.H @@ -74,6 +74,9 @@ public: //! Perform necessary steps to prepare for advection calculations virtual void pre_advection_actions(const FieldState fstate) = 0; + //! Perform necessary steps to prepare for nonlinear advection calculations + virtual void nonlinear_advection_actions(const FieldState fstate) = 0; + //! Compute the time derivative and advective term for the PDE system virtual void compute_advection_term(const FieldState fstate) = 0; diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index ba7f95605e..3b72339678 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -235,6 +235,205 @@ struct AdvectionOp } } + void preadvect_spatial( + const FieldState fstate, const amrex::Real dt, const amrex::Real time) + { + + const auto& repo = fields.repo; + const auto& geom = repo.mesh().Geom(); + + const auto& src_term = fields.src_term; + const auto& dof_field = fields.field.state(fstate); + auto bcrec_device = dof_field.bcrec_device(); + + for (int lev = 0; lev < repo.num_active_levels(); ++lev) { + u_mac(lev).setBndry(0.0); + v_mac(lev).setBndry(0.0); + w_mac(lev).setBndry(0.0); + } + + // + // Predict + // + for (int lev = 0; lev < repo.num_active_levels(); ++lev) { + if ((godunov_scheme == godunov::scheme::PPM_NOLIM) || + (godunov_scheme == godunov::scheme::WENOJS) || + (godunov_scheme == godunov::scheme::WENOZ) || + (godunov_scheme == godunov::scheme::CENTRAL) || + (godunov_scheme == godunov::scheme::WENOZ_CENTRAL)) { +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + { + amrex::FArrayBox scratch; + for (amrex::MFIter mfi( + dof_field(lev), amrex::TilingIfNotGPU()); + mfi.isValid(); ++mfi) { + amrex::Box const& bx = mfi.tilebox(); + amrex::Box const& bxg1 = amrex::grow(bx, 1); + amrex::Box const& xbx = mfi.nodaltilebox(0); + amrex::Box const& ybx = mfi.nodaltilebox(1); + amrex::Box const& zbx = mfi.nodaltilebox(2); + + amrex::Array4 const& a_umac = + u_mac(lev).array(mfi); + amrex::Array4 const& a_vmac = + v_mac(lev).array(mfi); + amrex::Array4 const& a_wmac = + w_mac(lev).array(mfi); + amrex::Array4 const& a_vel = + dof_field(lev).const_array(mfi); + amrex::Array4 const& a_f = + src_term(lev).const_array(mfi); + + scratch.resize(bxg1, ICNS::ndim * 12 + 3); + // Elixir eli = scratch.elixir(); // not + // needed because of streamSynchronize + // later + amrex::Real* p = scratch.dataPtr(); + + amrex::Array4 Imx = + makeArray4(p, bxg1, ICNS::ndim); + p += Imx.size(); + amrex::Array4 Ipx = + makeArray4(p, bxg1, ICNS::ndim); + p += Ipx.size(); + amrex::Array4 Imy = + makeArray4(p, bxg1, ICNS::ndim); + p += Imy.size(); + amrex::Array4 Ipy = + makeArray4(p, bxg1, ICNS::ndim); + p += Ipy.size(); + amrex::Array4 Imz = + makeArray4(p, bxg1, ICNS::ndim); + p += Imz.size(); + amrex::Array4 Ipz = + makeArray4(p, bxg1, ICNS::ndim); + p += Ipz.size(); + amrex::Array4 u_ad = makeArray4( + p, + amrex::Box(bx) + .grow(1, 1) + .grow(2, 1) + .surroundingNodes(0), + 1); + p += u_ad.size(); + amrex::Array4 v_ad = makeArray4( + p, + amrex::Box(bx) + .grow(0, 1) + .grow(2, 1) + .surroundingNodes(1), + 1); + p += v_ad.size(); + amrex::Array4 w_ad = makeArray4( + p, + amrex::Box(bx) + .grow(0, 1) + .grow(1, 1) + .surroundingNodes(2), + 1); + p += w_ad.size(); + + switch (godunov_scheme) { + /* + case godunov::scheme::PPM_NOLIM: { + godunov::predict_ppm( + lev, bxg1, ICNS::ndim, Imx, Ipx, Imy, Ipy, Imz, + Ipz, a_vel, a_vel, geom, dt, bcrec_device, + false); + break; + } + case godunov::scheme::WENOJS: { + godunov::predict_weno( + lev, bxg1, ICNS::ndim, Imx, Ipx, Imy, Ipy, Imz, + Ipz, a_vel, a_vel, geom, dt, bcrec_device, + true); + break; + } + */ + case godunov::scheme::WENOZ_CENTRAL: + case godunov::scheme::CENTRAL: + case godunov::scheme::WENOZ: { + godunov::calculate_weno( + lev, bxg1, ICNS::ndim, Imx, Ipx, Imy, Ipy, Imz, + Ipz, a_vel, a_vel, geom, bcrec_device, false); + break; + } + default: { + amrex::Abort( + "Only PPM_NOLIM, WENOZ, and WENOJS use this " + "code path"); + } + } + + godunov::make_trans_velocities( + lev, amrex::Box(u_ad), amrex::Box(v_ad), + amrex::Box(w_ad), u_ad, v_ad, w_ad, Imx, Ipx, Imy, + Ipy, Imz, Ipz, a_vel, a_f, geom, dt, bcrec_device, + godunov_use_forces_in_trans); + + godunov::calculate_godunov( + lev, bx, ICNS::ndim, xbx, ybx, zbx, a_umac, a_vmac, + a_wmac, a_vel, u_ad, v_ad, w_ad, Imx, Ipx, Imy, Ipy, + Imz, Ipz, a_f, p, geom, bcrec_device); + + amrex::Gpu::streamSynchronize(); // otherwise we might + // be using too much + // memory + } + } + } else if ( + (godunov_scheme == godunov::scheme::PPM) || + (godunov_scheme == godunov::scheme::PLM) || + (godunov_scheme == godunov::scheme::BDS)) { + amrex::Abort( + "AMReX-Hydro routines are not supported for non-linear " + "iterations"); + /* + const bool godunov_use_ppm = + godunov_scheme == godunov::scheme::PPM; + HydroUtils::ExtrapVelToFaces( + dof_field(lev), src_term(lev), u_mac(lev), + v_mac(lev), w_mac(lev), dof_field.bcrec(), + dof_field.bcrec_device().data(), + repo.mesh().Geom(lev), dt, godunov_use_ppm, + godunov_use_forces_in_trans, premac_advection_type); + */ + } else { + amrex::Abort("Invalid godunov scheme"); + } + } + + if (m_verbose > 2) { + diagnostics::PrintMaxMACVelLocations(repo, "before MAC projection"); + } + + // Populate boundaries (valid cells) using velocity BCs + m_macproj_op.set_inflow_velocity(time + 0.5 * dt); + + // MAC projection + m_macproj_op(fstate, dt); + + // Fill mac velocities (ghost cells) using velocity BCs + if (fvm::Godunov::nghost_state > 0) { + amrex::Array mac_vel = { + AMREX_D_DECL(&u_mac, &v_mac, &w_mac)}; + dof_field.fillpatch_sibling_fields( + time + 0.5 * dt, u_mac.num_grow(), mac_vel); + } + + for (int lev = 0; lev < repo.num_active_levels(); ++lev) { + u_mac(lev).FillBoundary(geom[lev].periodicity()); + v_mac(lev).FillBoundary(geom[lev].periodicity()); + w_mac(lev).FillBoundary(geom[lev].periodicity()); + } + + if (m_verbose > 2) { + diagnostics::PrintMaxMACVelLocations(repo, "after MAC projection"); + } + } + void operator()(const FieldState fstate, const amrex::Real dt) { const auto& repo = fields.repo; @@ -508,6 +707,12 @@ struct AdvectionOp m_macproj_op(fstate, dt); } + void preadvect_spatial( + const FieldState fstate, + const amrex::Real dt, + const amrex::Real /*time*/) + {} + void operator()(const FieldState fstate, const amrex::Real /*unused*/) { diff --git a/amr-wind/equation_systems/vof/vof_advection.H b/amr-wind/equation_systems/vof/vof_advection.H index 94eaab2af4..e51b6278c1 100644 --- a/amr-wind/equation_systems/vof/vof_advection.H +++ b/amr-wind/equation_systems/vof/vof_advection.H @@ -42,6 +42,12 @@ struct AdvectionOp const amrex::Real /*unused*/) {} + void preadvect_spatial( + const FieldState /*unused*/, + const amrex::Real /*unused*/, + const amrex::Real /*unused*/) + {} + void operator()(const FieldState /*unused*/, const amrex::Real dt) { static_assert( diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 843d91f4f8..65c8cc2edf 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -514,7 +514,7 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) } // Extrapolate and apply MAC projection for advection velocities - icns().pre_advection_actions(amr_wind::FieldState::Old); + icns().nonlinear_advection_actions(amr_wind::FieldState::Old); // For scalars only first // ************************************************************************************* @@ -581,8 +581,20 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) field.num_comp(), 1); } - // With scalars computed, compute advection of momentum - icns().compute_advection_term(amr_wind::FieldState::Old); + // Interpolate the NPH velocity from Old and New states + // Fillpatch the velocity + // change to new_time or remove + // icns().fields().field.fillpatch(0.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 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 @@ -596,19 +608,6 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) // ************************************************************************************* icns().compute_predictor_rhs(m_diff_type); - // Interpolate the NPH velocity from Old and New states - // Fillpatch the velocity - icns().fields().field.fillpatch(0.0); - // 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()); - - // Recompute advection of momentum - icns().compute_advection_term(amr_wind::FieldState::NPH); - // Recompute source and diffusion terms (if requested?) /* icns().compute_source_term(amr_wind::FieldState::NPH) icns().compute_diffusion_term(amr_wind::FieldState::NPH)*/ From de014108d0bed8b8a3f9b53589f16234d2bb7bbd Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Sat, 24 Aug 2024 16:50:08 -0600 Subject: [PATCH 10/47] Changing the non-linear iteration loop and adding amr_wind::io::print_nonlinear_residual() --- amr-wind/incflo.cpp | 7 ++-- amr-wind/incflo_advance.cpp | 17 ++++++++++ amr-wind/utilities/console_io.H | 3 ++ amr-wind/utilities/console_io.cpp | 55 +++++++++++++++++++++++++++++++ 4 files changed, 77 insertions(+), 5 deletions(-) diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index 3b1fac51d9..cdec9f4333 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -318,13 +318,10 @@ void incflo::do_advance() prescribe_advance(); } else { + advance(); for (int it_nl = 1; it_nl < m_adv_iters; ++it_nl) { - if (it_nl == 1) { - advance(); - } else { - advance_nonlinear(); - } + advance_nonlinear(); } } if (m_sim.has_overset()) { diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 65c8cc2edf..c4eaabc75a 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -515,6 +515,7 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) // Extrapolate and apply MAC projection for advection velocities icns().nonlinear_advection_actions(amr_wind::FieldState::Old); + // icns().nonlinear_advection_actions(amr_wind::FieldState::NPH); // For scalars only first // ************************************************************************************* @@ -662,6 +663,22 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) if (m_verbose > 2) { PrintMaxVelLocations("after nodal projection"); } + + // ScratchField to store the old np1 + // Should nghost be 0 here? + 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), -0.5, + 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); } // // Apply corrector: diff --git a/amr-wind/utilities/console_io.H b/amr-wind/utilities/console_io.H index 974c92b3dc..7b799ab279 100644 --- a/amr-wind/utilities/console_io.H +++ b/amr-wind/utilities/console_io.H @@ -3,6 +3,7 @@ #include #include "AMReX_MLMG.H" +#include "amr-wind/CFDSim.H" namespace amr_wind::io { @@ -20,6 +21,8 @@ void print_mlmg_info(const std::string& solve_name, const amrex::MLMG& mlmg); void print_tpls(std::ostream& /*out*/); +void print_nonlinear_residual(CFDSim& sim); + } // namespace amr_wind::io #endif /* CONSOLE_IO_H */ diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index bff2a9a267..49b6b38859 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -4,6 +4,7 @@ #include "amr-wind/AMRWindVersion.H" #include "AMReX.H" #include "AMReX_OpenMP.H" +#include "amr-wind/CFDSim.H" #ifdef AMR_WIND_USE_NETCDF #include "netcdf.h" @@ -214,4 +215,58 @@ void print_tpls(std::ostream& out) } } +void print_nonlinear_residual(CFDSim& sim) +{ + // using namespace amrex; + const int nlevels = sim.repo().num_active_levels(); + const auto& geom = sim.mesh().Geom(); + const auto& mesh = sim.mesh(); + + const auto& velocity_new = sim.pde_manager().icns().fields().field; + const auto& vel_np1_old = sim.repo().get_field("vel_np1_old"); + auto& vel_diff = sim.repo().get_field("vel_diff"); + + // amrex::Real total_volume_frac = 0.0; + amrex::Real rms_ucell = 0.0; + amrex::Real rms_vcell = 0.0; + amrex::Real rms_wcell = 0.0; + + for (int lev = 0; lev < nlevels; ++lev) { + + amrex::iMultiFab level_mask; + if (lev < nlevels - 1) { + level_mask = makeFineMask( + mesh.boxArray(lev), mesh.DistributionMap(lev), + mesh.boxArray(lev + 1), mesh.refRatio(lev), 1, 0); + } else { + level_mask.define( + mesh.boxArray(lev), mesh.DistributionMap(lev), 1, 0, + amrex::MFInfo()); + level_mask.setVal(1); + } + + const auto& vnew = velocity_new(lev); + const auto& velstar = vel_np1_old(lev); + auto& veldiff = vel_diff(lev); + + for (amrex::MFIter mfi(velocity_new(lev)); mfi.isValid(); ++mfi) { + const auto bx = mfi.tilebox(); + const auto& velstar_arr = velstar.const_array(mfi); + const auto& veldiff_arr = veldiff.array(mfi); + const auto& velnew_arr = vnew.const_array(mfi); + + amrex::ParallelFor( + bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + veldiff_arr(i, j, k, 0) = + velnew_arr(i, j, k, 0) - velstar_arr(i, j, k, 0); + veldiff_arr(i, j, k, 1) = + velnew_arr(i, j, k, 1) - velstar_arr(i, j, k, 1); + veldiff_arr(i, j, k, 2) = + velnew_arr(i, j, k, 2) - velstar_arr(i, j, k, 2); + }); + } + } + amrex::ParallelDescriptor::ReduceRealSum(rms_ucell); +} + } // namespace amr_wind::io From bdd1004117b6119718af89aa61a41b2e335b42fd Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Sat, 24 Aug 2024 19:56:48 -0600 Subject: [PATCH 11/47] Passing ScratchField instances as arguments for non-linear residual calculation --- amr-wind/incflo_advance.cpp | 2 +- amr-wind/utilities/console_io.H | 4 +++- amr-wind/utilities/console_io.cpp | 33 +++++++++++++++++++++---------- 3 files changed, 27 insertions(+), 12 deletions(-) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index c4eaabc75a..fe1bee9ba7 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -678,7 +678,7 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) icns().fields().field.state(amr_wind::FieldState::NPH), 0, 0, icns().fields().field.num_comp(), 1); - amr_wind::io::print_nonlinear_residual(m_sim); + amr_wind::io::print_nonlinear_residual(m_sim, *vel_diff, *vel_np1_old); } // // Apply corrector: diff --git a/amr-wind/utilities/console_io.H b/amr-wind/utilities/console_io.H index 7b799ab279..223ed95b12 100644 --- a/amr-wind/utilities/console_io.H +++ b/amr-wind/utilities/console_io.H @@ -21,7 +21,9 @@ void print_mlmg_info(const std::string& solve_name, const amrex::MLMG& mlmg); void print_tpls(std::ostream& /*out*/); -void print_nonlinear_residual(CFDSim& sim); +// void print_nonlinear_residual(CFDSim& sim); +void print_nonlinear_residual( + CFDSim& sim, ScratchField& vel_diff, ScratchField& vel_star); } // namespace amr_wind::io diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index 49b6b38859..bdce306af3 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -215,7 +215,8 @@ void print_tpls(std::ostream& out) } } -void print_nonlinear_residual(CFDSim& sim) +void print_nonlinear_residual( + CFDSim& sim, ScratchField& vel_diff, ScratchField& vel_star) { // using namespace amrex; const int nlevels = sim.repo().num_active_levels(); @@ -223,13 +224,8 @@ void print_nonlinear_residual(CFDSim& sim) const auto& mesh = sim.mesh(); const auto& velocity_new = sim.pde_manager().icns().fields().field; - const auto& vel_np1_old = sim.repo().get_field("vel_np1_old"); - auto& vel_diff = sim.repo().get_field("vel_diff"); - - // amrex::Real total_volume_frac = 0.0; - amrex::Real rms_ucell = 0.0; - amrex::Real rms_vcell = 0.0; - amrex::Real rms_wcell = 0.0; + // const auto& vel_np1_old = sim.repo().get_field("vel_np1_old"); + // auto& vel_diff = sim.repo().get_field("vel_diff"); for (int lev = 0; lev < nlevels; ++lev) { @@ -246,7 +242,9 @@ void print_nonlinear_residual(CFDSim& sim) } const auto& vnew = velocity_new(lev); - const auto& velstar = vel_np1_old(lev); + // const auto& velstar = vel_np1_old(lev); + // auto& veldiff = vel_diff(lev); + const auto& velstar = vel_star(lev); auto& veldiff = vel_diff(lev); for (amrex::MFIter mfi(velocity_new(lev)); mfi.isValid(); ++mfi) { @@ -266,7 +264,22 @@ void print_nonlinear_residual(CFDSim& sim) }); } } - amrex::ParallelDescriptor::ReduceRealSum(rms_ucell); + + amrex::Real rms_ucell = 0.0; + amrex::Real rms_vcell = 0.0; + amrex::Real rms_wcell = 0.0; + + for (int lev = 0; lev < nlevels; ++lev) { + rms_ucell += vel_diff(lev).norm2(0); + rms_vcell += vel_diff(lev).norm2(1); + rms_wcell += vel_diff(lev).norm2(2); + amrex::Print() << "Non-linear residual for u velocity in level" << lev + << "is " << rms_ucell << std::endl; + amrex::Print() << "Non-linear residual for v velocity in level" << lev + << "is " << rms_vcell << std::endl; + amrex::Print() << "Non-linear residual for w velocity in level" << lev + << "is " << rms_wcell << std::endl; + } } } // namespace amr_wind::io From bdb89dd16ce728b98d75cd68172367e20848ba6c Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Mon, 26 Aug 2024 13:48:34 -0600 Subject: [PATCH 12/47] Fixed print_nonlinear_residual() and moved the lincomb operation for NPH velocity before the scalar equations --- amr-wind/incflo.cpp | 3 +- amr-wind/incflo_advance.cpp | 32 +++++--- amr-wind/utilities/console_io.cpp | 22 +++--- .../abl_godunov_wenoz_nl.inp | 74 +++++++++++++++++++ 4 files changed, 112 insertions(+), 19 deletions(-) create mode 100644 test/test_files/abl_godunov_wenoz_nl/abl_godunov_wenoz_nl.inp diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index cdec9f4333..2f1519f1de 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -320,7 +320,8 @@ void incflo::do_advance() advance(); for (int it_nl = 1; it_nl < m_adv_iters; ++it_nl) { - + amrex::Print() << "Iteration " << it_nl << " in NL loop" + << std::endl; advance_nonlinear(); } } diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index fe1bee9ba7..0ffa7c2a09 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -436,6 +436,14 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) 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 // ************************************************************************************* @@ -514,8 +522,9 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) } // Extrapolate and apply MAC projection for advection velocities - icns().nonlinear_advection_actions(amr_wind::FieldState::Old); + // icns().nonlinear_advection_actions(amr_wind::FieldState::Old); // icns().nonlinear_advection_actions(amr_wind::FieldState::NPH); + icns().nonlinear_advection_actions(amr_wind::FieldState::New); // For scalars only first // ************************************************************************************* @@ -586,13 +595,16 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) // Fillpatch the velocity // change to new_time or remove // icns().fields().field.fillpatch(0.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()); + + /* +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 advection of momentum at NPH icns().compute_advection_term(amr_wind::FieldState::NPH); @@ -673,12 +685,14 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) "vel_diff", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL); // lincomb to get old np1 amr_wind::field_ops::lincomb( - (*vel_np1_old), -0.5, + (*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/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index bdce306af3..65f9b74699 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -252,15 +252,19 @@ void print_nonlinear_residual( const auto& velstar_arr = velstar.const_array(mfi); const auto& veldiff_arr = veldiff.array(mfi); const auto& velnew_arr = vnew.const_array(mfi); + const auto& levelmask_arr = level_mask.const_array(mfi); amrex::ParallelFor( bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { veldiff_arr(i, j, k, 0) = - velnew_arr(i, j, k, 0) - velstar_arr(i, j, k, 0); + (velnew_arr(i, j, k, 0) - velstar_arr(i, j, k, 0)) * + levelmask_arr(i, j, k); veldiff_arr(i, j, k, 1) = - velnew_arr(i, j, k, 1) - velstar_arr(i, j, k, 1); + (velnew_arr(i, j, k, 1) - velstar_arr(i, j, k, 1)) * + levelmask_arr(i, j, k); veldiff_arr(i, j, k, 2) = - velnew_arr(i, j, k, 2) - velstar_arr(i, j, k, 2); + (velnew_arr(i, j, k, 2) - velstar_arr(i, j, k, 2)) * + levelmask_arr(i, j, k); }); } } @@ -273,13 +277,13 @@ void print_nonlinear_residual( rms_ucell += vel_diff(lev).norm2(0); rms_vcell += vel_diff(lev).norm2(1); rms_wcell += vel_diff(lev).norm2(2); - amrex::Print() << "Non-linear residual for u velocity in level" << lev - << "is " << rms_ucell << std::endl; - amrex::Print() << "Non-linear residual for v velocity in level" << lev - << "is " << rms_vcell << std::endl; - amrex::Print() << "Non-linear residual for w velocity in level" << lev - << "is " << rms_wcell << std::endl; } + amrex::Print() << "Non-linear residual for u velocity " << rms_ucell + << std::endl; + amrex::Print() << "Non-linear residual for v velocity " << rms_vcell + << std::endl; + amrex::Print() << "Non-linear residual for w velocity " << rms_wcell + << std::endl; } } // namespace amr_wind::io diff --git a/test/test_files/abl_godunov_wenoz_nl/abl_godunov_wenoz_nl.inp b/test/test_files/abl_godunov_wenoz_nl/abl_godunov_wenoz_nl.inp new file mode 100644 index 0000000000..0a3d6c92bc --- /dev/null +++ b/test/test_files/abl_godunov_wenoz_nl/abl_godunov_wenoz_nl.inp @@ -0,0 +1,74 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +#time.stop_time = 22000.0 # Max (simulated) time to evolve +time.stop_time = 50.0 # Max (simulated) time to evolve +time.max_step = 100 # Max number of time steps + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.5 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +time.plot_interval = 10 # Steps between plot files +time.checkpoint_interval = -5 # Steps between checkpoint files + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0. 0. -9.81 # Gravitational force (3D) +incflo.density = 1.0 # Reference density +incflo.use_godunov = 1 +incflo.godunov_type = "weno_z" +incflo.diffusion_type = 2 +incflo.advection_iterations = 10 +transport.viscosity = 1.0e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.3333 +turbulence.model = Smagorinsky +Smagorinsky_coeffs.Cs = 0.135 + + +incflo.physics = ABL +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing ABLForcing +BoussinesqBuoyancy.reference_temperature = 300.0 +ABL.reference_temperature = 300.0 +CoriolisForcing.latitude = 41.3 +ABLForcing.abl_forcing_height = 90 + +incflo.velocity = 6.128355544951824 5.142300877492314 0.0 + +ABL.temperature_heights = 650.0 750.0 1000.0 +ABL.temperature_values = 300.0 308.0 308.75 + +ABL.kappa = .41 +ABL.surface_roughness_z0 = 0.15 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 48 48 48 # Grid cells at coarsest AMRlevel +amr.max_level = 0 # Max AMR level in hierarchy + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 1000. 1000. 1000. # Hi corner coordinates +geometry.is_periodic = 1 1 0 # Periodicity x y z (0/1) + +# Boundary conditions +zlo.type = "wall_model" + +zhi.type = "slip_wall" +zhi.temperature_type = "fixed_gradient" +zhi.temperature = 0.003 # tracer is used to specify potential temperature gradient + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# VERBOSITY # +#.......................................# +incflo.verbose = 0 # incflo_level From 53170ec02fcce3e32baf30ffde503fd0a8319f64 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Sun, 1 Sep 2024 17:52:47 -0600 Subject: [PATCH 13/47] Rebased on the main branch --- amr-wind/convection/Godunov.H | 1 - amr-wind/convection/incflo_godunov_central.H | 139 ------------ amr-wind/equation_systems/AdvOp_Godunov.H | 77 +++---- amr-wind/equation_systems/AdvOp_MOL.H | 6 - amr-wind/equation_systems/PDE.H | 8 - amr-wind/equation_systems/PDEBase.H | 3 - .../equation_systems/icns/icns_advection.H | 214 +----------------- amr-wind/equation_systems/vof/vof_advection.H | 6 - amr-wind/incflo_advance.cpp | 3 +- 9 files changed, 38 insertions(+), 419 deletions(-) delete mode 100644 amr-wind/convection/incflo_godunov_central.H diff --git a/amr-wind/convection/Godunov.H b/amr-wind/convection/Godunov.H index a93df095c4..dcbf32bdbb 100644 --- a/amr-wind/convection/Godunov.H +++ b/amr-wind/convection/Godunov.H @@ -10,7 +10,6 @@ namespace godunov { enum class scheme { PLM, PPM, PPM_NOLIM, BDS, WENO_JS, WENOZ, MINMOD, UPWIND }; - } // namespace godunov #endif /* Godunov_H */ diff --git a/amr-wind/convection/incflo_godunov_central.H b/amr-wind/convection/incflo_godunov_central.H deleted file mode 100644 index 58d1e96699..0000000000 --- a/amr-wind/convection/incflo_godunov_central.H +++ /dev/null @@ -1,139 +0,0 @@ -#ifndef GODUNOV_CENTRAL_H -#define GODUNOV_CENTRAL_H - -#include -#include - -/* This header file contains the inlined __host__ __device__ functions required - for the scalar advection routines for 3D Godunov. It also contains function - declarations for controlling host functions. */ - -namespace { - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void central( - const amrex::Real sm1, - const amrex::Real s0, - const amrex::Real sp1, - amrex::Real& dsm, - amrex::Real& dsp) -{ - // Calculate gradients on both sides - dsp = sp1 - s0; - dsm = s0 - sm1; -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_central_bc( - const int n, - const amrex::Real sm1, - const amrex::Real s0, - const amrex::Real sp1, - amrex::Real& dsm, - amrex::Real& dsp, - const int bclo, - const int bchi, - const int domlo, - const int domhi) -{ - using namespace amrex; - - if (bclo == BCType::ext_dir || bclo == BCType::hoextrap) { - if (n == domlo) { - // Ensure that left-side slope is used unchanged - dsm = s0 - sm1; - } - } - - if (bchi == BCType::ext_dir || bchi == BCType::hoextrap) { - if (n == domhi) { - // Ensure that the right-side slope is used unchanged - dsp = sp1 - s0; - } - } -} - -// !!! No upwinding for central scheme !!! -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_central_x( - const int i, - const int j, - const int k, - const int n, - const amrex::Real dx, - amrex::Real& Im, - amrex::Real& Ip, - const amrex::Array4& S, - const amrex::BCRec& bc, - const int domlo, - const int domhi) -{ - using namespace amrex; - amrex::Real sm1 = S(i - 1, j, k, n); - amrex::Real s0 = S(i, j, k, n); - amrex::Real sp1 = S(i + 1, j, k, n); - amrex::Real dsp = 0.0; - amrex::Real dsm = 0.0; - central(sm1, s0, sp1, dsm, dsp); - Godunov_central_bc( - i, sm1, s0, sp1, dsm, dsp, bc.lo(0), bc.hi(0), domlo, domhi); - - // Interpolate to edges - Ip = s0 + 0.5 * dsp; - Im = s0 - 0.5 * dsm; -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_central_y( - const int i, - const int j, - const int k, - const int n, - const amrex::Real dx, - amrex::Real& Im, - amrex::Real& Ip, - const amrex::Array4& S, - const amrex::BCRec& bc, - const int domlo, - const int domhi) -{ - using namespace amrex; - amrex::Real sm1 = S(i, j - 1, k, n); - amrex::Real s0 = S(i, j, k, n); - amrex::Real sp1 = S(i, j + 1, k, n); - amrex::Real dsp = 0.0; - amrex::Real dsm = 0.0; - central(sm1, s0, sp1, dsm, dsp); - Godunov_central_bc( - j, sm1, s0, sp1, dsm, dsp, bc.lo(1), bc.hi(1), domlo, domhi); - - Ip = s0 + 0.5 * dsp; - Im = s0 - 0.5 * dsm; -} - -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Godunov_central_z( - const int i, - const int j, - const int k, - const int n, - const amrex::Real dx, - amrex::Real& Im, - amrex::Real& Ip, - const amrex::Array4& S, - const amrex::BCRec& bc, - const int domlo, - const int domhi) -{ - using namespace amrex; - amrex::Real sm1 = S(i, j, k - 1, n); - amrex::Real s0 = S(i, j, k, n); - amrex::Real sp1 = S(i, j, k + 1, n); - amrex::Real dsp = 0.0; - amrex::Real dsm = 0.0; - central(sm1, s0, sp1, dsm, dsp); - Godunov_central_bc( - k, sm1, s0, sp1, dsm, dsp, bc.lo(2), bc.hi(2), domlo, domhi); - - Ip = s0 + 0.5 * dsp; - Im = s0 - 0.5 * dsm; -} - -} // namespace - -#endif diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 64b940aa27..6c6b07b9d4 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -80,12 +80,6 @@ struct AdvectionOp< const amrex::Real /*unused*/) {} - void preadvect_spatial( - const FieldState /*unused*/, - const amrex::Real /*unused*/, - const amrex::Real /*unused*/) - {} - void operator()(const FieldState fstate, const amrex::Real dt) { static_assert( @@ -145,7 +139,6 @@ struct AdvectionOp< rho_arr(i, j, k) * tra_arr(i, j, k, n); }); } -<<<<<<< HEAD if ((godunov_scheme == godunov::scheme::PPM) || (godunov_scheme == godunov::scheme::PPM_NOLIM) || @@ -168,7 +161,6 @@ struct AdvectionOp< limiter_type = PPM::WENOZ; } else if (godunov_scheme == godunov::scheme::WENO_JS) { limiter_type = PPM::WENO_JS; -======= } else { limiter_type = PPM::default_limiter; } @@ -196,50 +188,51 @@ struct AdvectionOp< } else { amrex::Abort("Invalid godunov scheme"); } + amrex::Gpu::streamSynchronize(); } - amrex::Gpu::streamSynchronize(); } - } - amrex::Vector> - fluxes(repo.num_active_levels()); - for (int lev = 0; lev < repo.num_active_levels(); ++lev) { - fluxes[lev][0] = &(*flux_x)(lev); - fluxes[lev][1] = &(*flux_y)(lev); - fluxes[lev][2] = &(*flux_z)(lev); - } + amrex::Vector> fluxes( + repo.num_active_levels()); + for (int lev = 0; lev < repo.num_active_levels(); ++lev) { + fluxes[lev][0] = &(*flux_x)(lev); + fluxes[lev][1] = &(*flux_y)(lev); + fluxes[lev][2] = &(*flux_z)(lev); + } - // In order to enforce conservation across coarse-fine boundaries we - // must be sure to average down the fluxes before we use them - for (int lev = repo.num_active_levels() - 1; lev > 0; --lev) { - amrex::IntVect rr = - geom[lev].Domain().size() / geom[lev - 1].Domain().size(); - amrex::average_down_faces( - GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1], rr, geom[lev - 1]); - } + // In order to enforce conservation across coarse-fine boundaries we + // must be sure to average down the fluxes before we use them + for (int lev = repo.num_active_levels() - 1; lev > 0; --lev) { + amrex::IntVect rr = + geom[lev].Domain().size() / geom[lev - 1].Domain().size(); + amrex::average_down_faces( + GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1], rr, + geom[lev - 1]); + } - for (int lev = 0; lev < repo.num_active_levels(); ++lev) { + for (int lev = 0; lev < repo.num_active_levels(); ++lev) { #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for (amrex::MFIter mfi(dof_field(lev), amrex::TilingIfNotGPU()); - mfi.isValid(); ++mfi) { - const auto& bx = mfi.tilebox(); - - HydroUtils::ComputeDivergence( - bx, conv_term(lev).array(mfi), (*flux_x)(lev).array(mfi), - (*flux_y)(lev).array(mfi), (*flux_z)(lev).array(mfi), PDE::ndim, - geom[lev], amrex::Real(-1.0), fluxes_are_area_weighted); + for (amrex::MFIter mfi(dof_field(lev), amrex::TilingIfNotGPU()); + mfi.isValid(); ++mfi) { + const auto& bx = mfi.tilebox(); + + HydroUtils::ComputeDivergence( + bx, conv_term(lev).array(mfi), (*flux_x)(lev).array(mfi), + (*flux_y)(lev).array(mfi), (*flux_z)(lev).array(mfi), + PDE::ndim, geom[lev], amrex::Real(-1.0), + fluxes_are_area_weighted); + } } } -} - -PDEFields& fields; -Field & density; -Field & u_mac; -Field & v_mac; -Field & w_mac; -amrex::Gpu::DeviceVector iconserv; + + PDEFields& fields; + Field& density; + Field& u_mac; + Field& v_mac; + Field& w_mac; + amrex::Gpu::DeviceVector iconserv; godunov::scheme godunov_scheme = godunov::scheme::WENOZ; std::string godunov_type{"weno_z"}; diff --git a/amr-wind/equation_systems/AdvOp_MOL.H b/amr-wind/equation_systems/AdvOp_MOL.H index 97fc278ab1..b276e39d78 100644 --- a/amr-wind/equation_systems/AdvOp_MOL.H +++ b/amr-wind/equation_systems/AdvOp_MOL.H @@ -43,12 +43,6 @@ struct AdvectionOp< const amrex::Real /*unused*/) {} - void preadvect_spatial( - const FieldState /*unused*/, - const amrex::Real /*unused*/, - const amrex::Real /*unused*/) - {} - void operator()(const FieldState fstate, const amrex::Real /*unused*/) { static_assert( diff --git a/amr-wind/equation_systems/PDE.H b/amr-wind/equation_systems/PDE.H index f74526c3ac..485616d680 100644 --- a/amr-wind/equation_systems/PDE.H +++ b/amr-wind/equation_systems/PDE.H @@ -140,14 +140,6 @@ public: m_adv_op->preadvect(fstate, m_time.delta_t(), m_time.current_time()); } - void nonlinear_advection_actions(const FieldState fstate) override - { - BL_PROFILE( - "amr-wind::" + this->identifier() + "::pre_advection_actions"); - m_adv_op->preadvect_spatial( - fstate, m_time.delta_t(), m_time.current_time()); - } - void compute_predictor_rhs(const DiffusionType difftype) override { BL_PROFILE( diff --git a/amr-wind/equation_systems/PDEBase.H b/amr-wind/equation_systems/PDEBase.H index 19206026c9..5abbd19278 100644 --- a/amr-wind/equation_systems/PDEBase.H +++ b/amr-wind/equation_systems/PDEBase.H @@ -74,9 +74,6 @@ public: //! Perform necessary steps to prepare for advection calculations virtual void pre_advection_actions(const FieldState fstate) = 0; - //! Perform necessary steps to prepare for nonlinear advection calculations - virtual void nonlinear_advection_actions(const FieldState fstate) = 0; - //! Compute the time derivative and advective term for the PDE system virtual void compute_advection_term(const FieldState fstate) = 0; diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 3b72339678..81b1512b7a 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -120,10 +120,6 @@ struct AdvectionOp godunov_scheme = godunov::scheme::WENO_JS; } else if (amrex::toLower(godunov_type) == "weno_z") { godunov_scheme = godunov::scheme::WENOZ; - } else if (amrex::toLower(godunov_type) == "central") { - godunov_scheme = godunov::scheme::CENTRAL; - } else if (amrex::toLower(godunov_type) == "weno_z_central") { - godunov_scheme = godunov::scheme::WENOZ_CENTRAL; } else { amrex::Abort( "Invalid godunov_type specified. For godunov_type select " @@ -235,205 +231,6 @@ struct AdvectionOp } } - void preadvect_spatial( - const FieldState fstate, const amrex::Real dt, const amrex::Real time) - { - - const auto& repo = fields.repo; - const auto& geom = repo.mesh().Geom(); - - const auto& src_term = fields.src_term; - const auto& dof_field = fields.field.state(fstate); - auto bcrec_device = dof_field.bcrec_device(); - - for (int lev = 0; lev < repo.num_active_levels(); ++lev) { - u_mac(lev).setBndry(0.0); - v_mac(lev).setBndry(0.0); - w_mac(lev).setBndry(0.0); - } - - // - // Predict - // - for (int lev = 0; lev < repo.num_active_levels(); ++lev) { - if ((godunov_scheme == godunov::scheme::PPM_NOLIM) || - (godunov_scheme == godunov::scheme::WENOJS) || - (godunov_scheme == godunov::scheme::WENOZ) || - (godunov_scheme == godunov::scheme::CENTRAL) || - (godunov_scheme == godunov::scheme::WENOZ_CENTRAL)) { -#ifdef AMREX_USE_OMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - { - amrex::FArrayBox scratch; - for (amrex::MFIter mfi( - dof_field(lev), amrex::TilingIfNotGPU()); - mfi.isValid(); ++mfi) { - amrex::Box const& bx = mfi.tilebox(); - amrex::Box const& bxg1 = amrex::grow(bx, 1); - amrex::Box const& xbx = mfi.nodaltilebox(0); - amrex::Box const& ybx = mfi.nodaltilebox(1); - amrex::Box const& zbx = mfi.nodaltilebox(2); - - amrex::Array4 const& a_umac = - u_mac(lev).array(mfi); - amrex::Array4 const& a_vmac = - v_mac(lev).array(mfi); - amrex::Array4 const& a_wmac = - w_mac(lev).array(mfi); - amrex::Array4 const& a_vel = - dof_field(lev).const_array(mfi); - amrex::Array4 const& a_f = - src_term(lev).const_array(mfi); - - scratch.resize(bxg1, ICNS::ndim * 12 + 3); - // Elixir eli = scratch.elixir(); // not - // needed because of streamSynchronize - // later - amrex::Real* p = scratch.dataPtr(); - - amrex::Array4 Imx = - makeArray4(p, bxg1, ICNS::ndim); - p += Imx.size(); - amrex::Array4 Ipx = - makeArray4(p, bxg1, ICNS::ndim); - p += Ipx.size(); - amrex::Array4 Imy = - makeArray4(p, bxg1, ICNS::ndim); - p += Imy.size(); - amrex::Array4 Ipy = - makeArray4(p, bxg1, ICNS::ndim); - p += Ipy.size(); - amrex::Array4 Imz = - makeArray4(p, bxg1, ICNS::ndim); - p += Imz.size(); - amrex::Array4 Ipz = - makeArray4(p, bxg1, ICNS::ndim); - p += Ipz.size(); - amrex::Array4 u_ad = makeArray4( - p, - amrex::Box(bx) - .grow(1, 1) - .grow(2, 1) - .surroundingNodes(0), - 1); - p += u_ad.size(); - amrex::Array4 v_ad = makeArray4( - p, - amrex::Box(bx) - .grow(0, 1) - .grow(2, 1) - .surroundingNodes(1), - 1); - p += v_ad.size(); - amrex::Array4 w_ad = makeArray4( - p, - amrex::Box(bx) - .grow(0, 1) - .grow(1, 1) - .surroundingNodes(2), - 1); - p += w_ad.size(); - - switch (godunov_scheme) { - /* - case godunov::scheme::PPM_NOLIM: { - godunov::predict_ppm( - lev, bxg1, ICNS::ndim, Imx, Ipx, Imy, Ipy, Imz, - Ipz, a_vel, a_vel, geom, dt, bcrec_device, - false); - break; - } - case godunov::scheme::WENOJS: { - godunov::predict_weno( - lev, bxg1, ICNS::ndim, Imx, Ipx, Imy, Ipy, Imz, - Ipz, a_vel, a_vel, geom, dt, bcrec_device, - true); - break; - } - */ - case godunov::scheme::WENOZ_CENTRAL: - case godunov::scheme::CENTRAL: - case godunov::scheme::WENOZ: { - godunov::calculate_weno( - lev, bxg1, ICNS::ndim, Imx, Ipx, Imy, Ipy, Imz, - Ipz, a_vel, a_vel, geom, bcrec_device, false); - break; - } - default: { - amrex::Abort( - "Only PPM_NOLIM, WENOZ, and WENOJS use this " - "code path"); - } - } - - godunov::make_trans_velocities( - lev, amrex::Box(u_ad), amrex::Box(v_ad), - amrex::Box(w_ad), u_ad, v_ad, w_ad, Imx, Ipx, Imy, - Ipy, Imz, Ipz, a_vel, a_f, geom, dt, bcrec_device, - godunov_use_forces_in_trans); - - godunov::calculate_godunov( - lev, bx, ICNS::ndim, xbx, ybx, zbx, a_umac, a_vmac, - a_wmac, a_vel, u_ad, v_ad, w_ad, Imx, Ipx, Imy, Ipy, - Imz, Ipz, a_f, p, geom, bcrec_device); - - amrex::Gpu::streamSynchronize(); // otherwise we might - // be using too much - // memory - } - } - } else if ( - (godunov_scheme == godunov::scheme::PPM) || - (godunov_scheme == godunov::scheme::PLM) || - (godunov_scheme == godunov::scheme::BDS)) { - amrex::Abort( - "AMReX-Hydro routines are not supported for non-linear " - "iterations"); - /* - const bool godunov_use_ppm = - godunov_scheme == godunov::scheme::PPM; - HydroUtils::ExtrapVelToFaces( - dof_field(lev), src_term(lev), u_mac(lev), - v_mac(lev), w_mac(lev), dof_field.bcrec(), - dof_field.bcrec_device().data(), - repo.mesh().Geom(lev), dt, godunov_use_ppm, - godunov_use_forces_in_trans, premac_advection_type); - */ - } else { - amrex::Abort("Invalid godunov scheme"); - } - } - - if (m_verbose > 2) { - diagnostics::PrintMaxMACVelLocations(repo, "before MAC projection"); - } - - // Populate boundaries (valid cells) using velocity BCs - m_macproj_op.set_inflow_velocity(time + 0.5 * dt); - - // MAC projection - m_macproj_op(fstate, dt); - - // Fill mac velocities (ghost cells) using velocity BCs - if (fvm::Godunov::nghost_state > 0) { - amrex::Array mac_vel = { - AMREX_D_DECL(&u_mac, &v_mac, &w_mac)}; - dof_field.fillpatch_sibling_fields( - time + 0.5 * dt, u_mac.num_grow(), mac_vel); - } - - for (int lev = 0; lev < repo.num_active_levels(); ++lev) { - u_mac(lev).FillBoundary(geom[lev].periodicity()); - v_mac(lev).FillBoundary(geom[lev].periodicity()); - w_mac(lev).FillBoundary(geom[lev].periodicity()); - } - - if (m_verbose > 2) { - diagnostics::PrintMaxMACVelLocations(repo, "after MAC projection"); - } - } - void operator()(const FieldState fstate, const amrex::Real dt) { const auto& repo = fields.repo; @@ -442,11 +239,8 @@ struct AdvectionOp const auto& src_term = fields.src_term; // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term; - // if fstate is NPH, then n and n+1 are known, need only spatial interp - bool only_spatial = - (fstate == FieldState::NPH || - godunov_scheme == godunov::scheme::CENTRAL); const auto& dof_field = fields.field.state(fstate); + bool only_spatial = (fstate == FieldState::NPH); auto flux_x = repo.create_scratch_field(ICNS::ndim, 0, amr_wind::FieldLoc::XFACE); @@ -707,12 +501,6 @@ struct AdvectionOp m_macproj_op(fstate, dt); } - void preadvect_spatial( - const FieldState fstate, - const amrex::Real dt, - const amrex::Real /*time*/) - {} - void operator()(const FieldState fstate, const amrex::Real /*unused*/) { diff --git a/amr-wind/equation_systems/vof/vof_advection.H b/amr-wind/equation_systems/vof/vof_advection.H index e51b6278c1..94eaab2af4 100644 --- a/amr-wind/equation_systems/vof/vof_advection.H +++ b/amr-wind/equation_systems/vof/vof_advection.H @@ -42,12 +42,6 @@ struct AdvectionOp const amrex::Real /*unused*/) {} - void preadvect_spatial( - const FieldState /*unused*/, - const amrex::Real /*unused*/, - const amrex::Real /*unused*/) - {} - void operator()(const FieldState /*unused*/, const amrex::Real dt) { static_assert( diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 0ffa7c2a09..f47c223b79 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -524,7 +524,8 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) // Extrapolate and apply MAC projection for advection velocities // icns().nonlinear_advection_actions(amr_wind::FieldState::Old); // icns().nonlinear_advection_actions(amr_wind::FieldState::NPH); - icns().nonlinear_advection_actions(amr_wind::FieldState::New); + //icns().nonlinear_advection_actions(amr_wind::FieldState::New); + icns().pre_advection_actions(amr_wind::FieldState::NPH); // For scalars only first // ************************************************************************************* From 0a7f355d5fd06c73cea8a242a2877ef2fb022ffc Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Mon, 2 Sep 2024 14:45:54 -0600 Subject: [PATCH 14/47] Replaced only_spatial statements with fstate checks for dt_arg --- amr-wind/equation_systems/AdvOp_Godunov.H | 4 +++- amr-wind/equation_systems/icns/icns_advection.H | 12 +++++++++--- amr-wind/incflo_advance.cpp | 3 --- 3 files changed, 12 insertions(+), 7 deletions(-) diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 6c6b07b9d4..3b5ade301a 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -168,7 +168,9 @@ struct AdvectionOp< (godunov_scheme == godunov::scheme::WENO_JS)) { m_allow_inflow_on_outflow = true; } - amrex::Real dt_arg = (only_spatial) ? 0.0 : dt; + // if state is NPH, then n and n+1 are known, and only + // spatial extrapolation is performed + amrex::Real dt_arg = (fstate == FieldState::NPH) ? 0.0 : dt; HydroUtils::ComputeFluxesOnBoxFromState( bx, PDE::ndim, mfi, diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 81b1512b7a..8bb7b1ebfa 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -190,11 +190,14 @@ struct AdvectionOp m_allow_inflow_on_outflow = true; } + // if state is NPH, then n and n+1 are known, and only + // spatial extrapolation is performed + amrex::Real dt_arg = (fstate == FieldState::NPH) ? 0.0 : dt; for (int lev = 0; lev < repo.num_active_levels(); ++lev) { HydroUtils::ExtrapVelToFaces( dof_field(lev), src_term(lev), u_mac(lev), v_mac(lev), w_mac(lev), dof_field.bcrec(), bcrec_device.data(), - repo.mesh().Geom(lev), dt, godunov_use_ppm, + repo.mesh().Geom(lev), dt_arg, godunov_use_ppm, godunov_use_forces_in_trans, premac_advection_type, limiter_type, m_allow_inflow_on_outflow); } @@ -240,7 +243,6 @@ struct AdvectionOp // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term; const auto& dof_field = fields.field.state(fstate); - bool only_spatial = (fstate == FieldState::NPH); auto flux_x = repo.create_scratch_field(ICNS::ndim, 0, amr_wind::FieldLoc::XFACE); @@ -333,6 +335,10 @@ struct AdvectionOp #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif + + // if state is NPH, then n and n+1 are known, and only + // spatial extrapolation is performed + amrex::Real dt_arg = (fstate == FieldState::NPH) ? 0.0 : dt; for (amrex::MFIter mfi(dof_field(lev), mfi_info); mfi.isValid(); ++mfi) { const auto& bx = mfi.tilebox(); @@ -347,7 +353,7 @@ struct AdvectionOp known_edge_state, u_mac(lev).const_array(mfi), v_mac(lev).const_array(mfi), w_mac(lev).const_array(mfi), divu, fq.const_array(mfi), - geom[lev], dt, dof_field.bcrec(), + geom[lev], dt_arg, dof_field.bcrec(), dof_field.bcrec_device().data(), iconserv.data(), godunov_use_ppm, godunov_use_forces_in_trans, is_velocity, fluxes_are_area_weighted, diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index f47c223b79..a8a5fc033e 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -522,9 +522,6 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) } // Extrapolate and apply MAC projection for advection velocities - // icns().nonlinear_advection_actions(amr_wind::FieldState::Old); - // icns().nonlinear_advection_actions(amr_wind::FieldState::NPH); - //icns().nonlinear_advection_actions(amr_wind::FieldState::New); icns().pre_advection_actions(amr_wind::FieldState::NPH); // For scalars only first From f839ff2970a3a72bd431e2e13db9f4111f66c8eb Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Mon, 9 Sep 2024 19:50:10 -0600 Subject: [PATCH 15/47] Moved the advection iteration loop outside of do_advance() to make it compatible with exawind-driver --- amr-wind/equation_systems/AdvOp_Godunov.H | 4 ++-- .../equation_systems/icns/icns_advection.H | 8 +++---- amr-wind/incflo.H | 2 +- amr-wind/incflo.cpp | 23 +++++++++++-------- 4 files changed, 20 insertions(+), 17 deletions(-) diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 3b5ade301a..6126c750e1 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -170,7 +170,7 @@ struct AdvectionOp< } // if state is NPH, then n and n+1 are known, and only // spatial extrapolation is performed - amrex::Real dt_arg = (fstate == FieldState::NPH) ? 0.0 : dt; + amrex::Real dt_extrap = (fstate == FieldState::NPH) ? 0.0 : dt; HydroUtils::ComputeFluxesOnBoxFromState( bx, PDE::ndim, mfi, @@ -181,7 +181,7 @@ struct AdvectionOp< known_edge_state, u_mac(lev).const_array(mfi), v_mac(lev).const_array(mfi), w_mac(lev).const_array(mfi), divu, - src_term(lev).const_array(mfi), geom[lev], dt_arg, + src_term(lev).const_array(mfi), geom[lev], dt_extrap, dof_field.bcrec(), dof_field.bcrec_device().data(), iconserv.data(), godunov_use_ppm, godunov_use_forces_in_trans, is_velocity, diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 8bb7b1ebfa..5a128cdbda 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -192,12 +192,12 @@ struct AdvectionOp // if state is NPH, then n and n+1 are known, and only // spatial extrapolation is performed - amrex::Real dt_arg = (fstate == FieldState::NPH) ? 0.0 : dt; + amrex::Real dt_extrap = (fstate == FieldState::NPH) ? 0.0 : dt; for (int lev = 0; lev < repo.num_active_levels(); ++lev) { HydroUtils::ExtrapVelToFaces( dof_field(lev), src_term(lev), u_mac(lev), v_mac(lev), w_mac(lev), dof_field.bcrec(), bcrec_device.data(), - repo.mesh().Geom(lev), dt_arg, godunov_use_ppm, + repo.mesh().Geom(lev), dt_extrap, godunov_use_ppm, godunov_use_forces_in_trans, premac_advection_type, limiter_type, m_allow_inflow_on_outflow); } @@ -338,7 +338,7 @@ struct AdvectionOp // if state is NPH, then n and n+1 are known, and only // spatial extrapolation is performed - amrex::Real dt_arg = (fstate == FieldState::NPH) ? 0.0 : dt; + amrex::Real dt_extrap = (fstate == FieldState::NPH) ? 0.0 : dt; for (amrex::MFIter mfi(dof_field(lev), mfi_info); mfi.isValid(); ++mfi) { const auto& bx = mfi.tilebox(); @@ -353,7 +353,7 @@ struct AdvectionOp known_edge_state, u_mac(lev).const_array(mfi), v_mac(lev).const_array(mfi), w_mac(lev).const_array(mfi), divu, fq.const_array(mfi), - geom[lev], dt_arg, dof_field.bcrec(), + geom[lev], dt_extrap, dof_field.bcrec(), dof_field.bcrec_device().data(), iconserv.data(), godunov_use_ppm, godunov_use_forces_in_trans, is_velocity, fluxes_are_area_weighted, diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 493ee8dd2e..0501197263 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -96,7 +96,7 @@ public: bool regrid_and_update(); void pre_advance_stage1(); void pre_advance_stage2(); - void do_advance(); + void do_advance(int inonlin); void advance(); void advance_nonlinear(); void prescribe_advance(); diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index 2f1519f1de..a5c30cc520 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -276,7 +276,8 @@ void incflo::Evolve() amrex::Real time1 = amrex::ParallelDescriptor::second(); // Advance to time t + dt - do_advance(); + for (int inonlin = 1; inonlin <= m_adv_iters; ++inonlin) + do_advance(inonlin); amrex::Print() << std::endl; amrex::Real time2 = amrex::ParallelDescriptor::second(); @@ -309,19 +310,19 @@ void incflo::Evolve() } } -void incflo::do_advance() +void incflo::do_advance(int inonlin) { - if (m_sim.has_overset()) { + if (m_sim.has_overset() && inonlin == 1) { m_ovst_ops.pre_advance_work(); } - if (m_prescribe_vel) { + if (m_prescribe_vel && inonlin == 1) { prescribe_advance(); } else { - - advance(); - for (int it_nl = 1; it_nl < m_adv_iters; ++it_nl) { - amrex::Print() << "Iteration " << it_nl << " in NL loop" - << std::endl; + if (inonlin == 1) { + advance(); + } else { + amrex::Print() << "Iteration " << inonlin + << " in advection iteration loop" << std::endl; advance_nonlinear(); } } @@ -381,7 +382,9 @@ void incflo::init_physics_and_pde() // Get number of advection iterations pp.query("advection_iterations", m_adv_iters); - m_adv_iters = std::max(1, m_adv_iters); + if (m_adv_iters < 1) + amrex::Abort( + "The number of advection iterations cannot be less than 1"); } auto& pde_mgr = m_sim.pde_manager(); From cd72385cb4a395980241464a271a0440141a9d24 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Wed, 11 Sep 2024 13:25:13 -0600 Subject: [PATCH 16/47] Added index iterations over velocity components into ParallelFor, and replaced the indivudal change-in-velocity norms with an Array instance --- amr-wind/utilities/console_io.cpp | 38 +++++++++++-------------------- 1 file changed, 13 insertions(+), 25 deletions(-) diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index 65f9b74699..32e200ba5e 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -224,8 +224,7 @@ void print_nonlinear_residual( const auto& mesh = sim.mesh(); const auto& velocity_new = sim.pde_manager().icns().fields().field; - // const auto& vel_np1_old = sim.repo().get_field("vel_np1_old"); - // auto& vel_diff = sim.repo().get_field("vel_diff"); + const int ncomp = AMREX_SPACEDIM; for (int lev = 0; lev < nlevels; ++lev) { @@ -242,8 +241,6 @@ void print_nonlinear_residual( } const auto& vnew = velocity_new(lev); - // const auto& velstar = vel_np1_old(lev); - // auto& veldiff = vel_diff(lev); const auto& velstar = vel_star(lev); auto& veldiff = vel_diff(lev); @@ -255,35 +252,26 @@ void print_nonlinear_residual( const auto& levelmask_arr = level_mask.const_array(mfi); amrex::ParallelFor( - bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - veldiff_arr(i, j, k, 0) = - (velnew_arr(i, j, k, 0) - velstar_arr(i, j, k, 0)) * - levelmask_arr(i, j, k); - veldiff_arr(i, j, k, 1) = - (velnew_arr(i, j, k, 1) - velstar_arr(i, j, k, 1)) * - levelmask_arr(i, j, k); - veldiff_arr(i, j, k, 2) = - (velnew_arr(i, j, k, 2) - velstar_arr(i, j, k, 2)) * + bx, AMREX_SPACEDIM, + [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { + veldiff_arr(i, j, k, n) = + (velnew_arr(i, j, k, n) - velstar_arr(i, j, k, n)) * levelmask_arr(i, j, k); }); } } - amrex::Real rms_ucell = 0.0; - amrex::Real rms_vcell = 0.0; - amrex::Real rms_wcell = 0.0; + amrex::Array rms_vel; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) rms_vel[idim] = 0.; for (int lev = 0; lev < nlevels; ++lev) { - rms_ucell += vel_diff(lev).norm2(0); - rms_vcell += vel_diff(lev).norm2(1); - rms_wcell += vel_diff(lev).norm2(2); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + rms_vel[idim] += vel_diff(lev).norm2(idim); + } } - amrex::Print() << "Non-linear residual for u velocity " << rms_ucell - << std::endl; - amrex::Print() << "Non-linear residual for v velocity " << rms_vcell - << std::endl; - amrex::Print() << "Non-linear residual for w velocity " << rms_wcell - << std::endl; + + amrex::Print() << "Non-linear residual for u " << rms_vel[0] << " v " + << rms_vel[1] << " w " << rms_vel[2] << std::endl; } } // namespace amr_wind::io From 55fc15b30124885d731fab950902409e92cca7a5 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Thu, 12 Sep 2024 13:35:03 -0600 Subject: [PATCH 17/47] Changing the initial advection iteration index from 1 to 0 --- amr-wind/incflo.H | 1 + amr-wind/incflo.cpp | 9 +++++---- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 0501197263..123904193f 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -97,6 +97,7 @@ public: void pre_advance_stage1(); void pre_advance_stage2(); void do_advance(int inonlin); + void post_overset_work(); void advance(); void advance_nonlinear(); void prescribe_advance(); diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index a5c30cc520..87042dcf3f 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -276,7 +276,7 @@ void incflo::Evolve() amrex::Real time1 = amrex::ParallelDescriptor::second(); // Advance to time t + dt - for (int inonlin = 1; inonlin <= m_adv_iters; ++inonlin) + for (int inonlin = 0; inonlin < m_adv_iters; ++inonlin) do_advance(inonlin); amrex::Print() << std::endl; @@ -312,13 +312,14 @@ void incflo::Evolve() void incflo::do_advance(int inonlin) { - if (m_sim.has_overset() && inonlin == 1) { + if (m_sim.has_overset()) { m_ovst_ops.pre_advance_work(); } - if (m_prescribe_vel && inonlin == 1) { + if (m_prescribe_vel && inonlin == 0) { prescribe_advance(); } else { - if (inonlin == 1) { + amrex::Print() << "INONLIN " << inonlin << std::endl; + if (inonlin == 0) { advance(); } else { amrex::Print() << "Iteration " << inonlin From 519dfdb30138e4c288918ba00ce272efb2b12cbf Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Thu, 12 Sep 2024 15:48:15 -0600 Subject: [PATCH 18/47] Conversion of regular MFIter to fused MFIter --- amr-wind/utilities/console_io.cpp | 30 +++++++++++------------------- 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index 32e200ba5e..2e52935dae 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -240,25 +240,17 @@ void print_nonlinear_residual( level_mask.setVal(1); } - const auto& vnew = velocity_new(lev); - const auto& velstar = vel_star(lev); - auto& veldiff = vel_diff(lev); - - for (amrex::MFIter mfi(velocity_new(lev)); mfi.isValid(); ++mfi) { - const auto bx = mfi.tilebox(); - const auto& velstar_arr = velstar.const_array(mfi); - const auto& veldiff_arr = veldiff.array(mfi); - const auto& velnew_arr = vnew.const_array(mfi); - const auto& levelmask_arr = level_mask.const_array(mfi); - - amrex::ParallelFor( - bx, AMREX_SPACEDIM, - [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept { - veldiff_arr(i, j, k, n) = - (velnew_arr(i, j, k, n) - velstar_arr(i, j, k, n)) * - levelmask_arr(i, j, k); - }); - } + const auto& velnew_arr = velocity_new(lev).const_arrays(); + const auto& velstar_arr = vel_star(lev).const_arrays(); + const auto& veldiff_arr = vel_diff(lev).arrays(); + const auto& levelmask_arr = level_mask.const_arrays(); + amrex::ParallelFor( + velocity_new(lev), amrex::IntVect(0), AMREX_SPACEDIM, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { + veldiff_arr[nbx](i, j, k, n) = (velnew_arr[nbx](i, j, k, n) - + velstar_arr[nbx](i, j, k, n)) * + levelmask_arr[nbx](i, j, k); + }); } amrex::Array rms_vel; From f951253abc61268b9dee6af3ac86f7f287b23999 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Fri, 13 Sep 2024 11:50:28 -0600 Subject: [PATCH 19/47] Clean-up --- amr-wind/incflo_advance.cpp | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index a8a5fc033e..7a2f7dffd6 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -589,21 +589,6 @@ void incflo::ApplyPredictorNonLinear(bool incremental_projection) field.num_comp(), 1); } - // Interpolate the NPH velocity from Old and New states - // Fillpatch the velocity - // change to new_time or remove - // icns().fields().field.fillpatch(0.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 advection of momentum at NPH icns().compute_advection_term(amr_wind::FieldState::NPH); @@ -675,7 +660,6 @@ amr_wind::field_ops::lincomb( } // ScratchField to store the old np1 - // Should nghost be 0 here? auto vel_np1_old = m_repo.create_scratch_field( "vel_np1_old", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL); From a3dfe627aa372175cc3a15ff0f81d8b7c7a82ebe Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Tue, 24 Sep 2024 12:57:43 -0600 Subject: [PATCH 20/47] 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 123904193f..5e111a1d9c 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 87042dcf3f..799b595026 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 7a2f7dffd6..925e2d0906 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; From 2712334e6e775c9ac37bbdf1d50ff55cd01d7fd4 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Tue, 24 Sep 2024 16:24:14 -0600 Subject: [PATCH 21/47] Added the overset mask (mask_cell) to non-linear error calculation --- amr-wind/utilities/console_io.cpp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index 2e52935dae..a908efb56d 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -224,6 +224,8 @@ void print_nonlinear_residual( const auto& mesh = sim.mesh(); const auto& velocity_new = sim.pde_manager().icns().fields().field; + const auto& oset_mask = sim.repo().get_int_field("mask_cell"); + const int ncomp = AMREX_SPACEDIM; for (int lev = 0; lev < nlevels; ++lev) { @@ -244,12 +246,19 @@ void print_nonlinear_residual( const auto& velstar_arr = vel_star(lev).const_arrays(); const auto& veldiff_arr = vel_diff(lev).arrays(); const auto& levelmask_arr = level_mask.const_arrays(); + const auto& osetmask_arr = oset_mask(lev).const_arrays(); + amrex::ParallelFor( velocity_new(lev), amrex::IntVect(0), AMREX_SPACEDIM, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k, int n) noexcept { - veldiff_arr[nbx](i, j, k, n) = (velnew_arr[nbx](i, j, k, n) - - velstar_arr[nbx](i, j, k, n)) * - levelmask_arr[nbx](i, j, k); + if (osetmask_arr[nbx](i, j, k) == 0) { + veldiff_arr[nbx](i, j, k, n) = 0.; + } else { + veldiff_arr[nbx](i, j, k, n) = + (velnew_arr[nbx](i, j, k, n) - + velstar_arr[nbx](i, j, k, n)) * + levelmask_arr[nbx](i, j, k); + } }); } From 2c9d87a3dc77f7391be0b27295eb46186cc0ec4b Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Mon, 30 Sep 2024 14:13:54 -0600 Subject: [PATCH 22/47] Changing output text in print_nonlinear_residual() --- amr-wind/utilities/console_io.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index a908efb56d..5e3f215057 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -271,7 +271,7 @@ void print_nonlinear_residual( } } - amrex::Print() << "Non-linear residual for u " << rms_vel[0] << " v " + amrex::Print() << "Norm of change over advection iterations in u " << rms_vel[0] << " v " << rms_vel[1] << " w " << rms_vel[2] << std::endl; } From 02125fb84626f94e8a349798fbb4e6d7afd1986c Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Wed, 9 Oct 2024 14:13:15 -0600 Subject: [PATCH 23/47] Formatting --- amr-wind/equation_systems/AdvOp_Godunov.H | 15 ++++++++------- amr-wind/setup/init.cpp | 2 +- amr-wind/utilities/console_io.cpp | 5 +++-- 3 files changed, 12 insertions(+), 10 deletions(-) diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 6126c750e1..196f136ee3 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -170,7 +170,8 @@ struct AdvectionOp< } // if state is NPH, then n and n+1 are known, and only // spatial extrapolation is performed - amrex::Real dt_extrap = (fstate == FieldState::NPH) ? 0.0 : dt; + amrex::Real dt_extrap = + (fstate == FieldState::NPH) ? 0.0 : dt; HydroUtils::ComputeFluxesOnBoxFromState( bx, PDE::ndim, mfi, @@ -236,12 +237,12 @@ struct AdvectionOp< Field& w_mac; amrex::Gpu::DeviceVector iconserv; -godunov::scheme godunov_scheme = godunov::scheme::WENOZ; -std::string godunov_type{"weno_z"}; -const bool fluxes_are_area_weighted{false}; -bool godunov_use_forces_in_trans{false}; -bool m_allow_inflow_on_outflow{false}; -std::string advection_type{"Godunov"}; + godunov::scheme godunov_scheme = godunov::scheme::WENOZ; + std::string godunov_type{"weno_z"}; + const bool fluxes_are_area_weighted{false}; + bool godunov_use_forces_in_trans{false}; + bool m_allow_inflow_on_outflow{false}; + std::string advection_type{"Godunov"}; }; } // namespace amr_wind::pde diff --git a/amr-wind/setup/init.cpp b/amr-wind/setup/init.cpp index a5ed713735..5ad166e110 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,0); + ApplyPredictor(true, 0); { auto& vel = icns().fields().field; diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index 5e3f215057..5f73e115f9 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -271,8 +271,9 @@ void print_nonlinear_residual( } } - amrex::Print() << "Norm of change over advection iterations in u " << rms_vel[0] << " v " - << rms_vel[1] << " w " << rms_vel[2] << std::endl; + amrex::Print() << "Norm of change over advection iterations in u " + << rms_vel[0] << " v " << rms_vel[1] << " w " << rms_vel[2] + << std::endl; } } // namespace amr_wind::io From c5dd75b9958e5771b1e620d7bdac08fa3071df35 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Wed, 9 Oct 2024 14:19:26 -0600 Subject: [PATCH 24/47] Changing dt_extrap to const --- amr-wind/equation_systems/AdvOp_Godunov.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 196f136ee3..43231eb05b 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -170,7 +170,7 @@ struct AdvectionOp< } // if state is NPH, then n and n+1 are known, and only // spatial extrapolation is performed - amrex::Real dt_extrap = + const amrex::Real dt_extrap = (fstate == FieldState::NPH) ? 0.0 : dt; HydroUtils::ComputeFluxesOnBoxFromState( From ae81b846d54590fe2ad352c5cba5be4c62a88bfa Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Wed, 9 Oct 2024 14:27:07 -0600 Subject: [PATCH 25/47] amrex::Abort for MOL within the non-linear iteration loop --- amr-wind/incflo_advance.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 925e2d0906..5a9bbc5ee9 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -55,6 +55,9 @@ void incflo::advance(int inonlin) ApplyPredictor(false, inonlin); if (!m_use_godunov) { + if (inonlin > 0) + amrex::Abort("Advection iterations are not supported for MOL"); + ApplyCorrector(); } } From 834f4a26191415dc42ca29f7cb38c2ad82d932af Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Wed, 9 Oct 2024 15:11:26 -0600 Subject: [PATCH 26/47] Fixing scope error for dt_extrap --- amr-wind/equation_systems/icns/icns_advection.H | 8 ++++---- amr-wind/incflo_advance.cpp | 7 +++++-- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 5a128cdbda..76bb6e6184 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -332,13 +332,13 @@ struct AdvectionOp m_allow_inflow_on_outflow = true; } + // if state is NPH, then n and n+1 are known, and only + // spatial extrapolation is performed + const amrex::Real dt_extrap = + (fstate == FieldState::NPH) ? 0.0 : dt; #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - - // if state is NPH, then n and n+1 are known, and only - // spatial extrapolation is performed - amrex::Real dt_extrap = (fstate == FieldState::NPH) ? 0.0 : dt; for (amrex::MFIter mfi(dof_field(lev), mfi_info); mfi.isValid(); ++mfi) { const auto& bx = mfi.tilebox(); diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 5a9bbc5ee9..bd8925a3d5 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -50,13 +50,16 @@ void incflo::pre_advance_stage2() void incflo::advance(int inonlin) { BL_PROFILE("amr-wind::incflo::advance"); - if (inonlin == 0) m_sim.pde_manager().advance_states(); + if (inonlin == 0) { + m_sim.pde_manager().advance_states() + }; ApplyPredictor(false, inonlin); if (!m_use_godunov) { - if (inonlin > 0) + if (inonlin > 0) { amrex::Abort("Advection iterations are not supported for MOL"); + } ApplyCorrector(); } From 5d05f0eeff422fe289bfe1af5d9b17f9a666cfbb Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Wed, 9 Oct 2024 15:18:36 -0600 Subject: [PATCH 27/47] Fixing CI warnings --- amr-wind/incflo_advance.cpp | 2 +- amr-wind/utilities/console_io.H | 1 - amr-wind/utilities/console_io.cpp | 2 -- 3 files changed, 1 insertion(+), 4 deletions(-) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index bd8925a3d5..d87bd7d841 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -51,7 +51,7 @@ 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(false, inonlin); diff --git a/amr-wind/utilities/console_io.H b/amr-wind/utilities/console_io.H index 223ed95b12..ed57676b94 100644 --- a/amr-wind/utilities/console_io.H +++ b/amr-wind/utilities/console_io.H @@ -21,7 +21,6 @@ void print_mlmg_info(const std::string& solve_name, const amrex::MLMG& mlmg); void print_tpls(std::ostream& /*out*/); -// void print_nonlinear_residual(CFDSim& sim); void print_nonlinear_residual( CFDSim& sim, ScratchField& vel_diff, ScratchField& vel_star); diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index 5f73e115f9..e33653016a 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -218,9 +218,7 @@ void print_tpls(std::ostream& out) void print_nonlinear_residual( CFDSim& sim, ScratchField& vel_diff, ScratchField& vel_star) { - // using namespace amrex; const int nlevels = sim.repo().num_active_levels(); - const auto& geom = sim.mesh().Geom(); const auto& mesh = sim.mesh(); const auto& velocity_new = sim.pde_manager().icns().fields().field; From e7d4e2ad9d09c935f0041a3aa0672c1dd65fa253 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Wed, 9 Oct 2024 15:32:53 -0600 Subject: [PATCH 28/47] Conditional operator for fstate --- amr-wind/equation_systems/icns/icns.H | 3 ++- amr-wind/incflo_advance.cpp | 21 ++++++++++++--------- amr-wind/utilities/console_io.cpp | 2 -- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/amr-wind/equation_systems/icns/icns.H b/amr-wind/equation_systems/icns/icns.H index 703cdbf5c4..564c222fa4 100644 --- a/amr-wind/equation_systems/icns/icns.H +++ b/amr-wind/equation_systems/icns/icns.H @@ -43,7 +43,8 @@ struct ICNS : VectorTransport static constexpr bool has_diffusion = true; // No n+1/2 state for velocity for now - static constexpr bool need_nph_state = true; + // static constexpr bool need_nph_state = true; + static constexpr bool need_nph_state = false; }; } // namespace amr_wind::pde diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index d87bd7d841..7bb05f75fc 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -354,15 +354,18 @@ void incflo::ApplyPredictor(bool incremental_projection, int inonlin) } // With scalars computed, compute advection of momentum - 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 - // independently for the right hand side, without 1/rho term + // if (inonlin == 0) { + // icns().compute_advection_term(amr_wind::FieldState::Old); + //} else { + // icns().compute_advection_term(amr_wind::FieldState::NPH); + //} + const auto fstate = + (inonlin == 0) ? amr_wind::FieldState::Old : amr_wind::FieldState::NPH; + icns().compute_advection_term(fstate); + + // ************************************************************************************* + // 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); diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index e33653016a..0f37287794 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -224,8 +224,6 @@ void print_nonlinear_residual( const auto& velocity_new = sim.pde_manager().icns().fields().field; const auto& oset_mask = sim.repo().get_int_field("mask_cell"); - const int ncomp = AMREX_SPACEDIM; - for (int lev = 0; lev < nlevels; ++lev) { amrex::iMultiFab level_mask; From b0a8b7da6ea8645ec8af4ffef9f682e797e49a9b Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Wed, 9 Oct 2024 16:45:53 -0600 Subject: [PATCH 29/47] Increasing the number of fields in the equation_systems/test_pde.cpp unit test by 1 due to the addition of the NPH state in icns.H, which is used for advection iterations --- amr-wind/equation_systems/icns/icns.H | 5 ++--- unit_tests/equation_systems/test_pde.cpp | 4 ++-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/amr-wind/equation_systems/icns/icns.H b/amr-wind/equation_systems/icns/icns.H index 564c222fa4..f6e11b1fd8 100644 --- a/amr-wind/equation_systems/icns/icns.H +++ b/amr-wind/equation_systems/icns/icns.H @@ -42,9 +42,8 @@ struct ICNS : VectorTransport static constexpr bool multiply_rho = true; static constexpr bool has_diffusion = true; - // No n+1/2 state for velocity for now - // static constexpr bool need_nph_state = true; - static constexpr bool need_nph_state = false; + // n+1/2 velocity for advection iterations + static constexpr bool need_nph_state = true; }; } // namespace amr_wind::pde diff --git a/unit_tests/equation_systems/test_pde.cpp b/unit_tests/equation_systems/test_pde.cpp index af23217c34..907e0ac0b9 100644 --- a/unit_tests/equation_systems/test_pde.cpp +++ b/unit_tests/equation_systems/test_pde.cpp @@ -19,7 +19,7 @@ TEST_F(PDETest, test_pde_create_godunov) EXPECT_EQ(pde_mgr.scalar_eqns().size(), 1); - EXPECT_EQ(mesh().field_repo().num_fields(), 21); + EXPECT_EQ(mesh().field_repo().num_fields(), 22); } TEST_F(PDETest, test_pde_create_mol) @@ -35,7 +35,7 @@ TEST_F(PDETest, test_pde_create_mol) EXPECT_EQ(pde_mgr.scalar_eqns().size(), 1); - EXPECT_EQ(mesh().field_repo().num_fields(), 25); + EXPECT_EQ(mesh().field_repo().num_fields(), 26); } } // namespace amr_wind_tests From fb7d9fc9fb981d35d6d669aa3bac0d70b6f646b8 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Thu, 10 Oct 2024 09:33:35 -0600 Subject: [PATCH 30/47] Update amr-wind/equation_systems/icns/icns_advection.H Co-authored-by: Marc T. Henry de Frahan --- amr-wind/equation_systems/icns/icns_advection.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 76bb6e6184..67b284a567 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -192,7 +192,7 @@ struct AdvectionOp // if state is NPH, then n and n+1 are known, and only // spatial extrapolation is performed - amrex::Real dt_extrap = (fstate == FieldState::NPH) ? 0.0 : dt; + const amrex::Real dt_extrap = (fstate == FieldState::NPH) ? 0.0 : dt; for (int lev = 0; lev < repo.num_active_levels(); ++lev) { HydroUtils::ExtrapVelToFaces( dof_field(lev), src_term(lev), u_mac(lev), v_mac(lev), From 710940fa1cf935f59d398e70da99f1cb56ca3bf0 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Thu, 10 Oct 2024 09:34:02 -0600 Subject: [PATCH 31/47] Update amr-wind/incflo_advance.cpp Co-authored-by: Marc T. Henry de Frahan --- amr-wind/incflo_advance.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 7bb05f75fc..28e4821516 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -47,7 +47,7 @@ void incflo::pre_advance_stage2() * * \callgraph */ -void incflo::advance(int inonlin) +void incflo::advance(const int inonlin) { BL_PROFILE("amr-wind::incflo::advance"); if (inonlin == 0) { From 2c8762e11e3a9b06bc0733a0124c0eedaced8e3e Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Thu, 10 Oct 2024 09:34:29 -0600 Subject: [PATCH 32/47] Update amr-wind/incflo.H Co-authored-by: Marc T. Henry de Frahan --- amr-wind/incflo.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 5e111a1d9c..05f12c424e 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -96,7 +96,7 @@ public: bool regrid_and_update(); void pre_advance_stage1(); void pre_advance_stage2(); - void do_advance(int inonlin); + void do_advance(const int inonlin); void post_overset_work(); void advance(int inonlin); void prescribe_advance(); From a7a4bf79f794fd46747210f45cca55de20265ea5 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Thu, 10 Oct 2024 09:37:12 -0600 Subject: [PATCH 33/47] Update amr-wind/incflo.H Co-authored-by: Marc T. Henry de Frahan --- amr-wind/incflo.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 05f12c424e..40373fd7e6 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -134,7 +134,7 @@ public: void ComputeDt(bool explicit_diffusion); void ComputePrescribeDt(); - void ApplyPredictor(bool incremental_projection = false, int inonlin = 0); + void ApplyPredictor(bool incremental_projection = false, const int inonlin = 0); void ApplyCorrector(); void ApplyPrescribeStep(); From fed93e2f9ba5807cd952256cc0efe5608f2651c9 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Thu, 10 Oct 2024 09:37:39 -0600 Subject: [PATCH 34/47] Update amr-wind/incflo.cpp Co-authored-by: Marc T. Henry de Frahan --- amr-wind/incflo.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index 799b595026..a1688f819b 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -310,7 +310,7 @@ void incflo::Evolve() } } -void incflo::do_advance(int inonlin) +void incflo::do_advance(const int inonlin) { if (m_sim.has_overset()) { m_ovst_ops.pre_advance_work(); From ec7d4b29ed7e4c624c6375d54841b9f89b0d5c51 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Thu, 10 Oct 2024 09:37:52 -0600 Subject: [PATCH 35/47] Update amr-wind/incflo_advance.cpp Co-authored-by: Marc T. Henry de Frahan --- amr-wind/incflo_advance.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 28e4821516..9aaba8db8d 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -52,7 +52,7 @@ void incflo::advance(const int inonlin) BL_PROFILE("amr-wind::incflo::advance"); if (inonlin == 0) { m_sim.pde_manager().advance_states(); - }; + } ApplyPredictor(false, inonlin); From d2dadbf89114f7c0b8f298355b8ff084eebd8424 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Thu, 10 Oct 2024 09:38:41 -0600 Subject: [PATCH 36/47] Update amr-wind/incflo_advance.cpp Co-authored-by: Marc T. Henry de Frahan --- amr-wind/incflo_advance.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 9aaba8db8d..c301d2016c 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -170,7 +170,7 @@ void incflo::advance(const int inonlin) *
  • \ref incflo::ApplyProjection "Apply projection" * */ -void incflo::ApplyPredictor(bool incremental_projection, int inonlin) +void incflo::ApplyPredictor(bool incremental_projection, const int inonlin) { BL_PROFILE("amr-wind::incflo::ApplyPredictor"); // We use the new time value for things computed on the "*" state From 6b5298b12673baeebf39834f76bc769b440cc531 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Thu, 10 Oct 2024 09:39:59 -0600 Subject: [PATCH 37/47] Update amr-wind/utilities/console_io.H Co-authored-by: Marc T. Henry de Frahan --- amr-wind/utilities/console_io.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/utilities/console_io.H b/amr-wind/utilities/console_io.H index ed57676b94..6b3cdabdc6 100644 --- a/amr-wind/utilities/console_io.H +++ b/amr-wind/utilities/console_io.H @@ -22,7 +22,7 @@ void print_mlmg_info(const std::string& solve_name, const amrex::MLMG& mlmg); void print_tpls(std::ostream& /*out*/); void print_nonlinear_residual( - CFDSim& sim, ScratchField& vel_diff, ScratchField& vel_star); + const CFDSim& sim, ScratchField& vel_diff, const ScratchField& vel_star); } // namespace amr_wind::io From 3472ccadbf2627d717f8f87ef5385d329a2d672c Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Thu, 10 Oct 2024 09:46:08 -0600 Subject: [PATCH 38/47] const arguments for print_nonlinear_residual() --- amr-wind/utilities/console_io.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index 0f37287794..3047450cdb 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -216,7 +216,7 @@ void print_tpls(std::ostream& out) } void print_nonlinear_residual( - CFDSim& sim, ScratchField& vel_diff, ScratchField& vel_star) + const CFDSim& sim, ScratchField& vel_diff, const ScratchField& vel_star) { const int nlevels = sim.repo().num_active_levels(); const auto& mesh = sim.mesh(); From c17d8df4e37b571b5a5a5b4b948be24ead0f5611 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Thu, 10 Oct 2024 09:56:18 -0600 Subject: [PATCH 39/47] Ternary expression for the field state that goes into pre_advection_actions() --- amr-wind/incflo_advance.cpp | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index c301d2016c..87cbd2723c 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -282,11 +282,9 @@ void incflo::ApplyPredictor(bool incremental_projection, const int inonlin) } // Extrapolate and apply MAC projection for advection velocities - if (inonlin == 0) { - icns().pre_advection_actions(amr_wind::FieldState::Old); - } else { - icns().pre_advection_actions(amr_wind::FieldState::NPH); - } + const auto fstate_preadv = + (inonlin == 0) ? amr_wind::FieldState::Old : amr_wind::FieldState::NPH; + icns().pre_advection_actions(fstate_preadv); // For scalars only first // ************************************************************************************* @@ -354,11 +352,6 @@ void incflo::ApplyPredictor(bool incremental_projection, const int inonlin) } // With scalars computed, compute advection of momentum - // if (inonlin == 0) { - // icns().compute_advection_term(amr_wind::FieldState::Old); - //} else { - // icns().compute_advection_term(amr_wind::FieldState::NPH); - //} const auto fstate = (inonlin == 0) ? amr_wind::FieldState::Old : amr_wind::FieldState::NPH; icns().compute_advection_term(fstate); From cd9d165d50d41b33059e9a12233cd12641e856b7 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Thu, 10 Oct 2024 10:05:31 -0600 Subject: [PATCH 40/47] Array initialization for rms_vel in console_io.cpp --- amr-wind/equation_systems/icns/icns_advection.H | 3 ++- amr-wind/incflo.H | 3 ++- amr-wind/utilities/console_io.cpp | 3 +-- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 67b284a567..77d61fabab 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -192,7 +192,8 @@ struct AdvectionOp // if state is NPH, then n and n+1 are known, and only // spatial extrapolation is performed - const amrex::Real dt_extrap = (fstate == FieldState::NPH) ? 0.0 : dt; + const amrex::Real dt_extrap = + (fstate == FieldState::NPH) ? 0.0 : dt; for (int lev = 0; lev < repo.num_active_levels(); ++lev) { HydroUtils::ExtrapVelToFaces( dof_field(lev), src_term(lev), u_mac(lev), v_mac(lev), diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 40373fd7e6..694b00a931 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -134,7 +134,8 @@ public: void ComputeDt(bool explicit_diffusion); void ComputePrescribeDt(); - void ApplyPredictor(bool incremental_projection = false, const int inonlin = 0); + void + ApplyPredictor(bool incremental_projection = false, const int inonlin = 0); void ApplyCorrector(); void ApplyPrescribeStep(); diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index 3047450cdb..4eb3f9920a 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -258,9 +258,8 @@ void print_nonlinear_residual( }); } - amrex::Array rms_vel; + amrex::Array rms_vel = {0.0}; - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) rms_vel[idim] = 0.; for (int lev = 0; lev < nlevels; ++lev) { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { rms_vel[idim] += vel_diff(lev).norm2(idim); From d78118da19279bac03719371042180961a033ce0 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Thu, 10 Oct 2024 10:17:11 -0600 Subject: [PATCH 41/47] Update amr-wind/utilities/console_io.cpp Co-authored-by: Marc T. Henry de Frahan --- amr-wind/utilities/console_io.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index 4eb3f9920a..a61dda7e16 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -266,9 +266,9 @@ void print_nonlinear_residual( } } - amrex::Print() << "Norm of change over advection iterations in u " - << rms_vel[0] << " v " << rms_vel[1] << " w " << rms_vel[2] - << std::endl; + amrex::Print() << "Norm of change over advection iterations in u: " + << rms_vel[0] << ", v: " << rms_vel[1] + << ", w: " << rms_vel[2] << std::endl; } } // namespace amr_wind::io From 2d87f103940ed72b08cfbb8549d21fef720b5d92 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Thu, 10 Oct 2024 10:31:33 -0600 Subject: [PATCH 42/47] Replacing amrex::Abort() instances with AMREX_ALWAYS_ASSERT_WITH_MESSAGE() --- amr-wind/incflo.cpp | 6 +++--- amr-wind/incflo_advance.cpp | 5 ++--- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index a1688f819b..f29f000e25 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -378,9 +378,9 @@ void incflo::init_physics_and_pde() // Get number of advection iterations pp.query("advection_iterations", m_adv_iters); - if (m_adv_iters < 1) - amrex::Abort( - "The number of advection iterations cannot be less than 1"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + m_adv_iters > 0, + "The number of advection iterations cannot be less than 1"); } auto& pde_mgr = m_sim.pde_manager(); diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 87cbd2723c..904fa54eea 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -57,9 +57,8 @@ void incflo::advance(const int inonlin) ApplyPredictor(false, inonlin); if (!m_use_godunov) { - if (inonlin > 0) { - amrex::Abort("Advection iterations are not supported for MOL"); - } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + inonlin == 0, "Advection iterations are not supported for MOL"); ApplyCorrector(); } From 195046e97d1fe9ab66468f1aec465fc051f92af0 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Fri, 11 Oct 2024 11:13:54 -0600 Subject: [PATCH 43/47] Removing outdated function prototype --- amr-wind/incflo.H | 1 - 1 file changed, 1 deletion(-) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 694b00a931..68a1b4ee19 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -97,7 +97,6 @@ public: void pre_advance_stage1(); void pre_advance_stage2(); void do_advance(const int inonlin); - void post_overset_work(); void advance(int inonlin); void prescribe_advance(); void post_advance_work(); From 56d6e63b674db6d08eb4a25b1572b5e204c3b525 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu <114435459+itopcuoglu@users.noreply.github.com> Date: Fri, 11 Oct 2024 13:46:37 -0600 Subject: [PATCH 44/47] Update amr-wind/incflo_advance.cpp Co-authored-by: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> --- amr-wind/incflo_advance.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index 904fa54eea..c7d78529da 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -415,7 +415,7 @@ void incflo::ApplyPredictor(bool incremental_projection, const int inonlin) PrintMaxVelLocations("after nodal projection"); } // ScratchField to store the old np1 - if (inonlin > 0) { + if (inonlin > 0 && m_verbose > 0) { auto vel_np1_old = m_repo.create_scratch_field( "vel_np1_old", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL); From c2634598b18779421870f4a23ae43aabc8aa7deb Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Fri, 11 Oct 2024 17:06:06 -0600 Subject: [PATCH 45/47] Renaming fixed point iteration indices --- amr-wind/incflo.H | 7 ++++--- amr-wind/incflo.cpp | 13 +++++++------ amr-wind/incflo_advance.cpp | 26 +++++++++++++++----------- amr-wind/utilities/console_io.cpp | 2 +- 4 files changed, 27 insertions(+), 21 deletions(-) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 68a1b4ee19..c35f15c307 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -96,7 +96,7 @@ public: bool regrid_and_update(); void pre_advance_stage1(); void pre_advance_stage2(); - void do_advance(const int inonlin); + void do_advance(const int ifixed_point_iteration); void advance(int inonlin); void prescribe_advance(); void post_advance_work(); @@ -133,8 +133,9 @@ public: void ComputeDt(bool explicit_diffusion); void ComputePrescribeDt(); - void - ApplyPredictor(bool incremental_projection = false, const int inonlin = 0); + void ApplyPredictor( + bool incremental_projection = false, + const int ifixed_point_iteration = 0); void ApplyCorrector(); void ApplyPrescribeStep(); diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index f29f000e25..7f27c0687c 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -276,8 +276,9 @@ void incflo::Evolve() amrex::Real time1 = amrex::ParallelDescriptor::second(); // Advance to time t + dt - for (int inonlin = 0; inonlin < m_adv_iters; ++inonlin) - do_advance(inonlin); + for (int iadvection_iteration = 0; iadvection_iteration < m_adv_iters; + ++iadvection_iteration) + do_advance(iadvection_iteration); amrex::Print() << std::endl; amrex::Real time2 = amrex::ParallelDescriptor::second(); @@ -310,17 +311,17 @@ void incflo::Evolve() } } -void incflo::do_advance(const int inonlin) +void incflo::do_advance(const int ifixed_point_iteration) { if (m_sim.has_overset()) { m_ovst_ops.pre_advance_work(); } - if (m_prescribe_vel && inonlin == 0) { + if (m_prescribe_vel && ifixed_point_iteration == 0) { prescribe_advance(); } else { - amrex::Print() << "Iteration " << inonlin + amrex::Print() << "Iteration " << ifixed_point_iteration << " in advection iteration loop" << std::endl; - advance(inonlin); + advance(ifixed_point_iteration); } 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 c7d78529da..cbedc0625c 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -47,18 +47,19 @@ void incflo::pre_advance_stage2() * * \callgraph */ -void incflo::advance(const int inonlin) +void incflo::advance(const int ifixed_point_iteration) { BL_PROFILE("amr-wind::incflo::advance"); - if (inonlin == 0) { + if (ifixed_point_iteration == 0) { m_sim.pde_manager().advance_states(); } - ApplyPredictor(false, inonlin); + ApplyPredictor(false, ifixed_point_iteration); if (!m_use_godunov) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - inonlin == 0, "Advection iterations are not supported for MOL"); + ifixed_point_iteration == 0, + "Advection iterations are not supported for MOL"); ApplyCorrector(); } @@ -169,7 +170,8 @@ void incflo::advance(const int inonlin) *
  • \ref incflo::ApplyProjection "Apply projection" * */ -void incflo::ApplyPredictor(bool incremental_projection, const int inonlin) +void incflo::ApplyPredictor( + bool incremental_projection, const int ifixed_point_iteration) { BL_PROFILE("amr-wind::incflo::ApplyPredictor"); // We use the new time value for things computed on the "*" state @@ -193,7 +195,7 @@ void incflo::ApplyPredictor(bool incremental_projection, const int inonlin) const auto& density_old = density_new.state(amr_wind::FieldState::Old); auto& density_nph = density_new.state(amr_wind::FieldState::NPH); - if (inonlin > 0) { + if (ifixed_point_iteration > 0) { icns().fields().field.fillpatch(m_time.current_time()); // Get n + 1/2 velocity amr_wind::field_ops::lincomb( @@ -281,8 +283,9 @@ void incflo::ApplyPredictor(bool incremental_projection, const int inonlin) } // Extrapolate and apply MAC projection for advection velocities - const auto fstate_preadv = - (inonlin == 0) ? amr_wind::FieldState::Old : amr_wind::FieldState::NPH; + const auto fstate_preadv = (ifixed_point_iteration == 0) + ? amr_wind::FieldState::Old + : amr_wind::FieldState::NPH; icns().pre_advection_actions(fstate_preadv); // For scalars only first @@ -351,8 +354,9 @@ void incflo::ApplyPredictor(bool incremental_projection, const int inonlin) } // With scalars computed, compute advection of momentum - const auto fstate = - (inonlin == 0) ? amr_wind::FieldState::Old : amr_wind::FieldState::NPH; + const auto fstate = (ifixed_point_iteration == 0) + ? amr_wind::FieldState::Old + : amr_wind::FieldState::NPH; icns().compute_advection_term(fstate); // ************************************************************************************* @@ -415,7 +419,7 @@ void incflo::ApplyPredictor(bool incremental_projection, const int inonlin) PrintMaxVelLocations("after nodal projection"); } // ScratchField to store the old np1 - if (inonlin > 0 && m_verbose > 0) { + if (ifixed_point_iteration > 0 && m_verbose > 0) { auto vel_np1_old = m_repo.create_scratch_field( "vel_np1_old", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL); diff --git a/amr-wind/utilities/console_io.cpp b/amr-wind/utilities/console_io.cpp index a61dda7e16..2c16191c20 100644 --- a/amr-wind/utilities/console_io.cpp +++ b/amr-wind/utilities/console_io.cpp @@ -266,7 +266,7 @@ void print_nonlinear_residual( } } - amrex::Print() << "Norm of change over advection iterations in u: " + amrex::Print() << "Norm of change over fixed point iterations in u: " << rms_vel[0] << ", v: " << rms_vel[1] << ", w: " << rms_vel[2] << std::endl; } From 91fefe16e405f16e950fcff619f97e0e92662902 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Mon, 14 Oct 2024 11:20:28 -0600 Subject: [PATCH 46/47] Changing the entire naming convention to fixed_point_iterations --- amr-wind/incflo.H | 4 ++-- amr-wind/incflo.cpp | 17 +++++++++-------- .../abl_godunov_wenoz_fixedpt.inp} | 2 +- 3 files changed, 12 insertions(+), 11 deletions(-) rename test/test_files/{abl_godunov_wenoz_nl/abl_godunov_wenoz_nl.inp => abl_godunov_wenoz_fixedpt/abl_godunov_wenoz_fixedpt.inp} (98%) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index c35f15c307..e2ec8a54d2 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -174,8 +174,8 @@ private: // Prescribe advection velocity bool m_prescribe_vel = false; - // Advection iterations every timestep - int m_adv_iters{1}; + // Fixed point iterations every timestep + int m_fixed_pt_iters{1}; //! number of cells on all levels including covered cells amrex::Long m_cell_count{-1}; diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index 7f27c0687c..a71cfdfbdb 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -276,9 +276,10 @@ void incflo::Evolve() amrex::Real time1 = amrex::ParallelDescriptor::second(); // Advance to time t + dt - for (int iadvection_iteration = 0; iadvection_iteration < m_adv_iters; - ++iadvection_iteration) - do_advance(iadvection_iteration); + for (int ifixed_point_iteration = 0; + ifixed_point_iteration < m_fixed_pt_iters; + ++ifixed_point_iteration) + do_advance(ifixed_point_iteration); amrex::Print() << std::endl; amrex::Real time2 = amrex::ParallelDescriptor::second(); @@ -319,8 +320,8 @@ void incflo::do_advance(const int ifixed_point_iteration) if (m_prescribe_vel && ifixed_point_iteration == 0) { prescribe_advance(); } else { - amrex::Print() << "Iteration " << ifixed_point_iteration - << " in advection iteration loop" << std::endl; + amrex::Print() << "Fixed point iteration " << ifixed_point_iteration + << std::endl; advance(ifixed_point_iteration); } if (m_sim.has_overset()) { @@ -378,10 +379,10 @@ void incflo::init_physics_and_pde() } // Get number of advection iterations - pp.query("advection_iterations", m_adv_iters); + pp.query("fixed_point_iterations", m_fixed_pt_iters); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_adv_iters > 0, - "The number of advection iterations cannot be less than 1"); + m_fixed_pt_iters > 0, + "The number of fixed point iterations cannot be less than 1"); } auto& pde_mgr = m_sim.pde_manager(); diff --git a/test/test_files/abl_godunov_wenoz_nl/abl_godunov_wenoz_nl.inp b/test/test_files/abl_godunov_wenoz_fixedpt/abl_godunov_wenoz_fixedpt.inp similarity index 98% rename from test/test_files/abl_godunov_wenoz_nl/abl_godunov_wenoz_nl.inp rename to test/test_files/abl_godunov_wenoz_fixedpt/abl_godunov_wenoz_fixedpt.inp index 0a3d6c92bc..6e8ab161cf 100644 --- a/test/test_files/abl_godunov_wenoz_nl/abl_godunov_wenoz_nl.inp +++ b/test/test_files/abl_godunov_wenoz_fixedpt/abl_godunov_wenoz_fixedpt.inp @@ -25,7 +25,7 @@ incflo.density = 1.0 # Reference density incflo.use_godunov = 1 incflo.godunov_type = "weno_z" incflo.diffusion_type = 2 -incflo.advection_iterations = 10 +incflo.fixed_point_iterations = 10 transport.viscosity = 1.0e-5 transport.laminar_prandtl = 0.7 transport.turbulent_prandtl = 0.3333 From 4b993ca7b0b74d94f7f1910a6ffa87b419e04ac8 Mon Sep 17 00:00:00 2001 From: Ilker Topcuoglu Date: Mon, 14 Oct 2024 13:17:45 -0600 Subject: [PATCH 47/47] More changes in names of variables related to fixed point iteration --- amr-wind/incflo.H | 10 +++++----- amr-wind/incflo.cpp | 21 ++++++++++----------- amr-wind/incflo_advance.cpp | 20 ++++++++++---------- 3 files changed, 25 insertions(+), 26 deletions(-) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index e2ec8a54d2..dc8174c4b6 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -96,8 +96,8 @@ public: bool regrid_and_update(); void pre_advance_stage1(); void pre_advance_stage2(); - void do_advance(const int ifixed_point_iteration); - void advance(int inonlin); + void do_advance(const int fixed_point_iteration); + void advance(const int fixed_point_iteration); void prescribe_advance(); void post_advance_work(); @@ -134,8 +134,8 @@ public: void ComputePrescribeDt(); void ApplyPredictor( - bool incremental_projection = false, - const int ifixed_point_iteration = 0); + const bool incremental_projection = false, + const int fixed_point_iteration = 0); void ApplyCorrector(); void ApplyPrescribeStep(); @@ -175,7 +175,7 @@ private: bool m_prescribe_vel = false; // Fixed point iterations every timestep - int m_fixed_pt_iters{1}; + int m_fixed_point_iterations{1}; //! number of cells on all levels including covered cells amrex::Long m_cell_count{-1}; diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index a71cfdfbdb..fbcfa4b656 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -276,10 +276,10 @@ void incflo::Evolve() amrex::Real time1 = amrex::ParallelDescriptor::second(); // Advance to time t + dt - for (int ifixed_point_iteration = 0; - ifixed_point_iteration < m_fixed_pt_iters; - ++ifixed_point_iteration) - do_advance(ifixed_point_iteration); + for (int fixed_point_iteration = 0; + fixed_point_iteration < m_fixed_point_iterations; + ++fixed_point_iteration) + do_advance(fixed_point_iteration); amrex::Print() << std::endl; amrex::Real time2 = amrex::ParallelDescriptor::second(); @@ -312,17 +312,17 @@ void incflo::Evolve() } } -void incflo::do_advance(const int ifixed_point_iteration) +void incflo::do_advance(const int fixed_point_iteration) { if (m_sim.has_overset()) { m_ovst_ops.pre_advance_work(); } - if (m_prescribe_vel && ifixed_point_iteration == 0) { + if (m_prescribe_vel && fixed_point_iteration == 0) { prescribe_advance(); } else { - amrex::Print() << "Fixed point iteration " << ifixed_point_iteration + amrex::Print() << "Fixed point iteration " << fixed_point_iteration << std::endl; - advance(ifixed_point_iteration); + advance(fixed_point_iteration); } if (m_sim.has_overset()) { m_ovst_ops.post_advance_work(); @@ -378,10 +378,9 @@ void incflo::init_physics_and_pde() m_sim.activate_overset(); } - // Get number of advection iterations - pp.query("fixed_point_iterations", m_fixed_pt_iters); + pp.query("fixed_point_iterations", m_fixed_point_iterations); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - m_fixed_pt_iters > 0, + m_fixed_point_iterations > 0, "The number of fixed point iterations cannot be less than 1"); } diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index cbedc0625c..a78480f1d2 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -47,19 +47,19 @@ void incflo::pre_advance_stage2() * * \callgraph */ -void incflo::advance(const int ifixed_point_iteration) +void incflo::advance(const int fixed_point_iteration) { BL_PROFILE("amr-wind::incflo::advance"); - if (ifixed_point_iteration == 0) { + if (fixed_point_iteration == 0) { m_sim.pde_manager().advance_states(); } - ApplyPredictor(false, ifixed_point_iteration); + ApplyPredictor(false, fixed_point_iteration); if (!m_use_godunov) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - ifixed_point_iteration == 0, - "Advection iterations are not supported for MOL"); + fixed_point_iteration == 0, + "Fixed point iterations are not supported for MOL"); ApplyCorrector(); } @@ -171,7 +171,7 @@ void incflo::advance(const int ifixed_point_iteration) * */ void incflo::ApplyPredictor( - bool incremental_projection, const int ifixed_point_iteration) + const bool incremental_projection, const int fixed_point_iteration) { BL_PROFILE("amr-wind::incflo::ApplyPredictor"); // We use the new time value for things computed on the "*" state @@ -195,7 +195,7 @@ void incflo::ApplyPredictor( const auto& density_old = density_new.state(amr_wind::FieldState::Old); auto& density_nph = density_new.state(amr_wind::FieldState::NPH); - if (ifixed_point_iteration > 0) { + if (fixed_point_iteration > 0) { icns().fields().field.fillpatch(m_time.current_time()); // Get n + 1/2 velocity amr_wind::field_ops::lincomb( @@ -283,7 +283,7 @@ void incflo::ApplyPredictor( } // Extrapolate and apply MAC projection for advection velocities - const auto fstate_preadv = (ifixed_point_iteration == 0) + const auto fstate_preadv = (fixed_point_iteration == 0) ? amr_wind::FieldState::Old : amr_wind::FieldState::NPH; icns().pre_advection_actions(fstate_preadv); @@ -354,7 +354,7 @@ void incflo::ApplyPredictor( } // With scalars computed, compute advection of momentum - const auto fstate = (ifixed_point_iteration == 0) + const auto fstate = (fixed_point_iteration == 0) ? amr_wind::FieldState::Old : amr_wind::FieldState::NPH; icns().compute_advection_term(fstate); @@ -419,7 +419,7 @@ void incflo::ApplyPredictor( PrintMaxVelLocations("after nodal projection"); } // ScratchField to store the old np1 - if (ifixed_point_iteration > 0 && m_verbose > 0) { + if (fixed_point_iteration > 0 && m_verbose > 0) { auto vel_np1_old = m_repo.create_scratch_field( "vel_np1_old", AMREX_SPACEDIM, 1, amr_wind::FieldLoc::CELL);