Vof fillpatch fixes (#1425)

* for postsolve, fillpatch should be at new_time. remove unnecessary density operation here, too * this version should pass the assert, but it is not fully accurate temporally * remove fillpatches within set_density calls. Assume input variable has already been fillpatched in some way * new reg test for using boundary planes with waves
Exawind · Dec 31, 2024 · e13a6a4 · e13a6a4
1 parent d905dce
commit e13a6a4
Show file tree

Hide file tree

Showing 9 changed files with 133 additions and 16 deletions.
diff --git a/amr-wind/equation_systems/vof/SplitAdvection.H b/amr-wind/equation_systems/vof/SplitAdvection.H
@@ -18,6 +18,7 @@ void split_advection_step(
     Field const& w_mac,
     amrex::GpuArray<BC, AMREX_SPACEDIM * 2> BCs,
     amrex::Vector<amrex::Geometry> geom,
+    amrex::Real time,
     amrex::Real dt,
     bool rm_debris);
 

diff --git a/amr-wind/equation_systems/vof/SplitAdvection.cpp b/amr-wind/equation_systems/vof/SplitAdvection.cpp
@@ -21,6 +21,7 @@ void multiphase::split_advection_step(
     Field const& w_mac,
     amrex::GpuArray<BC, AMREX_SPACEDIM * 2> BCs,
     amrex::Vector<amrex::Geometry> geom,
+    amrex::Real time,
     amrex::Real dt,
     bool rm_debris)
 {
@@ -99,7 +100,7 @@ void multiphase::split_advection_step(
     }
 
     // Average down vof and communicate
-    dof_field.fillpatch(0.0);
+    dof_field.fillpatch(time);
 }
 
 void multiphase::split_compute_fluxes(

diff --git a/amr-wind/equation_systems/vof/vof_advection.H b/amr-wind/equation_systems/vof/vof_advection.H
@@ -5,6 +5,7 @@
 #include "amr-wind/equation_systems/vof/vof.H"
 #include "amr-wind/equation_systems/vof/vof_hybridsolver_ops.H"
 #include "amr-wind/equation_systems/vof/SplitAdvection.H"
+#include "amr-wind/core/SimTime.H"
 
 namespace amr_wind::pde {
 
@@ -15,13 +16,14 @@ template <>
 struct AdvectionOp<VOF, fvm::Godunov>
 {
     AdvectionOp(
-        CFDSim& /* sim */,
+        CFDSim& sim,
         PDEFields& fields_in,
         bool /*unused*/,
         bool /*unused*/,
         bool /*unused*/,
         bool /*unused*/)
-        : fields(fields_in)
+        : m_time(sim.time())
+        , fields(fields_in)
         , u_mac(fields_in.repo.get_field("u_mac"))
         , v_mac(fields_in.repo.get_field("v_mac"))
         , w_mac(fields_in.repo.get_field("w_mac"))
@@ -105,15 +107,20 @@ struct AdvectionOp<VOF, fvm::Godunov>
         // Split advection step 1, with cmask calculation
         multiphase::split_advection_step(
             isweep, 0, nlevels, dof_field, fluxes, (*fluxC), advas, u_mac,
-            v_mac, w_mac, dof_field.bc_type(), geom, dt, m_rm_debris);
+            v_mac, w_mac, dof_field.bc_type(), geom, m_time.new_time(), dt,
+            m_rm_debris);
+        // (copy old boundaries to working state)
         // Split advection step 2
         multiphase::split_advection_step(
             isweep, 1, nlevels, dof_field, fluxes, (*fluxC), advas, u_mac,
-            v_mac, w_mac, dof_field.bc_type(), geom, dt, m_rm_debris);
+            v_mac, w_mac, dof_field.bc_type(), geom, m_time.new_time(), dt,
+            m_rm_debris);
+        // (copy old boundaries to working state)
         // Split advection step 3
         multiphase::split_advection_step(
             isweep, 2, nlevels, dof_field, fluxes, (*fluxC), advas, u_mac,
-            v_mac, w_mac, dof_field.bc_type(), geom, dt, m_rm_debris);
+            v_mac, w_mac, dof_field.bc_type(), geom, m_time.new_time(), dt,
+            m_rm_debris);
 
         // Replace masked cells using overset
         if (repo.int_field_exists("iblank_cell") && m_replace_mask) {
@@ -130,6 +137,7 @@ struct AdvectionOp<VOF, fvm::Godunov>
         }
     }
 
+    const SimTime& m_time;
     PDEFields& fields;
     Field& u_mac;
     Field& v_mac;

diff --git a/amr-wind/equation_systems/vof/vof_ops.H b/amr-wind/equation_systems/vof/vof_ops.H
@@ -70,13 +70,9 @@ struct PostSolveOp<VOF>
 
     void operator()(const amrex::Real /*unused*/)
     {
-        m_fields.field.fillpatch(m_sim.time().current_time());
+        m_fields.field.fillpatch(m_sim.time().new_time());
         auto& multiphase = m_sim.physics_manager().get<MultiPhase>();
         multiphase.set_density_via_vof();
-        // Do this for old state too if modified by vof sharpening
-        if (m_sim.repo().int_field_exists("iblank_cell")) {
-            multiphase.set_density_via_vof(amr_wind::FieldState::Old);
-        }
 
         // Calculate n+1/2 density
         multiphase.set_nph_density();

diff --git a/amr-wind/physics/multiphase/MultiPhase.cpp b/amr-wind/physics/multiphase/MultiPhase.cpp
@@ -305,7 +305,7 @@ void MultiPhase::set_density_via_levelset()
         auto& levelset = (*m_levelset)(lev);
 
         for (amrex::MFIter mfi(density); mfi.isValid(); ++mfi) {
-            const auto& vbx = mfi.validbox();
+            const auto& vbx = mfi.growntilebox();
             const auto& dx = geom[lev].CellSizeArray();
 
             const amrex::Array4<amrex::Real>& phi = levelset.array(mfi);
@@ -331,7 +331,6 @@ void MultiPhase::set_density_via_levelset()
                 });
         }
     }
-    m_density.fillpatch(m_sim.time().current_time());
 }
 
 void MultiPhase::set_density_via_vof(amr_wind::FieldState fstate)
@@ -343,7 +342,7 @@ void MultiPhase::set_density_via_vof(amr_wind::FieldState fstate)
         auto& vof = (*m_vof).state(fstate)(lev);
 
         for (amrex::MFIter mfi(density); mfi.isValid(); ++mfi) {
-            const auto& vbx = mfi.validbox();
+            const auto& vbx = mfi.growntilebox();
             const amrex::Array4<amrex::Real>& F = vof.array(mfi);
             const amrex::Array4<amrex::Real>& rho = density.array(mfi);
             const amrex::Real captured_rho1 = m_rho1;
@@ -355,7 +354,6 @@ void MultiPhase::set_density_via_vof(amr_wind::FieldState fstate)
                 });
         }
     }
-    m_density.fillpatch(m_sim.time().current_time());
 }
 
 void MultiPhase::set_nph_density()

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
@@ -375,6 +375,7 @@ if(AMR_WIND_TEST_WITH_FCOMPARE)
   add_test_red(abl_bndry_input_amr_native_xhi abl_bndry_output_native)
   add_test_red(abl_bndry_input_amr_native_mlbc abl_bndry_output_amr_native)
   add_test_red(abl_godunov_forcetimetable abl_godunov_timetable)
+  add_test_red(abl_multiphase_laminar_input abl_multiphase_laminar)
 endif()
 
 if(AMR_WIND_ENABLE_NETCDF)

diff --git a/test/test_files/abl_multiphase_laminar/abl_multiphase_laminar.inp b/test/test_files/abl_multiphase_laminar/abl_multiphase_laminar.inp
@@ -12,7 +12,7 @@ time.cfl              =   0.95       # CFL factor
 #            INPUT AND OUTPUT           #
 #.......................................#
 time.plot_interval            =  10       # Steps between plot files
-time.checkpoint_interval      =  -1         # Steps between checkpoint files
+time.checkpoint_interval      =  5        # Steps between checkpoint files
 io.outputs = velocity_mueff temperature_mueff
 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
 #               PHYSICS                 #
@@ -64,6 +64,13 @@ OceanWaves.Wave1.water_depth=400
 OceanWaves.Wave1.relax_zone_gen_length=40.0
 OceanWaves.Wave1.relax_zone_out_length=40.0
 
+ABL.bndry_file = "bndry_files"
+ABL.bndry_io_mode = 0
+ABL.bndry_planes = xlo
+ABL.bndry_output_start_time = 0.2
+ABL.bndry_var_names = velocity temperature vof density
+ABL.bndry_output_format = native
+
 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
 #        ADAPTIVE MESH REFINEMENT       #
 #.......................................#

diff --git a/test/test_files/abl_multiphase_laminar_input/abl_multiphase_laminar_input.inp b/test/test_files/abl_multiphase_laminar_input/abl_multiphase_laminar_input.inp
@@ -0,0 +1,100 @@
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            SIMULATION STOP            #
+#.......................................#
+time.stop_time               = 28800.0     # Max (simulated) time to evolve
+time.max_step                =  10          # Max number of time steps
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#         TIME STEP COMPUTATION         #
+#.......................................#
+time.fixed_dt         =   0.05        # Use this constant dt if > 0
+time.cfl              =   0.95       # CFL factor
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#            INPUT AND OUTPUT           #
+#.......................................#
+io.restart_file = "../abl_multiphase_laminar/chk00005"
+time.plot_interval            =  10       # Steps between plot files
+time.checkpoint_interval      =  -1         # Steps between checkpoint files
+io.outputs = velocity_mueff temperature_mueff
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#               PHYSICS                 #
+#.......................................#
+incflo.gravity        =  0.0  0.0 -9.81  # Gravitational force (3D)
+incflo.use_godunov = 1
+incflo.diffusion_type = 1
+incflo.godunov_type = "weno_z"
+incflo.mflux_type = "minmod"
+transport.model = TwoPhaseTransport
+transport.viscosity_fluid1 = 1e-3
+transport.viscosity_fluid2 = 1e-5
+transport.laminar_prandtl_fluid1 = 7.2
+transport.laminar_prandtl_fluid2 = 0.7
+transport.turbulent_prandtl = 0.3333
+turbulence.model = Laminar
+incflo.physics = MultiPhase ABL OceanWaves
+ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing GravityForcing
+ICNS.use_perturb_pressure = true
+MultiPhase.density_fluid1=1020.
+MultiPhase.density_fluid2=1.3223
+BoussinesqBuoyancy.reference_temperature = 263.5
+CoriolisForcing.east_vector = 1.0 0.0 0.0
+CoriolisForcing.north_vector = 0.0 1.0 0.0
+CoriolisForcing.latitude = 90.0
+CoriolisForcing.rotational_time_period = 90405.5439881955
+GeostrophicForcing.geostrophic_wind = 8.0 0.0 0.0
+GeostrophicForcing.wind_forcing_off_height = 3.0
+GeostrophicForcing.wind_forcing_ramp_height = 7.0
+incflo.velocity = 8.0 0.0 0.0
+ABL.reference_temperature = 263.5
+ABL.temperature_heights = 0.0 100 400.0
+ABL.temperature_values = 265.0 265.0 268.0
+ABL.perturb_temperature = false
+ABL.cutoff_height = 50.0
+ABL.perturb_velocity = true
+ABL.perturb_ref_height = 50.0
+ABL.Uperiods = 4.0
+ABL.Vperiods = 4.0
+ABL.deltaU = 1e-5
+ABL.deltaV = 1e-5
+ABL.normal_direction = 2
+# frequency is rapid for the sake of reg test
+ABL.stats_output_frequency = 5
+OceanWaves.label = Wave1
+OceanWaves.Wave1.type = LinearWaves
+OceanWaves.Wave1.wave_height=2.0
+OceanWaves.Wave1.wave_length=40.0
+OceanWaves.Wave1.water_depth=400
+OceanWaves.Wave1.relax_zone_gen_length=40.0
+OceanWaves.Wave1.relax_zone_out_length=40.0
+
+ABL.bndry_file = "../abl_multiphase_laminar/bndry_files"
+ABL.bndry_io_mode = 1
+ABL.bndry_planes = xlo
+ABL.bndry_output_start_time = 0.2
+ABL.bndry_var_names = velocity temperature vof density
+ABL.bndry_output_format = native
+
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#        ADAPTIVE MESH REFINEMENT       #
+#.......................................#
+amr.n_cell              = 48 48 96 # Grid cells at coarsest AMRlevel
+amr.max_level           = 2           # Max AMR level in hierarchy
+tagging.labels = static
+tagging.static.type = CartBoxRefinement
+tagging.static.static_refinement_def = static_box.txt
+#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
+#              GEOMETRY                 #
+#.......................................#
+geometry.prob_lo        =   0.       0.  -400.  # Lo corner coordinates
+geometry.prob_hi        =   400.   400.   400.  # Hi corner coordinates
+geometry.is_periodic    =   0   1   0   # Periodicity x y z (0/1)
+# Boundary conditions
+xlo.type = "mass_inflow"
+xlo.density = 1.0
+xlo.vof = 0.0
+xlo.temperature = 0.0
+xhi.type = "pressure_outflow"
+zlo.type =   "slip_wall"
+zhi.type =   "slip_wall"
+zhi.temperature_type = "fixed_gradient"
+zhi.temperature = 0.01
+incflo.verbose          =   0
diff --git a/test/test_files/abl_multiphase_laminar_input/static_box.txt b/test/test_files/abl_multiphase_laminar_input/static_box.txt
@@ -0,0 +1,5 @@
+2
+1
+0.0 0.0 -20. 400. 400. 20.
+1
+0.0 0.0 -8. 400. 400. 8.
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,5 @@
		2
		1
		0.0 0.0 -20. 400. 400. 20.
		1
		0.0 0.0 -8. 400. 400. 8.