Skip to content

Commit 97349ae

Browse files
mbkuhnmbkuhn
authored
Ocean Waves: relaxation zones in y (#1715)
* relaxation zones in y: compiling, passing existing unit tests * waves2amr lateral check in y * modify nonperiodic reg test * modify zones check * fix mixup between problo and probhi * proof of concept unit test * remove unnecessary has_outprofile variable * remove bad const
1 parent bb9646e commit 97349ae

File tree

5 files changed

+121
-98
lines changed

5 files changed

+121
-98
lines changed

‎amr-wind/ocean_waves/relaxation_zones/RelaxationZones.H‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,14 +27,15 @@ struct RelaxZonesBaseData
2727
amrex::Real beach_length{8.0};
2828
amrex::Real beach_length_factor{1.0};
2929

30+
amrex::Real zone_length_y{0.0};
31+
3032
amrex::Real current{0.0};
3133

3234
bool init_wave_field{false};
3335

3436
bool has_ramp{false};
3537

3638
bool has_beach{true};
37-
bool has_outprofile{false};
3839

3940
amrex::Real ramp_period;
4041

‎amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp‎

Lines changed: 61 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,6 @@ void read_inputs(
2525
pp.query("relax_zone_gen_length", wdata.gen_length);
2626
if (pp.contains("relax_zone_out_length")) {
2727
wdata.has_beach = false;
28-
wdata.has_outprofile = true;
2928
wdata.init_wave_field = true;
3029
pp.query("relax_zone_out_length", wdata.beach_length);
3130
} else {
@@ -34,6 +33,7 @@ void read_inputs(
3433
pp.query("numerical_beach_length_factor", wdata.beach_length_factor);
3534
pp.query("initialize_wave_field", wdata.init_wave_field);
3635
}
36+
pp.query("relax_zone_length_y", wdata.zone_length_y);
3737

3838
pp.query("current", wdata.current);
3939

@@ -120,17 +120,21 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata)
120120
const amrex::Real gen_length = wdata.gen_length;
121121
const amrex::Real beach_length = wdata.beach_length;
122122
const amrex::Real beach_length_factor = wdata.beach_length_factor;
123+
const amrex::Real zone_length_y = wdata.zone_length_y;
124+
const bool has_zone_y = zone_length_y > constants::EPS;
123125
const amrex::Real zsl = wdata.zsl;
124126
const amrex::Real current = wdata.current;
125127
const bool has_beach = wdata.has_beach;
126-
const bool has_outprofile = wdata.has_outprofile;
127128

128129
amrex::ParallelFor(
129130
velocity(lev), amrex::IntVect(0),
130131
[=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
131132
const amrex::Real x = amrex::min(
132133
amrex::max(problo[0] + (i + 0.5) * dx[0], problo[0]),
133134
probhi[0]);
135+
const amrex::Real y = amrex::min(
136+
amrex::max(problo[1] + (j + 0.5) * dx[1], problo[1]),
137+
probhi[1]);
134138
const amrex::Real z = amrex::min(
135139
amrex::max(problo[2] + (k + 0.5) * dx[2], problo[2]),
136140
probhi[2]);
@@ -141,128 +145,96 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata)
141145
const auto target_volfrac = target_volfrac_arrs[nbx];
142146
const auto target_vel = target_vel_arrs[nbx];
143147

148+
// Skip if in or near terrain
144149
bool in_or_near_terrain{false};
145150
if (terrain_exists) {
146151
in_or_near_terrain =
147152
(terrain_blank_flags[nbx](i, j, k) == 1 ||
148153
terrain_drag_flags[nbx](i, j, k) == 1);
149154
}
150155

151-
// Generation region
152-
if (x <= problo[0] + gen_length && !in_or_near_terrain) {
153-
const amrex::Real Gamma =
154-
utils::gamma_generate(x - problo[0], gen_length);
155-
// Get bounded new vof, incorporate with increment
156+
// Get gamma for each possible direction
157+
const amrex::Real Gamma_xlo =
158+
utils::gamma_generate(x - problo[0], gen_length);
159+
const amrex::Real Gamma_xhi = utils::gamma_absorb(
160+
x - (probhi[0] - beach_length), beach_length,
161+
beach_length_factor);
162+
amrex::Real Gamma_ylo = 1.;
163+
amrex::Real Gamma_yhi = 1.;
164+
if (has_zone_y) {
165+
Gamma_ylo =
166+
utils::gamma_generate(y - problo[1], zone_length_y);
167+
Gamma_yhi = utils::gamma_absorb(
168+
y - (probhi[1] - zone_length_y), zone_length_y, 1.0);
169+
}
170+
const amrex::Real Gamma = std::min(
171+
std::min(Gamma_xhi, Gamma_xlo),
172+
std::min(Gamma_yhi, Gamma_ylo));
173+
174+
// Skip if Gamma is close enough to 1
175+
bool outside_zones = Gamma + constants::EPS >= 1.;
176+
177+
if (!(outside_zones || in_or_near_terrain)) {
178+
// Create wave vector for generation, numerical beach
179+
const utils::WaveVec wave_sol{
180+
target_vel(i, j, k, 0), target_vel(i, j, k, 1),
181+
target_vel(i, j, k, 2), target_volfrac(i, j, k)};
182+
const utils::WaveVec quiescent{
183+
current, 0.0, 0.0,
184+
utils::free_surface_to_vof(zsl, z, dx[2])};
185+
const auto outlet = has_beach ? quiescent : wave_sol;
186+
187+
// Check for in beach, needed for velocity forcing
188+
const bool in_beach =
189+
has_beach && x + beach_length >= probhi[0];
190+
191+
// Harmonize between inlet/bulk profile and outlet profile
192+
const auto target_profile = utils::harmonize_profiles_1d(
193+
x, problo[0], gen_length, probhi[0], beach_length,
194+
wave_sol, wave_sol, outlet);
195+
196+
// Nudge solution toward target
156197
amrex::Real new_vof = utils::combine_linear(
157-
Gamma, target_volfrac(i, j, k), volfrac(i, j, k));
198+
Gamma, target_profile[3], volfrac(i, j, k));
158199
new_vof = (new_vof > 1. - vof_tiny)
159200
? 1.0
160201
: (new_vof < vof_tiny ? 0.0 : new_vof);
161202
const amrex::Real dvf = new_vof - volfrac(i, j, k);
162203
volfrac(i, j, k) += rampf * dvf;
163-
// Force liquid velocity only if target vof present
164-
const amrex::Real fvel_liq =
204+
// Liquid velocity forced only where velocity is known
205+
// - in most of domain, that is where target vof is nonzero
206+
// - in numerical beach, that is anywhere
207+
amrex::Real fvel_liq =
165208
(target_volfrac(i, j, k) > vof_tiny) ? 1.0 : 0.0;
209+
fvel_liq = in_beach ? 1.0 : fvel_liq;
166210
amrex::Real rho_ = rho1 * volfrac(i, j, k) +
167211
rho2 * (1.0 - volfrac(i, j, k));
168212
for (int n = 0; n < vel.ncomp; ++n) {
169213
// Get updated liquid velocity
170214
amrex::Real vel_liq = vel(i, j, k, n);
171215
const amrex::Real dvel_liq =
172216
utils::combine_linear(
173-
Gamma, target_vel(i, j, k, n), vel_liq) -
217+
Gamma, target_profile[n], vel_liq) -
174218
vel_liq;
175219
vel_liq += rampf * fvel_liq * dvel_liq;
176220
// If liquid was added, that liquid has target_vel
177221
amrex::Real integrated_vel_liq =
178222
volfrac(i, j, k) * vel_liq;
179223
integrated_vel_liq +=
180224
rampf * fvel_liq * amrex::max(0.0, dvf) *
181-
(target_vel(i, j, k, n) - vel(i, j, k, n));
225+
(target_profile[n] - vel(i, j, k, n));
182226
// Update overall velocity using momentum
183227
vel(i, j, k, n) =
184228
(rho1 * integrated_vel_liq +
185229
rho2 * (1. - volfrac(i, j, k)) * vel(i, j, k, n)) /
186230
rho_;
187231
}
188-
}
189-
// Outlet region
190-
if (x + beach_length >= probhi[0] && !in_or_near_terrain) {
191-
const amrex::Real Gamma = utils::gamma_absorb(
192-
x - (probhi[0] - beach_length), beach_length,
193-
beach_length_factor);
194-
// Numerical beach (sponge layer)
195-
if (has_beach) {
196-
// Get bounded new vof, save increment
197-
amrex::Real new_vof = utils::combine_linear(
198-
Gamma, utils::free_surface_to_vof(zsl, z, dx[2]),
199-
volfrac(i, j, k));
200-
new_vof = (new_vof > 1. - vof_tiny)
201-
? 1.0
202-
: (new_vof < vof_tiny ? 0.0 : new_vof);
203-
const amrex::Real dvf = new_vof - volfrac(i, j, k);
204-
volfrac(i, j, k) = new_vof;
205-
// Conserve momentum when density changes
206-
amrex::Real rho_ = rho1 * volfrac(i, j, k) +
207-
rho2 * (1.0 - volfrac(i, j, k));
208-
// Target vel is current, 0, 0
209-
amrex::Real out_vel{0.0};
210-
for (int n = 0; n < vel.ncomp; ++n) {
211-
if (n == 0) {
212-
out_vel = current;
213-
} else {
214-
out_vel = 0.0;
215-
}
216-
const amrex::Real vel_liq = utils::combine_linear(
217-
Gamma, out_vel, vel(i, j, k, n));
218-
amrex::Real integrated_vel_liq =
219-
volfrac(i, j, k) * vel_liq;
220-
integrated_vel_liq += amrex::max(0.0, dvf) *
221-
(out_vel - vel(i, j, k, n));
222-
vel(i, j, k, n) = (rho1 * integrated_vel_liq +
223-
rho2 * (1. - volfrac(i, j, k)) *
224-
vel(i, j, k, n)) /
225-
rho_;
226-
}
227-
}
228-
// Forcing to wave profile instead
229-
if (has_outprofile) {
230-
// Same steps as in wave generation region
231-
amrex::Real new_vof = utils::combine_linear(
232-
Gamma, target_volfrac(i, j, k), volfrac(i, j, k));
233-
new_vof = (new_vof > 1. - vof_tiny)
234-
? 1.0
235-
: (new_vof < vof_tiny ? 0.0 : new_vof);
236-
const amrex::Real dvf = new_vof - volfrac(i, j, k);
237-
volfrac(i, j, k) += dvf;
238-
const amrex::Real fvel_liq =
239-
(target_volfrac(i, j, k) > vof_tiny) ? 1.0 : 0.0;
240-
amrex::Real rho_ = rho1 * volfrac(i, j, k) +
241-
rho2 * (1.0 - volfrac(i, j, k));
242-
for (int n = 0; n < vel.ncomp; ++n) {
243-
amrex::Real vel_liq = vel(i, j, k, n);
244-
const amrex::Real dvel_liq =
245-
utils::combine_linear(
246-
Gamma, target_vel(i, j, k, n), vel_liq) -
247-
vel_liq;
248-
vel_liq += rampf * fvel_liq * dvel_liq;
249-
amrex::Real integrated_vel_liq =
250-
volfrac(i, j, k) * vel_liq;
251-
integrated_vel_liq +=
252-
rampf * fvel_liq * amrex::max(0.0, dvf) *
253-
(target_vel(i, j, k, n) - vel(i, j, k, n));
254-
vel(i, j, k, n) = (rho1 * integrated_vel_liq +
255-
rho2 * (1. - volfrac(i, j, k)) *
256-
vel(i, j, k, n)) /
257-
rho_;
258-
}
259-
}
260-
}
261232

262-
// Make sure that density is updated before entering the
263-
// solution
264-
rho(i, j, k) =
265-
rho1 * volfrac(i, j, k) + rho2 * (1. - volfrac(i, j, k));
233+
// Make sure that density is updated before entering the
234+
// solution
235+
rho(i, j, k) = rho1 * volfrac(i, j, k) +
236+
rho2 * (1. - volfrac(i, j, k));
237+
}
266238
});
267239
}
268240
amrex::Gpu::streamSynchronize();

‎amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H‎

Lines changed: 17 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -813,6 +813,8 @@ struct UpdateTargetFieldsOp<W2AWaves>
813813
bool flag_z = false;
814814
bool flag_xlo = false;
815815
bool flag_xhi = false;
816+
bool flag_ylo = false;
817+
bool flag_yhi = false;
816818
// Get heights for this processor, check overlap in z
817819
flag_z =
818820
(interp_to_mfab::get_local_height_indices(
@@ -821,17 +823,27 @@ struct UpdateTargetFieldsOp<W2AWaves>
821823
// No overlap from heights definitely means no interp
822824

823825
// Check lateral bounds (in x)
824-
const int dir = 0;
825826
flag_xlo =
826827
(interp_to_mfab::check_lateral_overlap_lo(
827-
wdata.gen_length, dir, ow_velocity.vec_ptrs(), geom) ==
828+
wdata.gen_length, 0, ow_velocity.vec_ptrs(), geom) ==
828829
1);
829830
// No overlap with gen region means no interp, unless ...
830-
if (wdata.has_outprofile) {
831+
if (!wdata.has_beach) {
831832
// ... if overlap exists here, needing interp
832833
flag_xhi =
833834
(interp_to_mfab::check_lateral_overlap_hi(
834-
wdata.beach_length, dir, ow_velocity.vec_ptrs(),
835+
wdata.beach_length, 0, ow_velocity.vec_ptrs(),
836+
geom) == 1);
837+
}
838+
// Then in y, if y zones are present
839+
if (wdata.zone_length_y > constants::EPS) {
840+
flag_ylo =
841+
(interp_to_mfab::check_lateral_overlap_lo(
842+
wdata.zone_length_y, 1, ow_velocity.vec_ptrs(),
843+
geom) == 1);
844+
flag_yhi =
845+
(interp_to_mfab::check_lateral_overlap_hi(
846+
wdata.zone_length_y, 1, ow_velocity.vec_ptrs(),
835847
geom) == 1);
836848
}
837849

@@ -842,7 +854,7 @@ struct UpdateTargetFieldsOp<W2AWaves>
842854
flag_xhi = true;
843855
}
844856

845-
if (flag_z && (flag_xlo || flag_xhi)) {
857+
if (flag_z && (flag_xlo || flag_xhi || flag_ylo || flag_yhi)) {
846858
// Interpolation is needed
847859
wdata.do_interp = true;
848860
// Do resizing

‎test/test_files/ow_w2a_nonperiodic/ow_w2a_nonperiodic.inp‎

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@ OceanWaves.W2A1.HOS_modes_filename = ../ow_w2a/modes_HOS_SWENSE.dat
2222
OceanWaves.W2A1.relax_zone_gen_length=200
2323
OceanWaves.W2A1.numerical_beach_length=200
2424
OceanWaves.W2A1.numerical_beach_length_factor=2.0
25+
OceanWaves.W2A1.relax_zone_length_y=50
2526
# These variables should change with resolution in z
2627
OceanWaves.W2A1.number_interp_points_in_z = 35
2728
OceanWaves.W2A1.interp_spacing_at_surface = 1.
@@ -57,10 +58,12 @@ geometry.is_periodic = 0 0 0 # Periodicity x y z (0/1)
5758

5859
xlo.type = wave_generation
5960
xhi.type = pressure_outflow
60-
ylo.type = slip_wall
61-
yhi.type = slip_wall
61+
ylo.type = wave_generation
62+
yhi.type = wave_generation
6263
zlo.type = slip_wall
6364
zhi.type = slip_wall
6465

6566
# density at inflow condition must match the gas density specified above
66-
xlo.density = 1.25
67+
xlo.density = 1.25
68+
ylo.density = 1.25
69+
yhi.density = 1.25

‎unit_tests/ocean_waves/test_wave_utils.cpp‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -196,4 +196,39 @@ TEST_F(WaveUtilsTest, harmonize_profiles)
196196
}
197197
}
198198

199+
// Proof of concept test for multiple gammas in x and y
200+
TEST_F(WaveUtilsTest, gamma_xy)
201+
{
202+
constexpr amrex::Real tol = 1e-12;
203+
constexpr amrex::Real zone_length = 0.25;
204+
205+
const amrex::Vector<amrex::Real> x{
206+
0., 1., 0.5, 0.5, zone_length, zone_length, 0.5 * zone_length};
207+
const amrex::Vector<amrex::Real> y{
208+
0.5, 0.5, 0., 1., zone_length, 0.5, 0.5 * zone_length};
209+
const amrex::Real gm =
210+
1.0 -
211+
(std::exp(std::pow(1. - x[x.size() - 1] / zone_length, 3.5)) - 1.0) /
212+
(std::exp(1.0) - 1.0);
213+
const amrex::Vector<amrex::Real> gold_gamma{0., 0., 0., 0., 1., 1., gm};
214+
215+
for (int n = 0; n < x.size(); ++n) {
216+
const amrex::Real gamma_xlo =
217+
amr_wind::ocean_waves::utils::gamma_generate(
218+
x[n] - 0., zone_length);
219+
const amrex::Real gamma_xhi =
220+
amr_wind::ocean_waves::utils::gamma_absorb(
221+
x[n] - (1. - zone_length), zone_length, 1.);
222+
const amrex::Real gamma_ylo =
223+
amr_wind::ocean_waves::utils::gamma_generate(
224+
y[n] - 0., zone_length);
225+
const amrex::Real gamma_yhi =
226+
amr_wind::ocean_waves::utils::gamma_absorb(
227+
y[n] - (1. - zone_length), zone_length, 1.);
228+
const amrex::Real gamma = std::min(
229+
std::min(gamma_xhi, gamma_xlo), std::min(gamma_yhi, gamma_ylo));
230+
EXPECT_NEAR(gamma, gold_gamma[n], tol);
231+
}
232+
}
233+
199234
} // namespace amr_wind_tests

0 commit comments

Comments
 (0)