diff --git a/amr-wind/ocean_waves/relaxation_zones/RelaxationZones.H b/amr-wind/ocean_waves/relaxation_zones/RelaxationZones.H index 9f8962a7bf..786c41dd3f 100644 --- a/amr-wind/ocean_waves/relaxation_zones/RelaxationZones.H +++ b/amr-wind/ocean_waves/relaxation_zones/RelaxationZones.H @@ -34,7 +34,7 @@ struct RelaxZonesBaseData bool has_beach{true}; bool has_outprofile{false}; - amrex::Real ramp_period{2.0}; + amrex::Real ramp_period; bool regrid_occurred{false}; }; diff --git a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp index b872d44754..7deed62055 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp @@ -39,7 +39,7 @@ void read_inputs( wdata.has_ramp = pp.contains("timeramp_period"); if (wdata.has_ramp) { - pp.query("timeramp_period", wdata.ramp_period); + pp.get("timeramp_period", wdata.ramp_period); } amrex::ParmParse pp_multiphase("MultiPhase"); @@ -59,6 +59,7 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata) const auto& mphase = sim.physics_manager().get(); const amrex::Real rho1 = mphase.rho1(); const amrex::Real rho2 = mphase.rho2(); + constexpr amrex::Real vof_tiny = 1e-12; for (int lev = 0; lev < nlevels; ++lev) { auto& ls = m_ow_levelset(lev); @@ -119,40 +120,47 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata) if (x <= problo[0] + gen_length) { const amrex::Real Gamma = utils::gamma_generate(x - problo[0], gen_length); - const amrex::Real vf = - (1. - Gamma) * target_volfrac(i, j, k) * rampf + + // Get bounded new vof, incorporate with increment + amrex::Real new_vof = + (1. - Gamma) * target_volfrac(i, j, k) + Gamma * volfrac(i, j, k); - volfrac(i, j, k) = (vf > 1. - 1.e-10) ? 1.0 : vf; - // Force liquid velocity, update according to mom. + new_vof = (new_vof > 1. - vof_tiny) + ? 1.0 + : (new_vof < vof_tiny ? 0.0 : new_vof); + const amrex::Real dvf = new_vof - volfrac(i, j, k); + volfrac(i, j, k) += rampf * dvf; + // Force liquid velocity only if target vof present + const amrex::Real fvel_liq = + (target_volfrac(i, j, k) > vof_tiny) ? 1.0 : 0.0; amrex::Real rho_ = rho1 * volfrac(i, j, k) + rho2 * (1.0 - volfrac(i, j, k)); - vel(i, j, k, 0) = - (rho1 * volfrac(i, j, k) * - (rampf * (1. - Gamma) * - target_vel(i, j, k, 0) + - Gamma * vel(i, j, k, 0)) + - rho2 * (1. - volfrac(i, j, k)) * vel(i, j, k, 0)) / - rho_; - vel(i, j, k, 1) = - (rho1 * volfrac(i, j, k) * - (rampf * (1. - Gamma) * - target_vel(i, j, k, 1) + - Gamma * vel(i, j, k, 1)) + - rho2 * (1. - volfrac(i, j, k)) * vel(i, j, k, 1)) / - rho_; - vel(i, j, k, 2) = - (rho1 * volfrac(i, j, k) * - (rampf * (1. - Gamma) * - target_vel(i, j, k, 2) + - Gamma * vel(i, j, k, 2)) + - rho2 * (1. - volfrac(i, j, k)) * vel(i, j, k, 2)) / - rho_; + for (int n = 0; n < vel.ncomp; ++n) { + // Get updated liquid velocity + amrex::Real vel_liq = vel(i, j, k, n); + const amrex::Real dvel_liq = + ((1. - Gamma) * target_vel(i, j, k, n) + + Gamma * vel_liq) - + vel_liq; + vel_liq += rampf * fvel_liq * dvel_liq; + // If liquid was added, that liquid has target_vel + amrex::Real integrated_vel_liq = + volfrac(i, j, k) * vel_liq; + integrated_vel_liq += + rampf * fvel_liq * amrex::max(0.0, dvf) * + (target_vel(i, j, k, n) - vel(i, j, k, n)); + // Update overall velocity using momentum + vel(i, j, k, n) = (rho1 * integrated_vel_liq + + rho2 * (1. - volfrac(i, j, k)) * + vel(i, j, k, n)) / + rho_; + } } - // Numerical beach (sponge layer) + // Outlet region if (x + beach_length >= probhi[0]) { const amrex::Real Gamma = utils::gamma_absorb( x - (probhi[0] - beach_length), beach_length, beach_length_factor); + // Numerical beach (sponge layer) if (has_beach) { volfrac(i, j, k) = (1.0 - Gamma) * @@ -162,45 +170,48 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata) amrex::Real rho_ = rho1 * volfrac(i, j, k) + rho2 * (1.0 - volfrac(i, j, k)); // Target solution in liquid is vel = 0 - vel(i, j, k, 0) = (rho1 * volfrac(i, j, k) * Gamma + - rho2 * (1. - volfrac(i, j, k))) * - vel(i, j, k, 0) / rho_; - vel(i, j, k, 1) = (rho1 * volfrac(i, j, k) * Gamma + - rho2 * (1. - volfrac(i, j, k))) * - vel(i, j, k, 1) / rho_; - vel(i, j, k, 2) = (rho1 * volfrac(i, j, k) * Gamma + - rho2 * (1. - volfrac(i, j, k))) * - vel(i, j, k, 2) / rho_; + for (int n = 0; n < vel.ncomp; ++n) { + vel(i, j, k, n) = + (rho1 * volfrac(i, j, k) * Gamma + + rho2 * (1. - volfrac(i, j, k))) * + vel(i, j, k, n) / rho_; + } } + // Forcing to wave profile instead if (has_outprofile) { - const amrex::Real vf = - (1. - Gamma) * target_volfrac(i, j, k) * rampf + + // Same steps as in wave generation region + amrex::Real new_vof = + (1. - Gamma) * target_volfrac(i, j, k) + Gamma * volfrac(i, j, k); - volfrac(i, j, k) = (vf > 1. - 1.e-10) ? 1.0 : vf; - // Force liquid velocity, update according to mom. + new_vof = + (new_vof > 1. - vof_tiny) + ? 1.0 + : (new_vof < vof_tiny ? 0.0 : new_vof); + const amrex::Real dvf = new_vof - volfrac(i, j, k); + volfrac(i, j, k) += dvf; + const amrex::Real fvel_liq = + (target_volfrac(i, j, k) > vof_tiny) ? 1.0 + : 0.0; amrex::Real rho_ = rho1 * volfrac(i, j, k) + rho2 * (1.0 - volfrac(i, j, k)); - vel(i, j, k, 0) = (rho1 * volfrac(i, j, k) * - (rampf * (1. - Gamma) * - target_vel(i, j, k, 0) + - Gamma * vel(i, j, k, 0)) + - rho2 * (1. - volfrac(i, j, k)) * - vel(i, j, k, 0)) / - rho_; - vel(i, j, k, 1) = (rho1 * volfrac(i, j, k) * - (rampf * (1. - Gamma) * - target_vel(i, j, k, 1) + - Gamma * vel(i, j, k, 1)) + - rho2 * (1. - volfrac(i, j, k)) * - vel(i, j, k, 1)) / - rho_; - vel(i, j, k, 2) = (rho1 * volfrac(i, j, k) * - (rampf * (1. - Gamma) * - target_vel(i, j, k, 2) + - Gamma * vel(i, j, k, 2)) + - rho2 * (1. - volfrac(i, j, k)) * - vel(i, j, k, 2)) / - rho_; + for (int n = 0; n < vel.ncomp; ++n) { + amrex::Real vel_liq = vel(i, j, k, n); + const amrex::Real dvel_liq = + ((1. - Gamma) * target_vel(i, j, k, n) + + Gamma * vel_liq) - + vel_liq; + vel_liq += rampf * fvel_liq * dvel_liq; + amrex::Real integrated_vel_liq = + volfrac(i, j, k) * vel_liq; + integrated_vel_liq += + rampf * fvel_liq * amrex::max(0.0, dvf) * + (target_vel(i, j, k, n) - vel(i, j, k, n)); + vel(i, j, k, n) = + (rho1 * integrated_vel_liq + + rho2 * (1. - volfrac(i, j, k)) * + vel(i, j, k, n)) / + rho_; + } } } diff --git a/docs/sphinx/user/inputs_ocean_waves.rst b/docs/sphinx/user/inputs_ocean_waves.rst index fbfd229d9c..76611f978d 100644 --- a/docs/sphinx/user/inputs_ocean_waves.rst +++ b/docs/sphinx/user/inputs_ocean_waves.rst @@ -76,7 +76,7 @@ This section is for setting up wave forcing and relaxation zones. .. input_param:: OceanWaves.label.timeramp_period - **type:** Real, optional, default = 2.0 + **type:** Real, optional An initial ramp-up period for the wave forcing. Without specifying a period, the wave forcing will begin at full strength. diff --git a/test/test_files/ow_linear/ow_linear.inp b/test/test_files/ow_linear/ow_linear.inp index 28d1e782c8..5b8461f781 100644 --- a/test/test_files/ow_linear/ow_linear.inp +++ b/test/test_files/ow_linear/ow_linear.inp @@ -20,9 +20,6 @@ time.checkpoint_interval = -1 # Steps between checkpoint files #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # PHYSICS # #.......................................# -incflo.use_godunov = 1 -incflo.diffusion_type = 2 -incflo.godunov_type="weno_z" transport.model = TwoPhaseTransport transport.viscosity_fluid1=0.0 transport.viscosity_fluid2=0.0 @@ -39,8 +36,6 @@ OceanWaves.Wave1.relax_zone_gen_length=4.0 OceanWaves.Wave1.numerical_beach_length=8.0 MultiPhase.density_fluid1=1000. MultiPhase.density_fluid2=1. -MultiPhase.interface_smoothing=0 -MultiPhase.interface_smoothing_frequency=1 ICNS.source_terms = GravityForcing ICNS.use_perturb_pressure = true ICNS.reconstruct_true_pressure = true @@ -57,15 +52,10 @@ geometry.prob_lo = 0.0 0. -1 # Lo corner coordinates geometry.prob_hi = 30. 1. 1 # Hi corner coordinates geometry.is_periodic = 0 1 0 # Periodicity x y z (0/1) -xlo.type = "pressure_inflow" +xlo.type = "wave_generation" xhi.type = "pressure_outflow" -xlo.vof_type = "zero_gradient" -xhi.vof_type = "zero_gradient" zlo.type = "slip_wall" -zhi.type = "pressure_outflow" -zlo.vof_type = "zero_gradient" -zhi.vof_type = "zero_gradient" - -incflo.verbose=1 +zhi.type = "slip_wall" +incflo.verbose=1 \ No newline at end of file diff --git a/test/test_files/ow_stokes/ow_stokes.inp b/test/test_files/ow_stokes/ow_stokes.inp index 9193f4b610..e411d9346c 100644 --- a/test/test_files/ow_stokes/ow_stokes.inp +++ b/test/test_files/ow_stokes/ow_stokes.inp @@ -20,9 +20,6 @@ time.checkpoint_interval = -1 # Steps between checkpoint files #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # PHYSICS # #.......................................# -incflo.use_godunov = 1 -incflo.diffusion_type = 2 -incflo.godunov_type="weno_z" transport.model = TwoPhaseTransport transport.viscosity_fluid1=0.0 transport.viscosity_fluid2=0.0 @@ -41,8 +38,6 @@ OceanWaves.Wave1.numerical_beach_length=2.0 OceanWaves.Wave1.numerical_beach_length_factor=2.0 MultiPhase.density_fluid1=1000. MultiPhase.density_fluid2=1. -MultiPhase.interface_smoothing=0 -MultiPhase.interface_smoothing_frequency=1 ICNS.source_terms = GravityForcing ICNS.use_perturb_pressure = true MultiPhase.verbose=1 @@ -58,15 +53,10 @@ geometry.prob_lo = 0.0 0. -1 # Lo corner coordinates geometry.prob_hi = 24. 1. 1 # Hi corner coordinates geometry.is_periodic = 0 1 0 # Periodicity x y z (0/1) -xlo.type = "pressure_inflow" +xlo.type = "wave_generation" xhi.type = "pressure_outflow" -xlo.vof_type = "zero_gradient" -xhi.vof_type = "zero_gradient" zlo.type = "slip_wall" -zhi.type = "pressure_outflow" -zlo.vof_type = "zero_gradient" -zhi.vof_type = "zero_gradient" - -incflo.verbose=1 +zhi.type = "slip_wall" +incflo.verbose=1 \ No newline at end of file diff --git a/test/test_files/ow_w2a/ow_w2a.inp b/test/test_files/ow_w2a/ow_w2a.inp index 001a335c70..60dad0c36e 100644 --- a/test/test_files/ow_w2a/ow_w2a.inp +++ b/test/test_files/ow_w2a/ow_w2a.inp @@ -14,7 +14,7 @@ time.cfl = 0.95 # CFL factor #.......................................# time.plot_interval = 10 # Steps between plot files time.checkpoint_interval = -1 # Steps between checkpoint files -io.outputs = density velocity p vof ow_levelset ow_velocity w2a_levelset w2a_velocity +io.outputs = density velocity p vof ow_levelset ow_velocity OceanWaves.label = W2A1 OceanWaves.W2A1.type = W2AWaves @@ -29,8 +29,6 @@ OceanWaves.W2A1.number_interp_above_surface = 5 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # PHYSICS # #.......................................# -incflo.use_godunov = 1 -incflo.godunov_type = weno_z incflo.mflux_type = minmod transport.model = TwoPhaseTransport transport.viscosity_fluid1=1.0e-3