Skip to content

Commit

Permalink
Improvements to wave generation - ocean waves (Exawind#1264)
Browse files Browse the repository at this point in the history 
* put loops where they should be already
* set up rampf properly for generation region
* duplicate changes for outprofile
* adding limitation to velocity forcing, rephrasing things
* tweak for when liquid is added
* remove false default, update docs
* correcting/updating some stuff in the reg tests while we're at it
  • Loading branch information
@mbkuhn
mbkuhn authored Oct 2, 2024
1 parent 091f6cc commit d4aa6fc
Show file tree
Hide file tree
Showing 6 changed files with 81 additions and 92 deletions.
2 changes: 1 addition & 1 deletion amr-wind/ocean_waves/relaxation_zones/RelaxationZones.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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};
};
Expand Down
133 changes: 72 additions & 61 deletions amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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");
Expand All @@ -59,6 +59,7 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata)
const auto& mphase = sim.physics_manager().get<MultiPhase>();
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);
Expand Down Expand Up @@ -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) *
Expand All @@ -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_;
}
}
}

Expand Down
2 changes: 1 addition & 1 deletion docs/sphinx/user/inputs_ocean_waves.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
16 changes: 3 additions & 13 deletions test/test_files/ow_linear/ow_linear.inp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
16 changes: 3 additions & 13 deletions test/test_files/ow_stokes/ow_stokes.inp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
4 changes: 1 addition & 3 deletions test/test_files/ow_w2a/ow_w2a.inp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down

0 comments on commit d4aa6fc

Please sign in to comment.