Fixed point iterations within a physical time step for the advection term

amr-wind/equation_systems/AdvOp_Godunov.H

@@ -168,6 +168,11 @@ struct AdvectionOp<
                         (godunov_scheme == godunov::scheme::WENO_JS)) {
                         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;
+
                     HydroUtils::ComputeFluxesOnBoxFromState(
                         bx, PDE::ndim, mfi,
                         (PDE::multiply_rho ? rhotrac : tra_arr),
@@ -177,7 +182,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_extrap,
                         dof_field.bcrec(), dof_field.bcrec_device().data(),
                         iconserv.data(), godunov_use_ppm,
                         godunov_use_forces_in_trans, is_velocity,

amr-wind/equation_systems/icns/icns.H

@@ -42,8 +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 = false;
+    // n+1/2 velocity for advection iterations
+    static constexpr bool need_nph_state = true;
 };
 
 } // namespace amr_wind::pde

amr-wind/equation_systems/icns/icns_advection.H

@@ -190,11 +190,15 @@ struct AdvectionOp<ICNS, fvm::Godunov>
                 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;
             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_extrap, godunov_use_ppm,
                     godunov_use_forces_in_trans, premac_advection_type,
                     limiter_type, m_allow_inflow_on_outflow);
             }
@@ -329,6 +333,10 @@ struct AdvectionOp<ICNS, fvm::Godunov>
                     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
@@ -346,7 +354,7 @@ struct AdvectionOp<ICNS, fvm::Godunov>
                         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_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,

amr-wind/incflo.H

@@ -96,8 +96,8 @@ public:
     bool regrid_and_update();
     void pre_advance_stage1();
     void pre_advance_stage2();
-    void do_advance();
-    void advance();
+    void do_advance(const int fixed_point_iteration);
+    void advance(const int fixed_point_iteration);
     void prescribe_advance();
     void post_advance_work();
 
@@ -133,7 +133,9 @@ public:
     void ComputeDt(bool explicit_diffusion);
     void ComputePrescribeDt();
 
-    void ApplyPredictor(bool incremental_projection = false);
+    void ApplyPredictor(
+        const bool incremental_projection = false,
+        const int fixed_point_iteration = 0);
     void ApplyCorrector();
     void ApplyPrescribeStep();
 
@@ -172,6 +174,9 @@ private:
     // Prescribe advection velocity
     bool m_prescribe_vel = false;
 
+    // Fixed point iterations every timestep
+    int m_fixed_point_iterations{1};
+
     //! number of cells on all levels including covered cells
     amrex::Long m_cell_count{-1};
 

amr-wind/incflo.cpp

@@ -276,7 +276,10 @@ void incflo::Evolve()
 
         amrex::Real time1 = amrex::ParallelDescriptor::second();
         // Advance to time t + dt
-        do_advance();
+        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();
@@ -309,15 +312,17 @@ void incflo::Evolve()
     }
 }
 
-void incflo::do_advance()
+void incflo::do_advance(const int fixed_point_iteration)
 {
     if (m_sim.has_overset()) {
         m_ovst_ops.pre_advance_work();
     }
-    if (m_prescribe_vel) {
+    if (m_prescribe_vel && fixed_point_iteration == 0) {
         prescribe_advance();
     } else {
-        advance();
+        amrex::Print() << "Fixed point iteration " << fixed_point_iteration
+                       << std::endl;
+        advance(fixed_point_iteration);
     }
     if (m_sim.has_overset()) {
         m_ovst_ops.post_advance_work();
@@ -372,6 +377,11 @@ void incflo::init_physics_and_pde()
         if (activate_overset) {
             m_sim.activate_overset();
         }
+
+        pp.query("fixed_point_iterations", m_fixed_point_iterations);
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+            m_fixed_point_iterations > 0,
+            "The number of fixed point iterations cannot be less than 1");
     }
 
     auto& pde_mgr = m_sim.pde_manager();

amr-wind/incflo_advance.cpp

@@ -47,15 +47,20 @@ void incflo::pre_advance_stage2()
  *
  * \callgraph
  */
-void incflo::advance()
+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();
+    }
 
-    m_sim.pde_manager().advance_states();
-
-    ApplyPredictor();
+    ApplyPredictor(false, fixed_point_iteration);
 
     if (!m_use_godunov) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+            fixed_point_iteration == 0,
+            "Fixed point iterations are not supported for MOL");
+
         ApplyCorrector();
     }
 }
@@ -165,10 +170,10 @@ void incflo::advance()
  *  <li> \ref incflo::ApplyProjection "Apply projection"
  *  </ol>
  */
-void incflo::ApplyPredictor(bool incremental_projection)
+void incflo::ApplyPredictor(
+    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
     Real new_time = m_time.new_time();
 
@@ -190,6 +195,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 (fixed_point_iteration > 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
     // *************************************************************************************
@@ -268,7 +283,10 @@ void incflo::ApplyPredictor(bool incremental_projection)
     }
 
     // Extrapolate and apply MAC projection for advection velocities
-    icns().pre_advection_actions(amr_wind::FieldState::Old);
+    const auto fstate_preadv = (fixed_point_iteration == 0)
+                                   ? amr_wind::FieldState::Old
+                                   : amr_wind::FieldState::NPH;
+    icns().pre_advection_actions(fstate_preadv);
 
     // For scalars only first
     // *************************************************************************************
@@ -336,11 +354,14 @@ void incflo::ApplyPredictor(bool incremental_projection)
     }
 
     // With scalars computed, compute advection of momentum
-    icns().compute_advection_term(amr_wind::FieldState::Old);
+    const auto fstate = (fixed_point_iteration == 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
+    // 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);
@@ -397,8 +418,25 @@ void incflo::ApplyPredictor(bool incremental_projection)
     if (m_verbose > 2) {
         PrintMaxVelLocations("after nodal projection");
     }
+    // ScratchField to store the old np1
+    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);
+
+        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:
 //

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;

amr-wind/utilities/console_io.H

@@ -3,6 +3,7 @@
 
 #include <iostream>
 #include "AMReX_MLMG.H"
+#include "amr-wind/CFDSim.H"
 
 namespace amr_wind::io {
 
@@ -20,6 +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(
+    const CFDSim& sim, ScratchField& vel_diff, const ScratchField& vel_star);
+
 } // namespace amr_wind::io
 
 #endif /* CONSOLE_IO_H */

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,60 @@ void print_tpls(std::ostream& out)
     }
 }
 
+void print_nonlinear_residual(
+    const CFDSim& sim, ScratchField& vel_diff, const ScratchField& vel_star)
+{
+    const int nlevels = sim.repo().num_active_levels();
+    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");
+
+    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& 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();
+        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 {
+                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);
+                }
+            });
+    }
+
+    amrex::Array<amrex::Real, AMREX_SPACEDIM> rms_vel = {0.0};
+
+    for (int lev = 0; lev < nlevels; ++lev) {
+        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
+            rms_vel[idim] += vel_diff(lev).norm2(idim);
+        }
+    }
+
+    amrex::Print() << "Norm of change over fixed point iterations in u:  "
+                   << rms_vel[0] << ", v: " << rms_vel[1]
+                   << ", w: " << rms_vel[2] << std::endl;
+}
+
 } // namespace amr_wind::io