From 935b17e63d282e89ea661660bb822030b66de804 Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Mon, 21 Oct 2024 11:00:48 -0600 Subject: [PATCH] Set up fill_boundaries as a step prior to ApplyPredictor (#1274) * stop setting mac bndry cells to zero a priori * do MAC fillpatch before step * set up pre_predictor_work to update data after first read * set up time switch for read_file * grouping calls at beginning of time step * add intended interp time * do a no-op for MAC fill the second time by passing some foextrap, some ext_dir BCs --- amr-wind/core/Field.cpp | 3 +- amr-wind/core/FieldFillPatchOps.H | 24 ++++---- amr-wind/core/Physics.H | 3 + amr-wind/equation_systems/AdvOp_Godunov.H | 11 ++++ amr-wind/equation_systems/PDEBase.H | 4 ++ amr-wind/equation_systems/PDEBase.cpp | 55 +++++++++++++++++++ amr-wind/equation_systems/icns/icns.H | 1 - .../equation_systems/icns/icns_advection.H | 55 ++++++++++--------- .../equation_systems/vof/vof_momentum_flux.H | 49 ++++++++++------- amr-wind/incflo.H | 1 + amr-wind/incflo.cpp | 4 ++ amr-wind/incflo_advance.cpp | 13 ++++- amr-wind/wind_energy/ABL.H | 2 + amr-wind/wind_energy/ABL.cpp | 6 ++ amr-wind/wind_energy/ABLBoundaryPlane.H | 12 +++- amr-wind/wind_energy/ABLBoundaryPlane.cpp | 19 ++++++- amr-wind/wind_energy/ABLFillInflow.H | 1 + amr-wind/wind_energy/ABLFillInflow.cpp | 36 ++++++++---- amr-wind/wind_energy/ABLFillMPL.H | 1 + amr-wind/wind_energy/ABLFillMPL.cpp | 3 +- 20 files changed, 227 insertions(+), 76 deletions(-) diff --git a/amr-wind/core/Field.cpp b/amr-wind/core/Field.cpp index 766c72aa14..889e66b5ed 100644 --- a/amr-wind/core/Field.cpp +++ b/amr-wind/core/Field.cpp @@ -230,7 +230,8 @@ void Field::fillpatch_sibling_fields( } fop.fillpatch_sibling_fields( - lev, time, mfabs, mfabs, cfabs, ng, m_info->m_bcrec, field_state()); + lev, time, mfabs, mfabs, cfabs, ng, m_info->m_bcrec, + m_info->m_bcrec, field_state()); } } diff --git a/amr-wind/core/FieldFillPatchOps.H b/amr-wind/core/FieldFillPatchOps.H index f6f58da440..5de870c86b 100644 --- a/amr-wind/core/FieldFillPatchOps.H +++ b/amr-wind/core/FieldFillPatchOps.H @@ -60,7 +60,8 @@ public: amrex::Array& ffabs, amrex::Array& cfabs, const amrex::IntVect& nghost, - const amrex::Vector& bcrec, + const amrex::Vector& fillpatch_bcrec, + const amrex::Vector& physbc_bcrec, const FieldState fstate = FieldState::New) = 0; //! Implementation that handles filling patches from a coarse to fine level @@ -121,6 +122,7 @@ public: amrex::Array& /*cfabs*/, const amrex::IntVect& /*nghost*/, const amrex::Vector& /*bcrec*/, + const amrex::Vector& /*bcrec*/, const FieldState /*fstate*/) override { for (const auto& mfab : mfabs) { @@ -272,14 +274,15 @@ public: amrex::Array& ffabs, amrex::Array& cfabs, const amrex::IntVect& nghost, - const amrex::Vector& bcrec, + const amrex::Vector& fillpatch_bcrec, + const amrex::Vector& physbc_bcrec, const FieldState /*fstate = FieldState::New*/) override { if (lev == 0) { for (int i = 0; i < static_cast(mfabs.size()); i++) { amrex::PhysBCFunct> physbc( - m_mesh.Geom(lev), bcrec, bc_functor_face(i)); + m_mesh.Geom(lev), physbc_bcrec, bc_functor_face(i)); amrex::FillPatchSingleLevel( *mfabs[i], nghost, time, {ffabs[i]}, {time}, 0, 0, 1, m_mesh.Geom(lev), physbc, i); @@ -287,19 +290,19 @@ public: } else { AMREX_D_TERM( amrex::PhysBCFunct> cphysbcx( - m_mesh.Geom(lev - 1), bcrec, bc_functor_face(0)); + m_mesh.Geom(lev - 1), physbc_bcrec, bc_functor_face(0)); , amrex::PhysBCFunct> cphysbcy( - m_mesh.Geom(lev - 1), bcrec, bc_functor_face(1)); + m_mesh.Geom(lev - 1), physbc_bcrec, bc_functor_face(1)); , amrex::PhysBCFunct> cphysbcz( - m_mesh.Geom(lev - 1), bcrec, bc_functor_face(2));); + m_mesh.Geom(lev - 1), physbc_bcrec, bc_functor_face(2));); AMREX_D_TERM( amrex::PhysBCFunct> fphysbcx( - m_mesh.Geom(lev), bcrec, bc_functor_face(0)); + m_mesh.Geom(lev), physbc_bcrec, bc_functor_face(0)); , amrex::PhysBCFunct> fphysbcy( - m_mesh.Geom(lev), bcrec, bc_functor_face(1)); + m_mesh.Geom(lev), physbc_bcrec, bc_functor_face(1)); , amrex::PhysBCFunct> fphysbcz( - m_mesh.Geom(lev), bcrec, bc_functor_face(2));); + m_mesh.Geom(lev), physbc_bcrec, bc_functor_face(2));); amrex::Array< amrex::PhysBCFunct>, @@ -313,7 +316,8 @@ public: amrex::Array idx = {AMREX_D_DECL(0, 1, 2)}; const amrex::Array, AMREX_SPACEDIM> - bcrec_arr = {AMREX_D_DECL(bcrec, bcrec, bcrec)}; + bcrec_arr = {AMREX_D_DECL( + fillpatch_bcrec, fillpatch_bcrec, fillpatch_bcrec)}; amrex::FillPatchTwoLevels( mfabs, nghost, time, {cfabs}, {time}, {ffabs}, {time}, 0, 0, 1, diff --git a/amr-wind/core/Physics.H b/amr-wind/core/Physics.H index ce8ed74d2e..d34fb21b0c 100644 --- a/amr-wind/core/Physics.H +++ b/amr-wind/core/Physics.H @@ -79,6 +79,9 @@ public: //! Perform tasks necessary before advancing timestep virtual void pre_advance_work() = 0; + //! Perform tasks necessary only once per timestep, after pre_advance + virtual void pre_predictor_work() {} + //! Perform tasks necessary after advancing timestep virtual void post_advance_work() = 0; diff --git a/amr-wind/equation_systems/AdvOp_Godunov.H b/amr-wind/equation_systems/AdvOp_Godunov.H index 43231eb05b..a43daf9b28 100644 --- a/amr-wind/equation_systems/AdvOp_Godunov.H +++ b/amr-wind/equation_systems/AdvOp_Godunov.H @@ -91,6 +91,7 @@ struct AdvectionOp< // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term; const auto& dof_field = fields.field.state(fstate); + const auto& dof_nph = fields.field.state(amr_wind::FieldState::NPH); auto flux_x = repo.create_scratch_field(PDE::ndim, 0, amr_wind::FieldLoc::XFACE); @@ -107,6 +108,7 @@ struct AdvectionOp< // only needed if multiplying by rho below const auto& den = density.state(fstate); + const auto& den_nph = density.state(amr_wind::FieldState::NPH); for (int lev = 0; lev < repo.num_active_levels(); ++lev) { amrex::MFItInfo mfi_info; @@ -121,15 +123,21 @@ struct AdvectionOp< ++mfi) { const auto& bx = mfi.tilebox(); auto rho_arr = den(lev).array(mfi); + auto rho_nph_arr = den_nph(lev).array(mfi); auto tra_arr = dof_field(lev).array(mfi); + auto tra_nph_arr = dof_nph(lev).array(mfi); amrex::FArrayBox rhotracfab; amrex::Array4 rhotrac; + amrex::FArrayBox rhotrac_nph_fab; + amrex::Array4 rhotrac_nph; if (PDE::multiply_rho) { auto rhotrac_box = amrex::grow(bx, fvm::Godunov::nghost_state); rhotracfab.resize(rhotrac_box, PDE::ndim); rhotrac = rhotracfab.array(); + rhotrac_nph_fab.resize(rhotrac_box, PDE::ndim); + rhotrac_nph = rhotrac_nph_fab.array(); amrex::ParallelFor( rhotrac_box, PDE::ndim, @@ -137,6 +145,8 @@ struct AdvectionOp< int i, int j, int k, int n) noexcept { rhotrac(i, j, k, n) = rho_arr(i, j, k) * tra_arr(i, j, k, n); + rhotrac_nph(i, j, k, n) = + rho_nph_arr(i, j, k) * tra_nph_arr(i, j, k, n); }); } @@ -176,6 +186,7 @@ struct AdvectionOp< HydroUtils::ComputeFluxesOnBoxFromState( bx, PDE::ndim, mfi, (PDE::multiply_rho ? rhotrac : tra_arr), + (PDE::multiply_rho ? rhotrac_nph : tra_nph_arr), (*flux_x)(lev).array(mfi), (*flux_y)(lev).array(mfi), (*flux_z)(lev).array(mfi), (*face_x)(lev).array(mfi), (*face_y)(lev).array(mfi), (*face_z)(lev).array(mfi), diff --git a/amr-wind/equation_systems/PDEBase.H b/amr-wind/equation_systems/PDEBase.H index 5abbd19278..f54d28324f 100644 --- a/amr-wind/equation_systems/PDEBase.H +++ b/amr-wind/equation_systems/PDEBase.H @@ -132,6 +132,10 @@ public: //! Advance states for all registered PDEs at the beginning of a timestep void advance_states(); + //! Prepare boundaries for registered PDEs at the beginning of a timestep: + //! -- intended for when boundaries rely on intermediate time states (nph) + void prepare_boundaries(); + //! Call fillpatch operator on state variables for all registered PDEs void fillpatch_state_fields( const amrex::Real time, const FieldState fstate = FieldState::New); diff --git a/amr-wind/equation_systems/PDEBase.cpp b/amr-wind/equation_systems/PDEBase.cpp index 10775cbb61..aa89872ff4 100644 --- a/amr-wind/equation_systems/PDEBase.cpp +++ b/amr-wind/equation_systems/PDEBase.cpp @@ -4,6 +4,7 @@ #include "amr-wind/core/FieldRepo.H" #include "amr-wind/equation_systems/PDEHelpers.H" #include "amr-wind/incflo_enums.H" +#include "amr-wind/core/field_ops.H" #include "AMReX_ParmParse.H" @@ -87,6 +88,60 @@ void PDEMgr::advance_states() } } +void PDEMgr::prepare_boundaries() +{ + // If state variables exist at NPH, fill their boundary cells + const auto nph_time = + m_sim.time().current_time() + 0.5 * m_sim.time().delta_t(); + if (m_constant_density && + m_sim.repo().field_exists("density", FieldState::NPH)) { + auto& nph_field = + m_sim.repo().get_field("density").state(FieldState::NPH); + auto& new_field = + m_sim.repo().get_field("density").state(FieldState::New); + // Need to check if this step is necessary + // Guessing that for no-op dirichlet bcs, the new needs to be copied to + // NPH + field_ops::copy( + nph_field, new_field, 0, 0, new_field.num_comp(), + new_field.num_grow()); + } + + if (m_sim.repo().field_exists( + icns().fields().field.base_name(), FieldState::NPH)) { + auto& nph_field = icns().fields().field.state(FieldState::NPH); + auto& new_field = icns().fields().field.state(FieldState::New); + field_ops::copy( + nph_field, new_field, 0, 0, new_field.num_comp(), + new_field.num_grow()); + // Insert fill physical boundary condition call here + } + for (auto& eqn : scalar_eqns()) { + eqn->fields().field.advance_states(); + if (m_sim.repo().field_exists( + eqn->fields().field.base_name(), FieldState::NPH)) { + auto& nph_field = eqn->fields().field.state(FieldState::NPH); + auto& new_field = eqn->fields().field.state(FieldState::New); + field_ops::copy( + nph_field, new_field, 0, 0, new_field.num_comp(), + new_field.num_grow()); + // Insert fill physical boundary condition call here + } + } + + // Fill mac velocities (ghost cells) using velocity BCs + auto& u_mac = m_sim.repo().get_field("u_mac"); + auto& v_mac = m_sim.repo().get_field("v_mac"); + auto& w_mac = m_sim.repo().get_field("w_mac"); + const amrex::IntVect zeros{0, 0, 0}; + if (u_mac.num_grow() > zeros) { + amrex::Array mac_vel = { + AMREX_D_DECL(&u_mac, &v_mac, &w_mac)}; + icns().fields().field.fillpatch_sibling_fields( + nph_time, u_mac.num_grow(), mac_vel); + } +} + void PDEMgr::fillpatch_state_fields( const amrex::Real time, const FieldState fstate) { diff --git a/amr-wind/equation_systems/icns/icns.H b/amr-wind/equation_systems/icns/icns.H index f6e11b1fd8..ec49d4e825 100644 --- a/amr-wind/equation_systems/icns/icns.H +++ b/amr-wind/equation_systems/icns/icns.H @@ -42,7 +42,6 @@ struct ICNS : VectorTransport static constexpr bool multiply_rho = true; static constexpr bool has_diffusion = true; - // n+1/2 velocity for advection iterations static constexpr bool need_nph_state = true; }; diff --git a/amr-wind/equation_systems/icns/icns_advection.H b/amr-wind/equation_systems/icns/icns_advection.H index 77d61fabab..c7032fec89 100644 --- a/amr-wind/equation_systems/icns/icns_advection.H +++ b/amr-wind/equation_systems/icns/icns_advection.H @@ -157,12 +157,6 @@ struct AdvectionOp 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 // @@ -244,6 +238,7 @@ struct AdvectionOp // cppcheck-suppress constVariableReference auto& conv_term = fields.conv_term; const auto& dof_field = fields.field.state(fstate); + const auto& dof_nph = fields.field.state(amr_wind::FieldState::NPH); auto flux_x = repo.create_scratch_field(ICNS::ndim, 0, amr_wind::FieldLoc::XFACE); @@ -260,6 +255,8 @@ struct AdvectionOp const auto& rho_o = repo.get_field("density").state(amr_wind::FieldState::Old); + const auto& rho_nph = + repo.get_field("density").state(amr_wind::FieldState::NPH); const bool mphase_vof = repo.field_exists("vof"); @@ -280,23 +277,28 @@ struct AdvectionOp ICNS::ndim, fvm::Godunov::nghost_src); amrex::MultiFab::Copy( fq, src_term(lev), 0, 0, ICNS::ndim, fvm::Godunov::nghost_src); + // form multifab for time-correct boundary condition of variable + amrex::MultiFab q_nph( + dof_field(lev).boxArray(), dof_field(lev).DistributionMap(), + ICNS::ndim, fvm::Godunov::nghost_state); + amrex::MultiFab::Copy( + q_nph, dof_nph(lev), 0, 0, ICNS::ndim, + fvm::Godunov::nghost_state); // Calculate fluxes using momentum directly if (!mphase_vof) { - amrex::MultiFab::Multiply( - q, rho_o(lev), 0, 0, 1, fvm::Godunov::nghost_state); - amrex::MultiFab::Multiply( - q, rho_o(lev), 0, 1, 1, fvm::Godunov::nghost_state); - amrex::MultiFab::Multiply( - q, rho_o(lev), 0, 2, 1, fvm::Godunov::nghost_state); - // Source terms are at old state during calculation of advection - // terms - amrex::MultiFab::Multiply( - fq, rho_o(lev), 0, 0, 1, fvm::Godunov::nghost_src); - amrex::MultiFab::Multiply( - fq, rho_o(lev), 0, 1, 1, fvm::Godunov::nghost_src); - amrex::MultiFab::Multiply( - fq, rho_o(lev), 0, 2, 1, fvm::Godunov::nghost_src); + + for (int idim = 0; idim < dof_field.num_comp(); ++idim) { + amrex::MultiFab::Multiply( + q, rho_o(lev), 0, idim, 1, fvm::Godunov::nghost_state); + // Source terms at old state during advection calculation + amrex::MultiFab::Multiply( + fq, rho_o(lev), 0, idim, 1, fvm::Godunov::nghost_src); + + amrex::MultiFab::Multiply( + q_nph, rho_nph(lev), 0, idim, 1, + fvm::Godunov::nghost_state); + } } amrex::MFItInfo mfi_info; @@ -348,10 +350,11 @@ struct AdvectionOp const auto& divu = tmpfab.array(); HydroUtils::ComputeFluxesOnBoxFromState( bx, ICNS::ndim, mfi, q.const_array(mfi), - (*flux_x)(lev).array(mfi), (*flux_y)(lev).array(mfi), - (*flux_z)(lev).array(mfi), (*face_x)(lev).array(mfi), - (*face_y)(lev).array(mfi), (*face_z)(lev).array(mfi), - known_edge_state, u_mac(lev).const_array(mfi), + q_nph.const_array(mfi), (*flux_x)(lev).array(mfi), + (*flux_y)(lev).array(mfi), (*flux_z)(lev).array(mfi), + (*face_x)(lev).array(mfi), (*face_y)(lev).array(mfi), + (*face_z)(lev).array(mfi), 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_extrap, dof_field.bcrec(), @@ -372,8 +375,8 @@ struct AdvectionOp // Loop levels multiphase::hybrid_fluxes( repo, ICNS::ndim, iconserv, (*flux_x), (*flux_y), (*flux_z), - dof_field, src_term, rho_o, u_mac, v_mac, w_mac, - dof_field.bcrec(), dof_field.bcrec_device().data(), + dof_field, dof_nph, src_term, rho_o, rho_nph, u_mac, v_mac, + w_mac, dof_field.bcrec(), dof_field.bcrec_device().data(), rho_o.bcrec(), rho_o.bcrec_device().data(), dt, mflux_scheme, godunov_use_forces_in_trans); } diff --git a/amr-wind/equation_systems/vof/vof_momentum_flux.H b/amr-wind/equation_systems/vof/vof_momentum_flux.H index 3c66aef57a..f082f06dde 100644 --- a/amr-wind/equation_systems/vof/vof_momentum_flux.H +++ b/amr-wind/equation_systems/vof/vof_momentum_flux.H @@ -13,8 +13,10 @@ static void hybrid_fluxes( ScratchField& flux_y, ScratchField& flux_z, const Field& dof_field, + const Field& dof_nph, const Field& src_term, const Field& rho_o, + const Field& rho_nph, const Field& u_mac, const Field& v_mac, const Field& w_mac, @@ -69,24 +71,30 @@ static void hybrid_fluxes( amrex::MultiFab::Copy( fq, src_term(lev), 0, 0, ncomp, fvm::Godunov::nghost_src); - amrex::MultiFab::Multiply( - q, rho_o(lev), 0, 0, 1, fvm::Godunov::nghost_state); - amrex::MultiFab::Multiply( - q, rho_o(lev), 0, 1, 1, fvm::Godunov::nghost_state); - amrex::MultiFab::Multiply( - q, rho_o(lev), 0, 2, 1, fvm::Godunov::nghost_state); - amrex::MultiFab::Multiply( - fq, rho_o(lev), 0, 0, 1, fvm::Godunov::nghost_src); - amrex::MultiFab::Multiply( - fq, rho_o(lev), 0, 1, 1, fvm::Godunov::nghost_src); - amrex::MultiFab::Multiply( - fq, rho_o(lev), 0, 2, 1, fvm::Godunov::nghost_src); + for (int idim = 0; idim < dof_field.num_comp(); ++idim) { + amrex::MultiFab::Multiply( + q, rho_o(lev), 0, idim, 1, fvm::Godunov::nghost_state); + amrex::MultiFab::Multiply( + fq, rho_o(lev), 0, idim, 1, fvm::Godunov::nghost_src); + } amrex::MultiFab frho( src_term(lev).boxArray(), src_term(lev).DistributionMap(), 1, fvm::Godunov::nghost_src); frho.setVal(0.0); + // Set up temporary arrays for NPH at boundaries + amrex::MultiFab q_nph( + dof_field(lev).boxArray(), dof_field(lev).DistributionMap(), ncomp, + fvm::Godunov::nghost_state); + amrex::MultiFab::Copy( + q_nph, dof_nph(lev), 0, 0, ncomp, fvm::Godunov::nghost_state); + // Density at NPH is fully calculated a priori + // Multiply to get momentum at NPH + for (int idim = 0; idim < dof_field.num_comp(); ++idim) { + amrex::MultiFab::Multiply(q_nph, rho_nph(lev), 0, idim, 1, 1); + } + std::string advection_type = "Godunov"; int limiter_type = PPM::UPWIND; // default if (mflux_scheme == godunov::scheme::MINMOD) { @@ -128,14 +136,15 @@ static void hybrid_fluxes( // Compute momentum flux quantities (multiplied by face velocity) HydroUtils::ComputeFluxesOnBoxFromState( - bx, ncomp, mfi, q.const_array(mfi), Fw_x, Fw_y, Fw_z, - (*face_x)(lev).array(mfi), (*face_y)(lev).array(mfi), - (*face_z)(lev).array(mfi), 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, velbc, velbc_d, iconserv.data(), godunov_use_ppm, - use_forces_in_trans, is_velocity, fluxes_are_area_weighted, - advection_type, limiter_type, allow_inflow_on_outflow); + bx, ncomp, mfi, q.const_array(mfi), q_nph.const_array(mfi), + Fw_x, Fw_y, Fw_z, (*face_x)(lev).array(mfi), + (*face_y)(lev).array(mfi), (*face_z)(lev).array(mfi), + 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, velbc, velbc_d, + iconserv.data(), godunov_use_ppm, use_forces_in_trans, + is_velocity, fluxes_are_area_weighted, advection_type, + limiter_type, allow_inflow_on_outflow); const auto& bxg1 = amrex::grow(bx, 1); const auto& xbx = amrex::surroundingNodes(bx, 0); diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 406af2077e..ae5d43978e 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -96,6 +96,7 @@ public: bool regrid_and_update(); void pre_advance_stage1(); void pre_advance_stage2(); + void prepare_time_step(); void do_advance(const int fixed_point_iteration); void advance(const int fixed_point_iteration); void prescribe_advance(); diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp index 2bfaa16e30..29ca8075d6 100644 --- a/amr-wind/incflo.cpp +++ b/amr-wind/incflo.cpp @@ -132,6 +132,10 @@ void incflo::prepare_for_time_integration() return; } + if (m_do_initial_proj || m_initial_iterations > 0) { + m_sim.pde_manager().prepare_boundaries(); + } + if (m_do_initial_proj) { InitialProjection(); } diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index a78480f1d2..74dc949da0 100644 --- a/amr-wind/incflo_advance.cpp +++ b/amr-wind/incflo_advance.cpp @@ -30,6 +30,15 @@ void incflo::pre_advance_stage2() m_sim.helics().pre_advance_work(); } +void incflo::prepare_time_step() +{ + m_sim.pde_manager().advance_states(); + m_sim.pde_manager().prepare_boundaries(); + for (auto& pp : m_sim.physics()) { + pp->pre_predictor_work(); + } +} + /** Advance simulation state by one timestep * * Performs the following actions at a given timestep @@ -51,7 +60,7 @@ void incflo::advance(const int fixed_point_iteration) { BL_PROFILE("amr-wind::incflo::advance"); if (fixed_point_iteration == 0) { - m_sim.pde_manager().advance_states(); + prepare_time_step(); } ApplyPredictor(false, fixed_point_iteration); @@ -656,7 +665,7 @@ void incflo::prescribe_advance() { BL_PROFILE("amr-wind::incflo::prescribe_advance"); - m_sim.pde_manager().advance_states(); + prepare_time_step(); ApplyPrescribeStep(); } diff --git a/amr-wind/wind_energy/ABL.H b/amr-wind/wind_energy/ABL.H index c7ed06768c..f03fa227f6 100644 --- a/amr-wind/wind_energy/ABL.H +++ b/amr-wind/wind_energy/ABL.H @@ -80,6 +80,8 @@ public: void pre_advance_work() override; + void pre_predictor_work() override; + void post_advance_work() override; void register_forcing_term(pde::icns::ABLForcing* forcing) const diff --git a/amr-wind/wind_energy/ABL.cpp b/amr-wind/wind_energy/ABL.cpp index c1873f36ba..7c50d78d4b 100644 --- a/amr-wind/wind_energy/ABL.cpp +++ b/amr-wind/wind_energy/ABL.cpp @@ -242,6 +242,12 @@ void ABL::pre_advance_work() m_abl_mpl->pre_advance_work(); } +/** Perform tasks at the beginning of a timestep, but after pre_advance + * + * For ABL simulations, this method updates plane data to new_time (n+1) + */ +void ABL::pre_predictor_work() { m_bndry_plane->pre_predictor_work(); } + /** Perform tasks at the end of a new timestep * * For ABL simulations, this method writes all plane-averaged profiles and diff --git a/amr-wind/wind_energy/ABLBoundaryPlane.H b/amr-wind/wind_energy/ABLBoundaryPlane.H index 0a6d8b5813..692849a3a6 100644 --- a/amr-wind/wind_energy/ABLBoundaryPlane.H +++ b/amr-wind/wind_energy/ABLBoundaryPlane.H @@ -118,6 +118,8 @@ public: void pre_advance_work(); + void pre_predictor_work(); + void post_advance_work(); void initialize_data(); @@ -128,7 +130,7 @@ public: void read_header(); - void read_file(); + void read_file(const bool /* nph_target_time*/); void populate_data( const int /*lev*/, @@ -156,6 +158,12 @@ public: const int /*lev*/, const amrex::Orientation /*ori*/) const; + bool is_data_newer_than(const amrex::Real time) const + { + // Will replace with m_in_data.tinterp() + return (m_intended_interp_time - time > 1e-12); + } + io_mode mode() const { return m_io_mode; } private: @@ -187,6 +195,8 @@ private: //! Start outputting after this time amrex::Real m_out_start_time{0.0}; + amrex::Real m_intended_interp_time{0.0}; + #ifdef AMR_WIND_USE_NETCDF //! NetCDF time output counter size_t m_out_counter{0}; diff --git a/amr-wind/wind_energy/ABLBoundaryPlane.cpp b/amr-wind/wind_energy/ABLBoundaryPlane.cpp index 5175cf33d9..9f220b5e55 100644 --- a/amr-wind/wind_energy/ABLBoundaryPlane.cpp +++ b/amr-wind/wind_energy/ABLBoundaryPlane.cpp @@ -301,7 +301,7 @@ void ABLBoundaryPlane::post_init_actions() write_header(); write_file(); read_header(); - read_file(); + read_file(false); } void ABLBoundaryPlane::pre_advance_work() @@ -309,7 +309,15 @@ void ABLBoundaryPlane::pre_advance_work() if (!m_is_initialized) { return; } - read_file(); + read_file(true); +} + +void ABLBoundaryPlane::pre_predictor_work() +{ + if (!m_is_initialized) { + return; + } + read_file(false); } void ABLBoundaryPlane::post_advance_work() @@ -751,7 +759,7 @@ void ABLBoundaryPlane::read_header() } } -void ABLBoundaryPlane::read_file() +void ABLBoundaryPlane::read_file(const bool nph_target_time) { BL_PROFILE("amr-wind::ABLBoundaryPlane::read_file"); if (m_io_mode != io_mode::input) { @@ -760,10 +768,14 @@ void ABLBoundaryPlane::read_file() // populate planes and interpolate const amrex::Real time = m_time.new_time(); + const amrex::Real conditional_time = + m_time.new_time() + (nph_target_time ? 0.5 : 0.0) * + (m_time.current_time() - m_time.new_time()); AMREX_ALWAYS_ASSERT((m_in_times[0] <= time) && (time < m_in_times.back())); // return early if current data files can still be interpolated in time if ((m_in_data.tn() <= time) && (time < m_in_data.tnp1())) { + m_intended_interp_time = conditional_time; m_in_data.interpolate(time); return; } @@ -859,6 +871,7 @@ void ABLBoundaryPlane::read_file() } } + m_intended_interp_time = conditional_time; m_in_data.interpolate(time); } diff --git a/amr-wind/wind_energy/ABLFillInflow.H b/amr-wind/wind_energy/ABLFillInflow.H index 1fb440212b..a71f831ac2 100644 --- a/amr-wind/wind_energy/ABLFillInflow.H +++ b/amr-wind/wind_energy/ABLFillInflow.H @@ -41,6 +41,7 @@ public: amrex::Array& cfabs, const amrex::IntVect& nghost, const amrex::Vector& bcrec, + const amrex::Vector& /* unused */, const FieldState fstate = FieldState::New) override; //! Implementation that handles filling patches from a coarse to fine level diff --git a/amr-wind/wind_energy/ABLFillInflow.cpp b/amr-wind/wind_energy/ABLFillInflow.cpp index 503f0c620f..8d20de5722 100644 --- a/amr-wind/wind_energy/ABLFillInflow.cpp +++ b/amr-wind/wind_energy/ABLFillInflow.cpp @@ -61,10 +61,14 @@ void ABLFillInflow::fillpatch_sibling_fields( amrex::Array& cfabs, const amrex::IntVect& nghost, const amrex::Vector& bcrec, + const amrex::Vector& /* unused */, const FieldState fstate) { + // Avoid trying to read after planes have already been populated + const bool plane_data_unchanged = m_bndry_plane.is_data_newer_than(time); // For an ABL fill, we just foextrap the mac velocities - amrex::Vector lbcrec(m_field.num_comp()); + amrex::Vector fp_bcrec(m_field.num_comp()); + amrex::Vector ph_bcrec(m_field.num_comp()); const auto& ibctype = m_field.bc_type(); for (amrex::OrientationIter oit; oit != nullptr; ++oit) { auto ori = oit(); @@ -74,28 +78,38 @@ void ABLFillInflow::fillpatch_sibling_fields( for (int i = 0; i < m_field.num_comp(); ++i) { if (bct == BC::mass_inflow) { if (side == amrex::Orientation::low) { - lbcrec[i].setLo(dir, amrex::BCType::foextrap); + ph_bcrec[i].setLo(dir, amrex::BCType::foextrap); + fp_bcrec[i].setLo( + dir, plane_data_unchanged ? amrex::BCType::ext_dir + : amrex::BCType::foextrap); } else { - lbcrec[i].setHi(dir, amrex::BCType::foextrap); + ph_bcrec[i].setHi(dir, amrex::BCType::foextrap); + fp_bcrec[i].setHi( + dir, plane_data_unchanged ? amrex::BCType::ext_dir + : amrex::BCType::foextrap); } } else { if (side == amrex::Orientation::low) { - lbcrec[i].setLo(dir, bcrec[i].lo(dir)); + ph_bcrec[i].setLo(dir, bcrec[i].lo(dir)); + fp_bcrec[i].setLo(dir, bcrec[i].lo(dir)); } else { - lbcrec[i].setHi(dir, bcrec[i].hi(dir)); + ph_bcrec[i].setHi(dir, bcrec[i].hi(dir)); + fp_bcrec[i].setHi(dir, bcrec[i].hi(dir)); } } } } FieldFillPatchOps::fillpatch_sibling_fields( - lev, time, mfabs, ffabs, cfabs, nghost, lbcrec, fstate); + lev, time, mfabs, ffabs, cfabs, nghost, fp_bcrec, ph_bcrec, fstate); - for (int i = 0; i < static_cast(mfabs.size()); i++) { - // use new_time to populate boundary data instead of half-time - // to avoid interpolating from precursor data - m_bndry_plane.populate_data( - lev, m_time.new_time(), m_field, *mfabs[i], 0, i); + if (!plane_data_unchanged) { + for (int i = 0; i < static_cast(mfabs.size()); i++) { + // use new_time to populate boundary data instead of half-time + // to avoid interpolating from precursor data + m_bndry_plane.populate_data( + lev, m_time.new_time(), m_field, *mfabs[i], 0, i); + } } } diff --git a/amr-wind/wind_energy/ABLFillMPL.H b/amr-wind/wind_energy/ABLFillMPL.H index 2d5414b1a9..39cbf1fe68 100644 --- a/amr-wind/wind_energy/ABLFillMPL.H +++ b/amr-wind/wind_energy/ABLFillMPL.H @@ -41,6 +41,7 @@ public: amrex::Array& cfabs, const amrex::IntVect& nghost, const amrex::Vector& bcrec, + const amrex::Vector& /* unused */, const FieldState fstate = FieldState::New) override; //! Implementation that handles filling patches from a coarse to fine level diff --git a/amr-wind/wind_energy/ABLFillMPL.cpp b/amr-wind/wind_energy/ABLFillMPL.cpp index bd620deb8d..29cce18ac9 100644 --- a/amr-wind/wind_energy/ABLFillMPL.cpp +++ b/amr-wind/wind_energy/ABLFillMPL.cpp @@ -73,6 +73,7 @@ void ABLFillMPL::fillpatch_sibling_fields( amrex::Array& cfabs, const amrex::IntVect& nghost, const amrex::Vector& bcrec, + const amrex::Vector& /* unused */, const FieldState fstate) { if (m_field.base_name() == "velocity") { @@ -102,7 +103,7 @@ void ABLFillMPL::fillpatch_sibling_fields( } FieldFillPatchOps::fillpatch_sibling_fields( - lev, time, mfabs, ffabs, cfabs, nghost, lbcrec, fstate); + lev, time, mfabs, ffabs, cfabs, nghost, lbcrec, lbcrec, fstate); for (int i = 0; i < static_cast(mfabs.size()); i++) { m_abl_mpl.set_velocity(lev, time, m_field, *mfabs[i], 0, i);