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

@@ -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());
     }
 }
 

amr-wind/core/FieldFillPatchOps.H

@@ -60,7 +60,8 @@ public:
         amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& ffabs,
         amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& cfabs,
         const amrex::IntVect& nghost,
-        const amrex::Vector<amrex::BCRec>& bcrec,
+        const amrex::Vector<amrex::BCRec>& fillpatch_bcrec,
+        const amrex::Vector<amrex::BCRec>& 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<amrex::MultiFab*, AMREX_SPACEDIM>& /*cfabs*/,
         const amrex::IntVect& /*nghost*/,
         const amrex::Vector<amrex::BCRec>& /*bcrec*/,
+        const amrex::Vector<amrex::BCRec>& /*bcrec*/,
         const FieldState /*fstate*/) override
     {
         for (const auto& mfab : mfabs) {
@@ -272,34 +274,35 @@ public:
         amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& ffabs,
         amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>& cfabs,
         const amrex::IntVect& nghost,
-        const amrex::Vector<amrex::BCRec>& bcrec,
+        const amrex::Vector<amrex::BCRec>& fillpatch_bcrec,
+        const amrex::Vector<amrex::BCRec>& physbc_bcrec,
         const FieldState /*fstate = FieldState::New*/) override
     {
 
         if (lev == 0) {
             for (int i = 0; i < static_cast<int>(mfabs.size()); i++) {
                 amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> 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);
             }
         } else {
             AMREX_D_TERM(
                 amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> cphysbcx(
-                    m_mesh.Geom(lev - 1), bcrec, bc_functor_face(0));
+                    m_mesh.Geom(lev - 1), physbc_bcrec, bc_functor_face(0));
                 , amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> cphysbcy(
-                      m_mesh.Geom(lev - 1), bcrec, bc_functor_face(1));
+                      m_mesh.Geom(lev - 1), physbc_bcrec, bc_functor_face(1));
                 , amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> 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<amrex::GpuBndryFuncFab<Functor>> fphysbcx(
-                    m_mesh.Geom(lev), bcrec, bc_functor_face(0));
+                    m_mesh.Geom(lev), physbc_bcrec, bc_functor_face(0));
                 , amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> fphysbcy(
-                      m_mesh.Geom(lev), bcrec, bc_functor_face(1));
+                      m_mesh.Geom(lev), physbc_bcrec, bc_functor_face(1));
                 , amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>> fphysbcz(
-                      m_mesh.Geom(lev), bcrec, bc_functor_face(2)););
+                      m_mesh.Geom(lev), physbc_bcrec, bc_functor_face(2)););
 
             amrex::Array<
                 amrex::PhysBCFunct<amrex::GpuBndryFuncFab<Functor>>,
@@ -313,7 +316,8 @@ public:
 
             amrex::Array<int, AMREX_SPACEDIM> idx = {AMREX_D_DECL(0, 1, 2)};
             const amrex::Array<amrex::Vector<amrex::BCRec>, 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,

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;
 

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,22 +123,30 @@ 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<amrex::Real> rhotrac;
+                amrex::FArrayBox rhotrac_nph_fab;
+                amrex::Array4<amrex::Real> 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,
                         [=] AMREX_GPU_DEVICE(
                             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),

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);

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<Field*, AMREX_SPACEDIM> 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)
 {

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;
 };
 

amr-wind/equation_systems/icns/icns_advection.H

@@ -157,12 +157,6 @@ struct AdvectionOp<ICNS, fvm::Godunov>
         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<ICNS, fvm::Godunov>
         // 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<ICNS, fvm::Godunov>
 
         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, fvm::Godunov>
                 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<ICNS, fvm::Godunov>
                     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<ICNS, fvm::Godunov>
             // 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);
         }