From 3dca78df13cce612065692ee47d576701a2edbbd Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 13:01:41 -0700 Subject: [PATCH 01/31] init w2a waves as intended --- .../relaxation_zones/waves2amr_ops.H | 87 ++++++++++++------- 1 file changed, 54 insertions(+), 33 deletions(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index 0c11fca116..0f94cd23c0 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -414,37 +414,19 @@ struct InitDataOp #ifdef AMR_WIND_USE_W2A auto& wdata = data.meta(); - bool init_waves = wdata.init_wave_field; auto& sim = data.sim(); + // Fill ow fields, then populate flow fields according to setup + auto& ow_levelset = sim.repo().get_field("ow_levelset"); + auto& ow_velocity = sim.repo().get_field("ow_velocity"); auto& levelset = sim.repo().get_field("levelset"); auto& velocity = sim.repo().get_field("velocity"); const auto& problo = geom.ProbLoArray(); + const auto& probhi = geom.ProbHiArray(); const auto& dx = geom.CellSizeArray(); - // Blank initialization if asked for - if (!init_waves) { - // Loop to populate field data - for (amrex::MFIter mfi(levelset(level)); mfi.isValid(); ++mfi) { - auto phi = levelset(level).array(mfi); - auto vel = velocity(level).array(mfi); - const auto& vbx = mfi.validbox(); - const auto& gbx = grow(vbx, 3); - const amrex::Real zsl = wdata.zsl; - amrex::ParallelFor( - gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; - phi(i, j, k) = zsl - z; - vel(i, j, k, 0) = 0.0; - vel(i, j, k, 1) = 0.0; - vel(i, j, k, 2) = 0.0; - }); - } - return; - } - // Set t_last to 0.0 to signify information read in wdata.t_last = 0.0; @@ -457,20 +439,59 @@ struct InitDataOp // Interpolate to MultiFabs (one level at a time) interp_to_mfab::interp_eta_to_levelset_multifab( wdata.n0, wdata.n1, wdata.dx0, wdata.dx1, wdata.zsl, - wdata.sp_eta_vec, levelset(level), problo, dx); + wdata.sp_eta_vec, ow_levelset(level), problo, dx); interp_to_mfab::interp_velocity_to_multifab( wdata.n0, wdata.n1, wdata.dx0, wdata.dx1, indvec, wdata.hvec, - wdata.sp_u_vec, wdata.sp_v_vec, wdata.sp_w_vec, velocity(level), + wdata.sp_u_vec, wdata.sp_v_vec, wdata.sp_w_vec, ow_velocity(level), problo, dx); - // Zero velocity in pure liquid cells - postprocess_velocity_mfab_liquid(velocity(level), levelset(level), dx); - - // Copy to ow fields for future interpolation - auto& ow_levelset = sim.repo().get_field("ow_levelset"); - auto& ow_velocity = sim.repo().get_field("ow_velocity"); - amrex::MultiFab::Copy(ow_levelset(level), levelset(level), 0, 0, 1, 3); - amrex::MultiFab::Copy( - ow_velocity(level), velocity(level), 0, 0, AMREX_SPACEDIM, 3); + // Zero velocity in pure air cells + postprocess_velocity_mfab_liquid( + ow_velocity(level), ow_levelset(level), dx); + + // Populate flow fields according to intended forcing and init setup + for (amrex::MFIter mfi(levelset(level)); mfi.isValid(); ++mfi) { + const auto& ow_phi = ow_levelset(level).array(mfi); + const auto& ow_vel = ow_velocity(level).array(mfi); + const auto& phi = levelset(level).array(mfi); + const auto& vel = velocity(level).array(mfi); + const bool has_beach = wdata.has_beach; + const bool init_wave_field = wdata.init_wave_field; + const auto& gbx3 = mfi.growntilebox(3); + + amrex::ParallelFor( + gbx3, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + + const amrex::Real gen_length = wdata.gen_length; + const amrex::Real beach_length = wdata.beach_length; + const amrex::Real zero_sea_level = wdata.zsl; + // Wave profile + const utils::WaveVec wave_sol{ + ow_vel(i, j, k, 0), ow_vel(i, j, k, 1), + ow_vel(i, j, k, 2), ow_phi(i, j, k) + z}; + // Quiescent profile + const utils::WaveVec quiescent{ + 0.0, 0.0, 0.0, zero_sea_level}; + + // Specify initial state for each region of domain + const auto bulk = init_wave_field ? wave_sol : quiescent; + const auto outlet = has_beach ? quiescent : wave_sol; + + const auto local_profile = utils::harmonize_profiles_1d( + x, problo[0], gen_length, probhi[0], beach_length, + wave_sol, bulk, outlet); + + phi(i, j, k) = local_profile[3] - z; + const amrex::Real cell_length_2D = + std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); + if (phi(i, j, k) + cell_length_2D >= 0) { + vel(i, j, k, 0) = local_profile[0]; + vel(i, j, k, 1) = local_profile[1]; + vel(i, j, k, 2) = local_profile[2]; + } + }); + } // Start w2a fields at 0, some areas will not be modified auto& w2a_levelset = sim.repo().get_field("w2a_levelset"); From d253836d294a5005fe7344d6876795ac36014188 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 13:16:53 -0700 Subject: [PATCH 02/31] set up abl multiphase reg test for w2a, shows obvious flaws --- test/CMakeLists.txt | 1 + .../abl_multiphase_w2a/abl_multiphase_w2a.inp | 88 +++++++++++++++++++ .../abl_multiphase_w2a/static_box.refine | 4 + 3 files changed, 93 insertions(+) create mode 100644 test/test_files/abl_multiphase_w2a/abl_multiphase_w2a.inp create mode 100644 test/test_files/abl_multiphase_w2a/static_box.refine diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5f8a103049..8799e0eeee 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -318,6 +318,7 @@ endif() if(AMR_WIND_ENABLE_W2A) add_test_re(ow_w2a) + add_test_re(abl_multiphase_w2a) endif() #============================================================================= diff --git a/test/test_files/abl_multiphase_w2a/abl_multiphase_w2a.inp b/test/test_files/abl_multiphase_w2a/abl_multiphase_w2a.inp new file mode 100644 index 0000000000..c2f134669a --- /dev/null +++ b/test/test_files/abl_multiphase_w2a/abl_multiphase_w2a.inp @@ -0,0 +1,88 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 12.25 # Max (simulated) time to evolve +time.max_step = 10 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 1.0 # Use this constant dt if > 0 +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 +io.outputs = density velocity p vof ow_levelset ow_velocity + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.mflux_type = minmod +transport.model = TwoPhaseTransport +transport.viscosity_fluid1=1.0e-3 +transport.viscosity_fluid2=1.8e-5 +turbulence.model = Laminar + +incflo.physics = MultiPhase ABL OceanWaves +MultiPhase.density_fluid1=1000 +MultiPhase.density_fluid2=1.25 +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing ABLForcing 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 +ABLForcing.abl_forcing_height = 200 +ABLForcing.abl_forcing_off_height = 10.0 +ABLForcing.abl_forcing_ramp_height = 10.0 +ABLForcing.abl_forcing_band = 1 +incflo.velocity = 6.128355544951824 5.142300877492314 0.0 + +BoussinesqBuoyancy.reference_temperature = 263.5 +ABL.reference_temperature = 263.5 +ABL.temperature_heights = -5.0 0.0 100 400.0 1001.0 +ABL.temperature_values = 263.5 265.0 265.0 268.0 270.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 + +OceanWaves.label = W2A1 +OceanWaves.W2A1.type = W2AWaves +OceanWaves.W2A1.HOS_modes_filename = ../ow_w2a/modes_HOS_SWENSE.dat +OceanWaves.W2A1.relax_zone_gen_length=900 +OceanWaves.W2A1.relax_zone_out_length=900 +# These variables should change with resolution in z +OceanWaves.W2A1.number_interp_points_in_z = 35 +OceanWaves.W2A1.interp_spacing_at_surface = 10. +OceanWaves.W2A1.number_interp_above_surface = 5 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 64 64 32 +amr.max_level = 1 +tagging.labels = sr +tagging.sr.type = CartBoxRefinement +tagging.sr.static_refinement_def = static_box.refine +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0.0 0.0 -1000 # Lo corner coordinates +geometry.prob_hi = 9.3717632253682877E+03 9.3717632253682877E+03 1000 # Hi corner coordinates +geometry.is_periodic = 1 1 0 # Periodicity x y z (0/1) + +zlo.type = slip_wall +zhi.type = slip_wall diff --git a/test/test_files/abl_multiphase_w2a/static_box.refine b/test/test_files/abl_multiphase_w2a/static_box.refine new file mode 100644 index 0000000000..09d360cf70 --- /dev/null +++ b/test/test_files/abl_multiphase_w2a/static_box.refine @@ -0,0 +1,4 @@ +1 +1 +0 0 -100 9.3717632253682877E+03 9.3717632253682877E+03 100 + From d824c952fac1af0a3a0d02ecedadbedaa8f51390 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 13:40:07 -0700 Subject: [PATCH 03/31] move variable assignment outside parallelfor --- amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index 0f94cd23c0..8fb10515c7 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -454,6 +454,11 @@ struct InitDataOp const auto& ow_vel = ow_velocity(level).array(mfi); const auto& phi = levelset(level).array(mfi); const auto& vel = velocity(level).array(mfi); + + const amrex::Real gen_length = wdata.gen_length; + const amrex::Real beach_length = wdata.beach_length; + const amrex::Real zero_sea_level = wdata.zsl; + const bool has_beach = wdata.has_beach; const bool init_wave_field = wdata.init_wave_field; const auto& gbx3 = mfi.growntilebox(3); @@ -463,9 +468,6 @@ struct InitDataOp const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; - const amrex::Real gen_length = wdata.gen_length; - const amrex::Real beach_length = wdata.beach_length; - const amrex::Real zero_sea_level = wdata.zsl; // Wave profile const utils::WaveVec wave_sol{ ow_vel(i, j, k, 0), ow_vel(i, j, k, 1), From cbce3524d63ffa265f7ebcfac35bbd797c62d097 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 13:52:15 -0700 Subject: [PATCH 04/31] tweak time stepping settings --- test/test_files/abl_multiphase_w2a/abl_multiphase_w2a.inp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_files/abl_multiphase_w2a/abl_multiphase_w2a.inp b/test/test_files/abl_multiphase_w2a/abl_multiphase_w2a.inp index c2f134669a..88351f32f5 100644 --- a/test/test_files/abl_multiphase_w2a/abl_multiphase_w2a.inp +++ b/test/test_files/abl_multiphase_w2a/abl_multiphase_w2a.inp @@ -7,13 +7,14 @@ time.max_step = 10 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# # TIME STEP COMPUTATION # #.......................................# -time.fixed_dt = 1.0 # Use this constant dt if > 0 +time.fixed_dt = 0.5 # Use this constant dt if > 0 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.use_force_cfl = false io.outputs = density velocity p vof ow_levelset ow_velocity #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# From c9c73b22a97bac2ea37d0491a222d4e7281ae77a Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 14:06:27 -0700 Subject: [PATCH 05/31] documentation (for arguments now in reg test) --- docs/sphinx/user/inputs_Momentum_Sources.rst | 27 ++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/docs/sphinx/user/inputs_Momentum_Sources.rst b/docs/sphinx/user/inputs_Momentum_Sources.rst index 8e70598ddc..e7d46b4024 100644 --- a/docs/sphinx/user/inputs_Momentum_Sources.rst +++ b/docs/sphinx/user/inputs_Momentum_Sources.rst @@ -121,6 +121,33 @@ Section: Momentum Sources The start time for writing to the forcing timetable output file. The default is 0. +.. input_param:: ABLForcing.abl_forcing_off_height + + **type:** Real, required for multiphase simulations with ABL + + This parameter indicates the vertical distance above the water level that the ABL + forcing term should be turned off. This tuning parameter is used to avoid applying + the ABL forcing to ocean waves. This is not used when the volume fraction field (vof) + is not present in the simulation. + +.. input_param:: ABLForcing.abl_forcing_off_height + + **type:** Real, required for multiphase simulations with ABL + + This parameter indicates the vertical distance above the water level and the "off height" + that the ABL forcing term should be ramped from zero to full strength. This is not used + when the volume fraction field (vof) is not present in the simulation. + +.. input_param:: ABLForcing.abl_forcing_band + + **type:** Real, optional for multiphase simulations with ABL + + This parameter is an additional safeguard against applying ABL forcing within the waves. + This specifies the number of computational cells in a band around the air-water interface + that the ABL forcing should be deactivated. While the other arguments relate to the height coordinate + within the domain, this argument is relative to the actual position of water in the simulation. + The default value is 2. + .. input_param:: BodyForce.type **type:** String, optional From 6fe11fea741dca554aa0f8a038c22817e30aa1de Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 14:18:32 -0700 Subject: [PATCH 06/31] docs fix --- docs/sphinx/user/inputs_Momentum_Sources.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/sphinx/user/inputs_Momentum_Sources.rst b/docs/sphinx/user/inputs_Momentum_Sources.rst index e7d46b4024..b75bb53648 100644 --- a/docs/sphinx/user/inputs_Momentum_Sources.rst +++ b/docs/sphinx/user/inputs_Momentum_Sources.rst @@ -130,7 +130,7 @@ Section: Momentum Sources the ABL forcing to ocean waves. This is not used when the volume fraction field (vof) is not present in the simulation. -.. input_param:: ABLForcing.abl_forcing_off_height +.. input_param:: ABLForcing.abl_forcing_ramp_height **type:** Real, required for multiphase simulations with ABL From 1eb1d53b97ee270da253260f9d2e2b1f0604ba83 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 15:23:48 -0700 Subject: [PATCH 07/31] get those MFIters outta here --- .../relaxation_zones/linear_waves_ops.H | 230 ++++++++--------- .../relaxation_zones/relaxation_zones_ops.cpp | 244 +++++++++--------- .../relaxation_zones/stokes_waves_ops.H | 171 ++++++------ .../relaxation_zones/waves2amr_ops.H | 160 ++++++------ 4 files changed, 388 insertions(+), 417 deletions(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H index f1137bafbd..8b624d3a5b 100644 --- a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H @@ -45,71 +45,68 @@ struct InitDataOp const auto& problo = geom.ProbLoArray(); const auto& probhi = geom.ProbHiArray(); const auto& dx = geom.CellSizeArray(); - for (amrex::MFIter mfi(m_levelset(level)); mfi.isValid(); ++mfi) { - const auto& phi = m_levelset(level).array(mfi); - const auto& vel = m_velocity(level).array(mfi); - - const amrex::Real zero_sea_level = wdata.zsl; - const amrex::Real gen_length = wdata.gen_length; - const amrex::Real beach_length = wdata.beach_length; - const amrex::Real g = wdata.g; - const bool has_beach = wdata.has_beach; - const bool init_wave_field = wdata.init_wave_field; - const auto& gbx3 = mfi.growntilebox(3); - - const amrex::Real wave_height = wdata.wave_height; - const amrex::Real wave_length = wdata.wave_length; - const amrex::Real water_depth = wdata.water_depth; - amrex::ParallelFor( - gbx3, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; - const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; - - const amrex::Real wave_number = 2. * M_PI / wave_length; - const amrex::Real omega = std::pow( - wave_number * g * std::tanh(wave_number * water_depth), - 0.5); - const amrex::Real phase = wave_number * x; - - // Wave profile - const amrex::Real eta_w = - wave_height / 2.0 * std::cos(phase) + zero_sea_level; - const amrex::Real u_w = - omega * wave_height / 2.0 * - std::cosh( - wave_number * (z - zero_sea_level + water_depth)) / - std::sinh(wave_number * water_depth) * std::cos(phase); - const amrex::Real v_w = 0.0; - const amrex::Real w_w = - omega * wave_height / 2.0 * - std::sinh( - wave_number * (z - zero_sea_level + water_depth)) / - std::sinh(wave_number * water_depth) * std::sin(phase); - const utils::WaveVec wave_sol{u_w, v_w, w_w, eta_w}; - - // Quiescent profile - const utils::WaveVec quiescent{ - 0.0, 0.0, 0.0, zero_sea_level}; - - // Specify initial state for each region of domain - const auto bulk = init_wave_field ? wave_sol : quiescent; - const auto outlet = has_beach ? quiescent : wave_sol; - - const auto local_profile = utils::harmonize_profiles_1d( - x, problo[0], gen_length, probhi[0], beach_length, - wave_sol, bulk, outlet); - - phi(i, j, k) = local_profile[3] - z; - const amrex::Real cell_length_2D = - std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi(i, j, k) + cell_length_2D >= 0) { - vel(i, j, k, 0) = local_profile[0]; - vel(i, j, k, 1) = local_profile[1]; - vel(i, j, k, 2) = local_profile[2]; - } - }); - } + const auto& phi = m_levelset(level).arrays(); + const auto& vel = m_velocity(level).arrays(); + + const amrex::Real zero_sea_level = wdata.zsl; + const amrex::Real gen_length = wdata.gen_length; + const amrex::Real beach_length = wdata.beach_length; + const amrex::Real g = wdata.g; + const bool has_beach = wdata.has_beach; + const bool init_wave_field = wdata.init_wave_field; + + const amrex::Real wave_height = wdata.wave_height; + const amrex::Real wave_length = wdata.wave_length; + const amrex::Real water_depth = wdata.water_depth; + amrex::ParallelFor( + m_velocity(level), amrex::IntVect(3), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + + const amrex::Real wave_number = 2. * M_PI / wave_length; + const amrex::Real omega = std::pow( + wave_number * g * std::tanh(wave_number * water_depth), + 0.5); + const amrex::Real phase = wave_number * x; + + // Wave profile + const amrex::Real eta_w = + wave_height / 2.0 * std::cos(phase) + zero_sea_level; + const amrex::Real u_w = + omega * wave_height / 2.0 * + std::cosh( + wave_number * (z - zero_sea_level + water_depth)) / + std::sinh(wave_number * water_depth) * std::cos(phase); + const amrex::Real v_w = 0.0; + const amrex::Real w_w = + omega * wave_height / 2.0 * + std::sinh( + wave_number * (z - zero_sea_level + water_depth)) / + std::sinh(wave_number * water_depth) * std::sin(phase); + const utils::WaveVec wave_sol{u_w, v_w, w_w, eta_w}; + + // Quiescent profile + const utils::WaveVec quiescent{0.0, 0.0, 0.0, zero_sea_level}; + + // Specify initial state for each region of domain + const auto bulk = init_wave_field ? wave_sol : quiescent; + const auto outlet = has_beach ? quiescent : wave_sol; + + const auto local_profile = utils::harmonize_profiles_1d( + x, problo[0], gen_length, probhi[0], beach_length, wave_sol, + bulk, outlet); + + phi[nbx](i, j, k) = local_profile[3] - z; + const amrex::Real cell_length_2D = + std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); + if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + vel[nbx](i, j, k, 0) = local_profile[0]; + vel[nbx](i, j, k, 1) = local_profile[1]; + vel[nbx](i, j, k, 2) = local_profile[2]; + } + }); } }; @@ -135,61 +132,56 @@ struct UpdateRelaxZonesOp const auto& problo = geom[lev].ProbLoArray(); const auto& dx = geom[lev].CellSizeArray(); - for (amrex::MFIter mfi(m_ow_levelset(lev)); mfi.isValid(); ++mfi) { - const auto& phi = m_ow_levelset(lev).array(mfi); - const auto& vel = m_ow_velocity(lev).array(mfi); - - const amrex::Real wave_height = wdata.wave_height; - const amrex::Real wave_length = wdata.wave_length; - const amrex::Real water_depth = wdata.water_depth; - const amrex::Real zero_sea_level = wdata.zsl; - const amrex::Real g = wdata.g; - - const auto& gbx = mfi.growntilebox(3); - amrex::ParallelFor( - gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; - const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; - - const amrex::Real wave_number = 2. * M_PI / wave_length; - const amrex::Real omega = std::pow( - wave_number * g * - std::tanh(wave_number * water_depth), - 0.5); - const amrex::Real phase = - wave_number * x - omega * time; - - const amrex::Real eta = - wave_height / 2.0 * std::cos(phase) + - zero_sea_level; - - phi(i, j, k) = eta - z; - - const amrex::Real cell_length_2D = - std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi(i, j, k) + cell_length_2D >= 0) { - vel(i, j, k, 0) = - omega * wave_height / 2.0 * - std::cosh( - wave_number * - (z - zero_sea_level + water_depth)) / - std::sinh(wave_number * water_depth) * - std::cos(phase); - vel(i, j, k, 1) = 0.0; - vel(i, j, k, 2) = - omega * wave_height / 2.0 * - std::sinh( - wave_number * - (z - zero_sea_level + water_depth)) / - std::sinh(wave_number * water_depth) * - std::sin(phase); - } else { - vel(i, j, k, 0) = 0.; - vel(i, j, k, 1) = 0.; - vel(i, j, k, 2) = 0.; - } - }); - } + const auto& phi = m_ow_levelset(lev).arrays(); + const auto& vel = m_ow_velocity(lev).arrays(); + + const amrex::Real wave_height = wdata.wave_height; + const amrex::Real wave_length = wdata.wave_length; + const amrex::Real water_depth = wdata.water_depth; + const amrex::Real zero_sea_level = wdata.zsl; + const amrex::Real g = wdata.g; + + amrex::ParallelFor( + m_ow_velocity(lev), amrex::IntVect(3), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + + const amrex::Real wave_number = 2. * M_PI / wave_length; + const amrex::Real omega = std::pow( + wave_number * g * std::tanh(wave_number * water_depth), + 0.5); + const amrex::Real phase = wave_number * x - omega * time; + + const amrex::Real eta = + wave_height / 2.0 * std::cos(phase) + zero_sea_level; + + phi[nbx](i, j, k) = eta - z; + + const amrex::Real cell_length_2D = + std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); + if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + vel[nbx](i, j, k, 0) = + omega * wave_height / 2.0 * + std::cosh( + wave_number * + (z - zero_sea_level + water_depth)) / + std::sinh(wave_number * water_depth) * + std::cos(phase); + vel[nbx](i, j, k, 1) = 0.0; + vel[nbx](i, j, k, 2) = + omega * wave_height / 2.0 * + std::sinh( + wave_number * + (z - zero_sea_level + water_depth)) / + std::sinh(wave_number * water_depth) * + std::sin(phase); + } else { + vel[nbx](i, j, k, 0) = 0.; + vel[nbx](i, j, k, 1) = 0.; + vel[nbx](i, j, k, 2) = 0.; + } + }); } } }; 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 d64474ea8c..4c98d83a0b 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp @@ -60,21 +60,16 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata) constexpr amrex::Real vof_tiny = 1e-12; for (int lev = 0; lev < nlevels; ++lev) { - auto& ls = m_ow_levelset(lev); - auto& target_vof = m_ow_vof(lev); const auto& dx = geom[lev].CellSizeArray(); - - for (amrex::MFIter mfi(ls); mfi.isValid(); ++mfi) { - const auto& gbx = mfi.growntilebox(2); - const amrex::Array4& phi = ls.array(mfi); - const amrex::Array4& volfrac = target_vof.array(mfi); - const amrex::Real eps = 2. * std::cbrt(dx[0] * dx[1] * dx[2]); - amrex::ParallelFor( - gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - volfrac(i, j, k) = - multiphase::levelset_to_vof(i, j, k, eps, phi); - }); - } + const auto target_phi = m_ow_levelset(lev).arrays(); + auto target_volfrac = m_ow_vof(lev).arrays(); + const amrex::Real eps = 2. * std::cbrt(dx[0] * dx[1] * dx[2]); + amrex::ParallelFor( + m_ow_vof(lev), amrex::IntVect(2), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + target_volfrac[nbx](i, j, k) = + multiphase::levelset_to_vof(i, j, k, eps, target_phi[nbx]); + }); } // Get time @@ -87,147 +82,144 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata) auto& density = sim.repo().get_field("density"); for (int lev = 0; lev < nlevels; ++lev) { - for (amrex::MFIter mfi(vof(lev)); mfi.isValid(); ++mfi) { - const auto& gbx = mfi.growntilebox(2); - const auto& dx = geom[lev].CellSizeArray(); - const auto& problo = geom[lev].ProbLoArray(); - const auto& probhi = geom[lev].ProbHiArray(); - auto vel = velocity(lev).array(mfi); - auto rho = density(lev).array(mfi); - auto volfrac = vof(lev).array(mfi); - auto target_volfrac = m_ow_vof(lev).array(mfi); - auto target_vel = m_ow_vel(lev).array(mfi); - - const amrex::Real gen_length = wdata.gen_length; - const amrex::Real beach_length = wdata.beach_length; - const amrex::Real beach_length_factor = wdata.beach_length_factor; - const amrex::Real zsl = wdata.zsl; - const bool has_beach = wdata.has_beach; - const bool has_outprofile = wdata.has_outprofile; - - amrex::ParallelFor( - gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real x = amrex::min( - amrex::max(problo[0] + (i + 0.5) * dx[0], problo[0]), - probhi[0]); - const amrex::Real z = amrex::min( - amrex::max(problo[2] + (k + 0.5) * dx[2], problo[2]), - probhi[2]); - - // Generation region - if (x <= problo[0] + gen_length) { - const amrex::Real Gamma = - utils::gamma_generate(x - problo[0], gen_length); - // Get bounded new vof, incorporate with increment + const auto& dx = geom[lev].CellSizeArray(); + const auto& problo = geom[lev].ProbLoArray(); + const auto& probhi = geom[lev].ProbHiArray(); + auto vel_arrs = velocity(lev).arrays(); + auto rho_arrs = density(lev).arrays(); + auto volfrac_arrs = vof(lev).arrays(); + const auto target_volfrac_arrs = m_ow_vof(lev).const_arrays(); + const auto target_vel_arrs = m_ow_vel(lev).const_arrays(); + + const amrex::Real gen_length = wdata.gen_length; + const amrex::Real beach_length = wdata.beach_length; + const amrex::Real beach_length_factor = wdata.beach_length_factor; + const amrex::Real zsl = wdata.zsl; + const bool has_beach = wdata.has_beach; + const bool has_outprofile = wdata.has_outprofile; + + amrex::ParallelFor( + velocity(lev), amrex::IntVect(2), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real x = amrex::min( + amrex::max(problo[0] + (i + 0.5) * dx[0], problo[0]), + probhi[0]); + const amrex::Real z = amrex::min( + amrex::max(problo[2] + (k + 0.5) * dx[2], problo[2]), + probhi[2]); + + auto vel = vel_arrs[nbx]; + auto rho = rho_arrs[nbx]; + auto volfrac = volfrac_arrs[nbx]; + const auto target_volfrac = target_volfrac_arrs[nbx]; + const auto target_vel = target_vel_arrs[nbx]; + + // Generation region + if (x <= problo[0] + gen_length) { + const amrex::Real Gamma = + utils::gamma_generate(x - problo[0], gen_length); + // Get bounded new vof, incorporate with increment + amrex::Real new_vof = utils::combine_linear( + Gamma, target_volfrac(i, j, k), volfrac(i, j, k)); + 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)); + 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 = + utils::combine_linear( + Gamma, target_vel(i, j, k, n), 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_; + } + } + // 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) { + // Get bounded new vof, save increment + amrex::Real new_vof = utils::combine_linear( + Gamma, utils::free_surface_to_vof(zsl, z, dx[2]), + volfrac(i, j, k)); + 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) = new_vof; + // Conserve momentum when density changes + amrex::Real rho_ = rho1 * volfrac(i, j, k) + + rho2 * (1.0 - volfrac(i, j, k)); + // Target solution in liquid is vel = 0, assume + // added liquid already has target velocity + for (int n = 0; n < vel.ncomp; ++n) { + vel(i, j, k, n) = + (rho1 * (volfrac(i, j, k) * Gamma - + amrex::max(0.0, dvf)) + + rho2 * (1. - volfrac(i, j, k))) * + vel(i, j, k, n) / rho_; + } + } + // Forcing to wave profile instead + if (has_outprofile) { + // Same steps as in wave generation region amrex::Real new_vof = utils::combine_linear( Gamma, target_volfrac(i, j, k), volfrac(i, j, k)); 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 + 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)); 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 = utils::combine_linear( Gamma, target_vel(i, j, k, n), 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_; } } - // 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) { - // Get bounded new vof, save increment - amrex::Real new_vof = utils::combine_linear( - Gamma, - utils::free_surface_to_vof(zsl, z, dx[2]), - volfrac(i, j, k)); - 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) = new_vof; - // Conserve momentum when density changes - amrex::Real rho_ = rho1 * volfrac(i, j, k) + - rho2 * (1.0 - volfrac(i, j, k)); - // Target solution in liquid is vel = 0, assume - // added liquid already has target velocity - for (int n = 0; n < vel.ncomp; ++n) { - vel(i, j, k, n) = - (rho1 * (volfrac(i, j, k) * Gamma - - amrex::max(0.0, dvf)) + - rho2 * (1. - volfrac(i, j, k))) * - vel(i, j, k, n) / rho_; - } - } - // Forcing to wave profile instead - if (has_outprofile) { - // Same steps as in wave generation region - amrex::Real new_vof = utils::combine_linear( - Gamma, target_volfrac(i, j, k), - volfrac(i, j, k)); - 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)); - for (int n = 0; n < vel.ncomp; ++n) { - amrex::Real vel_liq = vel(i, j, k, n); - const amrex::Real dvel_liq = - utils::combine_linear( - Gamma, target_vel(i, j, k, n), - 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_; - } - } - } + } - // Make sure that density is updated before entering the - // solution - rho(i, j, k) = rho1 * volfrac(i, j, k) + - rho2 * (1. - volfrac(i, j, k)); - }); - } + // Make sure that density is updated before entering the + // solution + rho(i, j, k) = + rho1 * volfrac(i, j, k) + rho2 * (1. - volfrac(i, j, k)); + }); } // This helps for having periodic boundaries, but will need to be addressed // for the general case diff --git a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H index 9f84227123..0ff477842c 100644 --- a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H @@ -84,57 +84,54 @@ struct InitDataOp const auto& problo = geom.ProbLoArray(); const auto& probhi = geom.ProbHiArray(); const auto& dx = geom.CellSizeArray(); - for (amrex::MFIter mfi(m_levelset(level)); mfi.isValid(); ++mfi) { - const auto& phi = m_levelset(level).array(mfi); - const auto& vel = velocity(level).array(mfi); - - const amrex::Real wave_height = wdata.wave_height; - const amrex::Real wave_length = wdata.wave_length; - const amrex::Real water_depth = wdata.water_depth; - const amrex::Real zero_sea_level = wdata.zsl; - const amrex::Real gen_length = wdata.gen_length; - const amrex::Real beach_length = wdata.beach_length; - const amrex::Real g = wdata.g; - const int order = wdata.order; - const bool has_beach = wdata.has_beach; - const bool init_wave_field = wdata.init_wave_field; - const auto& gbx3 = mfi.growntilebox(3); - - amrex::ParallelFor( - gbx3, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; - const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; - - // Wave profile - amrex::Real eta_w{0.0}, u_w{0.0}, v_w{0.0}, w_w{0.0}; - relaxation_zones::stokes_waves( - order, wave_length, water_depth, wave_height, - zero_sea_level, g, x, z, 0.0, eta_w, u_w, v_w, w_w); - const utils::WaveVec wave_sol{u_w, v_w, w_w, eta_w}; - - // Quiescent profile - const utils::WaveVec quiescent{ - 0.0, 0.0, 0.0, zero_sea_level}; - - // Specify initial state for each region of domain - const auto bulk = init_wave_field ? wave_sol : quiescent; - const auto outlet = has_beach ? quiescent : wave_sol; - - const auto local_profile = utils::harmonize_profiles_1d( - x, problo[0], gen_length, probhi[0], beach_length, - wave_sol, bulk, outlet); - - phi(i, j, k) = local_profile[3] - z; - const amrex::Real cell_length_2D = - std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi(i, j, k) + cell_length_2D >= 0) { - vel(i, j, k, 0) = local_profile[0]; - vel(i, j, k, 1) = local_profile[1]; - vel(i, j, k, 2) = local_profile[2]; - } - }); - } + const auto& phi = m_levelset(level).arrays(); + const auto& vel = velocity(level).arrays(); + + const amrex::Real wave_height = wdata.wave_height; + const amrex::Real wave_length = wdata.wave_length; + const amrex::Real water_depth = wdata.water_depth; + const amrex::Real zero_sea_level = wdata.zsl; + const amrex::Real gen_length = wdata.gen_length; + const amrex::Real beach_length = wdata.beach_length; + const amrex::Real g = wdata.g; + const int order = wdata.order; + const bool has_beach = wdata.has_beach; + const bool init_wave_field = wdata.init_wave_field; + + amrex::ParallelFor( + velocity(level), amrex::IntVect(3), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + + // Wave profile + amrex::Real eta_w{0.0}, u_w{0.0}, v_w{0.0}, w_w{0.0}; + relaxation_zones::stokes_waves( + order, wave_length, water_depth, wave_height, + zero_sea_level, g, x, z, 0.0, eta_w, u_w, v_w, w_w); + const utils::WaveVec wave_sol{u_w, v_w, w_w, eta_w}; + + // Quiescent profile + const utils::WaveVec quiescent{0.0, 0.0, 0.0, zero_sea_level}; + + // Specify initial state for each region of domain + const auto bulk = init_wave_field ? wave_sol : quiescent; + const auto outlet = has_beach ? quiescent : wave_sol; + + const auto local_profile = utils::harmonize_profiles_1d( + x, problo[0], gen_length, probhi[0], beach_length, wave_sol, + bulk, outlet); + + phi[nbx](i, j, k) = local_profile[3] - z; + const amrex::Real cell_length_2D = + std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); + if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + vel[nbx](i, j, k, 0) = local_profile[0]; + vel[nbx](i, j, k, 1) = local_profile[1]; + vel[nbx](i, j, k, 2) = local_profile[2]; + } + }); } }; @@ -160,44 +157,42 @@ struct UpdateRelaxZonesOp const auto& problo = geom[lev].ProbLoArray(); const auto& dx = geom[lev].CellSizeArray(); - for (amrex::MFIter mfi(m_ow_levelset(lev)); mfi.isValid(); ++mfi) { - const auto& phi = m_ow_levelset(lev).array(mfi); - const auto& vel = m_ow_velocity(lev).array(mfi); - - const amrex::Real wave_height = wdata.wave_height; - const amrex::Real wave_length = wdata.wave_length; - const amrex::Real water_depth = wdata.water_depth; - const amrex::Real zero_sea_level = wdata.zsl; - const amrex::Real g = wdata.g; - const int order = wdata.order; - - const auto& gbx = mfi.growntilebox(); - amrex::ParallelFor( - gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; - const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; - - amrex::Real eta{0.0}, u_w{0.0}, v_w{0.0}, w_w{0.0}; - - relaxation_zones::stokes_waves( - order, wave_length, water_depth, wave_height, - zero_sea_level, g, x, z, time, eta, u_w, v_w, w_w); - - phi(i, j, k) = eta - z; - const amrex::Real cell_length_2D = - std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi(i, j, k) + cell_length_2D >= 0) { - // Wave velocity within a cell of interface - vel(i, j, k, 0) = u_w; - vel(i, j, k, 1) = v_w; - vel(i, j, k, 2) = w_w; - } else { - vel(i, j, k, 0) = 0.; - vel(i, j, k, 1) = 0.; - vel(i, j, k, 2) = 0.; - } - }); - } + const auto& phi = m_ow_levelset(lev).arrays(); + const auto& vel = m_ow_velocity(lev).arrays(); + + const amrex::Real wave_height = wdata.wave_height; + const amrex::Real wave_length = wdata.wave_length; + const amrex::Real water_depth = wdata.water_depth; + const amrex::Real zero_sea_level = wdata.zsl; + const amrex::Real g = wdata.g; + const int order = wdata.order; + + amrex::ParallelFor( + m_ow_velocity(lev), amrex::IntVect(3), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + + amrex::Real eta{0.0}, u_w{0.0}, v_w{0.0}, w_w{0.0}; + + relaxation_zones::stokes_waves( + order, wave_length, water_depth, wave_height, + zero_sea_level, g, x, z, time, eta, u_w, v_w, w_w); + + phi[nbx](i, j, k) = eta - z; + const amrex::Real cell_length_2D = + std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); + if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + // Wave velocity within a cell of interface + vel[nbx](i, j, k, 0) = u_w; + vel[nbx](i, j, k, 1) = v_w; + vel[nbx](i, j, k, 2) = w_w; + } else { + vel[nbx](i, j, k, 0) = 0.; + vel[nbx](i, j, k, 1) = 0.; + vel[nbx](i, j, k, 2) = 0.; + } + }); } } }; diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index 8fb10515c7..259565d842 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -82,23 +82,20 @@ void postprocess_velocity_mfab_liquid( amrex::MultiFab& lvs_mfab, const amrex::GpuArray dx) { - for (amrex::MFIter mfi(vel_mfab); mfi.isValid(); ++mfi) { - auto vel = vel_mfab.array(mfi); - auto phi = lvs_mfab.const_array(mfi); - const auto& vbx = mfi.validbox(); - const auto& gbx = grow(vbx, 3); - amrex::ParallelFor( - gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // Set velocity to zero if no liquid present - const amrex::Real cell_length_2D = - std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi(i, j, k) + cell_length_2D >= 0) { - vel(i, j, k, 0) = 0.0; - vel(i, j, k, 1) = 0.0; - vel(i, j, k, 2) = 0.0; - } - }); - } + auto vel = vel_mfab.arrays(); + const auto phi = lvs_mfab.const_arrays(); + amrex::ParallelFor( + vel_mfab, amrex::IntVect(3), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + // Set velocity to zero if no liquid present + const amrex::Real cell_length_2D = + std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); + if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + vel[nbx](i, j, k, 0) = 0.0; + vel[nbx](i, j, k, 1) = 0.0; + vel[nbx](i, j, k, 2) = 0.0; + } + }); } void postprocess_velocity_field_liquid( @@ -449,51 +446,48 @@ struct InitDataOp ow_velocity(level), ow_levelset(level), dx); // Populate flow fields according to intended forcing and init setup - for (amrex::MFIter mfi(levelset(level)); mfi.isValid(); ++mfi) { - const auto& ow_phi = ow_levelset(level).array(mfi); - const auto& ow_vel = ow_velocity(level).array(mfi); - const auto& phi = levelset(level).array(mfi); - const auto& vel = velocity(level).array(mfi); + const auto& ow_phi = ow_levelset(level).const_arrays(); + const auto& ow_vel = ow_velocity(level).const_arrays(); + const auto& phi = levelset(level).arrays(); + const auto& vel = velocity(level).arrays(); - const amrex::Real gen_length = wdata.gen_length; - const amrex::Real beach_length = wdata.beach_length; - const amrex::Real zero_sea_level = wdata.zsl; + const amrex::Real gen_length = wdata.gen_length; + const amrex::Real beach_length = wdata.beach_length; + const amrex::Real zero_sea_level = wdata.zsl; - const bool has_beach = wdata.has_beach; - const bool init_wave_field = wdata.init_wave_field; - const auto& gbx3 = mfi.growntilebox(3); + const bool has_beach = wdata.has_beach; + const bool init_wave_field = wdata.init_wave_field; - amrex::ParallelFor( - gbx3, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; - const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; - - // Wave profile - const utils::WaveVec wave_sol{ - ow_vel(i, j, k, 0), ow_vel(i, j, k, 1), - ow_vel(i, j, k, 2), ow_phi(i, j, k) + z}; - // Quiescent profile - const utils::WaveVec quiescent{ - 0.0, 0.0, 0.0, zero_sea_level}; - - // Specify initial state for each region of domain - const auto bulk = init_wave_field ? wave_sol : quiescent; - const auto outlet = has_beach ? quiescent : wave_sol; - - const auto local_profile = utils::harmonize_profiles_1d( - x, problo[0], gen_length, probhi[0], beach_length, - wave_sol, bulk, outlet); - - phi(i, j, k) = local_profile[3] - z; - const amrex::Real cell_length_2D = - std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi(i, j, k) + cell_length_2D >= 0) { - vel(i, j, k, 0) = local_profile[0]; - vel(i, j, k, 1) = local_profile[1]; - vel(i, j, k, 2) = local_profile[2]; - } - }); - } + amrex::ParallelFor( + velocity(level), amrex::IntVect(3), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; + const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; + + // Wave profile + const utils::WaveVec wave_sol{ + ow_vel[nbx](i, j, k, 0), ow_vel[nbx](i, j, k, 1), + ow_vel[nbx](i, j, k, 2), ow_phi[nbx](i, j, k) + z}; + // Quiescent profile + const utils::WaveVec quiescent{0.0, 0.0, 0.0, zero_sea_level}; + + // Specify initial state for each region of domain + const auto bulk = init_wave_field ? wave_sol : quiescent; + const auto outlet = has_beach ? quiescent : wave_sol; + + const auto local_profile = utils::harmonize_profiles_1d( + x, problo[0], gen_length, probhi[0], beach_length, wave_sol, + bulk, outlet); + + phi[nbx](i, j, k) = local_profile[3] - z; + const amrex::Real cell_length_2D = + std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); + if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + vel[nbx](i, j, k, 0) = local_profile[0]; + vel[nbx](i, j, k, 1) = local_profile[1]; + vel[nbx](i, j, k, 2) = local_profile[2]; + } + }); // Start w2a fields at 0, some areas will not be modified auto& w2a_levelset = sim.repo().get_field("w2a_levelset"); @@ -676,31 +670,29 @@ struct UpdateRelaxZonesOp // Temporally interpolate at every timestep to get target solution for (int lev = 0; lev < nlevels; ++lev) { - for (amrex::MFIter mfi(m_ow_levelset(lev)); mfi.isValid(); ++mfi) { - auto phi = m_ow_levelset(lev).array(mfi); - auto vel = m_ow_velocity(lev).array(mfi); - auto W2A_phi = w2a_levelset(lev).array(mfi); - auto W2A_vel = w2a_velocity(lev).array(mfi); - - const amrex::Real W2A_t = wdata.t; - const auto& gbx = mfi.growntilebox(3); - amrex::ParallelFor( - gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // Interpolate temporally every time - phi(i, j, k) += (W2A_phi(i, j, k) - phi(i, j, k)) * - (newtime - t_last) / - (W2A_t - t_last + 1e-16); - vel(i, j, k, 0) += - (W2A_vel(i, j, k, 0) - vel(i, j, k, 0)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); - vel(i, j, k, 1) += - (W2A_vel(i, j, k, 1) - vel(i, j, k, 1)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); - vel(i, j, k, 2) += - (W2A_vel(i, j, k, 2) - vel(i, j, k, 2)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); - }); - } + auto phi = m_ow_levelset(lev).arrays(); + auto vel = m_ow_velocity(lev).arrays(); + auto W2A_phi = w2a_levelset(lev).arrays(); + auto W2A_vel = w2a_velocity(lev).arrays(); + + const amrex::Real W2A_t = wdata.t; + amrex::ParallelFor( + m_ow_levelset(lev), amrex::IntVect(3), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + // Interpolate temporally every time + phi[nbx](i, j, k) += + (W2A_phi[nbx](i, j, k) - phi[nbx](i, j, k)) * + (newtime - t_last) / (W2A_t - t_last + 1e-16); + vel[nbx](i, j, k, 0) += + (W2A_vel[nbx](i, j, k, 0) - vel[nbx](i, j, k, 0)) * + (newtime - t_last) / (W2A_t - t_last + 1e-16); + vel[nbx](i, j, k, 1) += + (W2A_vel[nbx](i, j, k, 1) - vel[nbx](i, j, k, 1)) * + (newtime - t_last) / (W2A_t - t_last + 1e-16); + vel[nbx](i, j, k, 2) += + (W2A_vel[nbx](i, j, k, 2) - vel[nbx](i, j, k, 2)) * + (newtime - t_last) / (W2A_t - t_last + 1e-16); + }); } #else amrex::ignore_unused(data); From e023a6f8f9c780c6acaa2d18e4efb9d20c2c7a4b Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 11:13:17 -0700 Subject: [PATCH 08/31] changes to use ow fields in terrain drag --- amr-wind/physics/TerrainDrag.H | 10 ++ amr-wind/physics/TerrainDrag.cpp | 180 ++++++++++++++++++------------- 2 files changed, 117 insertions(+), 73 deletions(-) diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index 6497be1d53..344334a8fc 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -34,6 +34,8 @@ public: void post_advance_work() override {} + std::string wave_velocity_field_name() { return m_wave_velocity_name; } + private: CFDSim& m_sim; const FieldRepo& m_repo; @@ -54,6 +56,14 @@ private: amrex::Vector m_xrough; amrex::Vector m_yrough; amrex::Vector m_z0rough; + + //! Terrain Drag Force Term + bool m_terrain_is_waves{false}; + Field* m_wave_volume_fraction{nullptr}; + Field* m_wave_negative_elevation{nullptr}; + const std::string m_wave_volume_fraction_name{"ow_vof"}; + const std::string m_wave_negative_elevation_name{"ow_levelset"}; + const std::string m_wave_velocity_name{"ow_velocity"}; }; } // namespace amr_wind::terraindrag diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 2b6b661757..78ff7a6481 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -22,22 +22,31 @@ TerrainDrag::TerrainDrag(CFDSim& sim) , m_terrainz0(sim.repo().declare_field("terrainz0", 1, 1, 1)) , m_terrain_height(sim.repo().declare_field("terrain_height", 1, 1, 1)) { - std::string terrain_file("terrain.amrwind"); - std::string roughness_file("terrain.roughness"); - amrex::ParmParse pp(identifier()); - pp.query("terrain_file", terrain_file); - pp.query("roughness_file", roughness_file); - - ioutils::read_flat_grid_file( - terrain_file, m_xterrain, m_yterrain, m_zterrain); - - // No checks for the file as it is optional currently - std::ifstream file(roughness_file, std::ios::in); - if (file.good()) { + + m_terrain_is_waves = sim.physics_manager().contains("OceanWaves"); + + if (!m_terrain_is_waves) { + std::string terrain_file("terrain.amrwind"); + std::string roughness_file("terrain.roughness"); + amrex::ParmParse pp(identifier()); + pp.query("terrain_file", terrain_file); + pp.query("roughness_file", roughness_file); + ioutils::read_flat_grid_file( - roughness_file, m_xrough, m_yrough, m_z0rough); + terrain_file, m_xterrain, m_yterrain, m_zterrain); + + // No checks for the file as it is optional currently + std::ifstream file(roughness_file, std::ios::in); + if (file.good()) { + ioutils::read_flat_grid_file( + roughness_file, m_xrough, m_yrough, m_z0rough); + } + file.close(); + } else { + m_wave_volume_fraction = &m_repo.get_field(m_wave_volume_fraction_name); + m_wave_negative_elevation = + &m_repo.get_field(m_wave_negative_elevation_name); } - file.close(); m_sim.io_manager().register_output_int_var("terrain_drag"); m_sim.io_manager().register_output_int_var("terrain_blank"); @@ -59,69 +68,94 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) auto& terrainz0 = m_terrainz0(level); auto& terrain_height = m_terrain_height(level); auto& drag = m_terrain_drag(level); - const auto xterrain_size = m_xterrain.size(); - const auto yterrain_size = m_yterrain.size(); - const auto zterrain_size = m_zterrain.size(); - amrex::Gpu::DeviceVector device_xterrain(xterrain_size); - amrex::Gpu::DeviceVector device_yterrain(yterrain_size); - amrex::Gpu::DeviceVector device_zterrain(zterrain_size); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), - device_xterrain.begin()); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), - device_yterrain.begin()); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), - device_zterrain.begin()); - const auto* xterrain_ptr = device_xterrain.data(); - const auto* yterrain_ptr = device_yterrain.data(); - const auto* zterrain_ptr = device_zterrain.data(); - // Copy Roughness to gpu - const auto xrough_size = m_xrough.size(); - const auto yrough_size = m_yrough.size(); - const auto z0rough_size = m_z0rough.size(); - amrex::Gpu::DeviceVector device_xrough(xrough_size); - amrex::Gpu::DeviceVector device_yrough(yrough_size); - amrex::Gpu::DeviceVector device_z0rough(z0rough_size); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_xrough.begin(), m_xrough.end(), - device_xrough.begin()); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_yrough.begin(), m_yrough.end(), - device_yrough.begin()); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_z0rough.begin(), m_z0rough.end(), - device_z0rough.begin()); - const auto* xrough_ptr = device_xrough.data(); - const auto* yrough_ptr = device_yrough.data(); - const auto* z0rough_ptr = device_z0rough.data(); + auto levelBlanking = blanking.arrays(); auto levelDrag = drag.arrays(); auto levelz0 = terrainz0.arrays(); auto levelheight = terrain_height.arrays(); - amrex::ParallelFor( - blanking, m_terrain_blank.num_grow(), - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; - const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; - const amrex::Real terrainHt = interp::bilinear( - xterrain_ptr, xterrain_ptr + xterrain_size, yterrain_ptr, - yterrain_ptr + yterrain_size, zterrain_ptr, x, y); - levelBlanking[nbx](i, j, k, 0) = - static_cast((z <= terrainHt) && (z > prob_lo[2])); - levelheight[nbx](i, j, k, 0) = - std::max(std::abs(z - terrainHt), 0.5 * dx[2]); - - amrex::Real roughz0 = 0.1; - if (xrough_size > 0) { - roughz0 = interp::bilinear( - xrough_ptr, xrough_ptr + xrough_size, yrough_ptr, - yrough_ptr + yrough_size, z0rough_ptr, x, y); - } - levelz0[nbx](i, j, k, 0) = roughz0; - }); + + if (!m_terrain_is_waves) { + const auto xterrain_size = m_xterrain.size(); + const auto yterrain_size = m_yterrain.size(); + const auto zterrain_size = m_zterrain.size(); + amrex::Gpu::DeviceVector device_xterrain(xterrain_size); + amrex::Gpu::DeviceVector device_yterrain(yterrain_size); + amrex::Gpu::DeviceVector device_zterrain(zterrain_size); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), + device_xterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), + device_yterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), + device_zterrain.begin()); + const auto* xterrain_ptr = device_xterrain.data(); + const auto* yterrain_ptr = device_yterrain.data(); + const auto* zterrain_ptr = device_zterrain.data(); + // Copy Roughness to gpu + const auto xrough_size = m_xrough.size(); + const auto yrough_size = m_yrough.size(); + const auto z0rough_size = m_z0rough.size(); + amrex::Gpu::DeviceVector device_xrough(xrough_size); + amrex::Gpu::DeviceVector device_yrough(yrough_size); + amrex::Gpu::DeviceVector device_z0rough(z0rough_size); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_xrough.begin(), m_xrough.end(), + device_xrough.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_yrough.begin(), m_yrough.end(), + device_yrough.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_z0rough.begin(), m_z0rough.end(), + device_z0rough.begin()); + const auto* xrough_ptr = device_xrough.data(); + const auto* yrough_ptr = device_yrough.data(); + const auto* z0rough_ptr = device_z0rough.data(); + + amrex::ParallelFor( + blanking, m_terrain_blank.num_grow(), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + const amrex::Real terrainHt = interp::bilinear( + xterrain_ptr, xterrain_ptr + xterrain_size, yterrain_ptr, + yterrain_ptr + yterrain_size, zterrain_ptr, x, y); + levelBlanking[nbx](i, j, k, 0) = + static_cast((z <= terrainHt) && (z > prob_lo[2])); + levelheight[nbx](i, j, k, 0) = + std::max(std::abs(z - terrainHt), 0.5 * dx[2]); + + amrex::Real roughz0 = 0.1; + if (xrough_size > 0) { + roughz0 = interp::bilinear( + xrough_ptr, xrough_ptr + xrough_size, yrough_ptr, + yrough_ptr + yrough_size, z0rough_ptr, x, y); + } + levelz0[nbx](i, j, k, 0) = roughz0; + }); + } else { + const auto negative_wave_elevation = + (*m_wave_negative_elevation)(level).const_arrays(); + const auto wave_vol_frac = + (*m_wave_volume_fraction)(level).const_arrays(); + + // Get terrain blanking from ocean waves fields + amrex::ParallelFor( + blanking, m_terrain_blank.num_grow(), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + levelBlanking[nbx](i, j, k, 0) = static_cast( + (wave_vol_frac[nbx](i, j, k) >= 0.5) && (z > prob_lo[2])); + levelheight[nbx](i, j, k, 0) = std::max( + std::abs(negative_wave_elevation[nbx](i, j, k)), + 0.5 * dx[2]); + + // Set uniform roughness (ocean surface) + levelz0[nbx](i, j, k, 0) = 1e-4; + }); + } amrex::ParallelFor( blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { if ((levelBlanking[nbx](i, j, k, 0) == 0) && (k > 0) && From ca7b1ca7a31362197969831bd49f2036da342d6d Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 11:14:19 -0700 Subject: [PATCH 09/31] create modes for ocean waves --- amr-wind/ocean_waves/OceanWaves.H | 5 +++++ amr-wind/ocean_waves/OceanWaves.cpp | 23 ++++++++++++++++++----- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/amr-wind/ocean_waves/OceanWaves.H b/amr-wind/ocean_waves/OceanWaves.H index a3850589e7..de32157347 100644 --- a/amr-wind/ocean_waves/OceanWaves.H +++ b/amr-wind/ocean_waves/OceanWaves.H @@ -68,6 +68,11 @@ private: //! Ocean waves pressure gradient // Field& m_ow_pressure; + + //! Multiphase mode forces the velocity and vof in the domain + //! When off, single-phase mode simply provides the target variables, and + //! forcing is handled in other parts of the code + bool m_multiphase_mode{true}; }; } // namespace ocean_waves diff --git a/amr-wind/ocean_waves/OceanWaves.cpp b/amr-wind/ocean_waves/OceanWaves.cpp index 86372b2c7f..47778089e9 100644 --- a/amr-wind/ocean_waves/OceanWaves.cpp +++ b/amr-wind/ocean_waves/OceanWaves.cpp @@ -17,8 +17,14 @@ OceanWaves::OceanWaves(CFDSim& sim) , m_ow_velocity( sim.repo().declare_field("ow_velocity", AMREX_SPACEDIM, 3, 1)) { - if (!sim.physics_manager().contains("MultiPhase")) { - amrex::Abort("OceanWaves requires Multiphase physics to be active"); + if ((!sim.physics_manager().contains("MultiPhase") || + sim.physics_manager().contains("TerrainDrag"))) { + amrex::Abort( + "OceanWaves requires MultiPhase or TerrainDrag physics to be " + "active"); + } + if (!sim.physics_manager.contains("MultiPhase")) { + m_multiphase_mode = false; } m_ow_levelset.set_default_fillpatch_bc(sim.time()); m_ow_vof.set_default_fillpatch_bc(sim.time()); @@ -54,13 +60,15 @@ void OceanWaves::pre_init_actions() void OceanWaves::initialize_fields(int level, const amrex::Geometry& geom) { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::initialize_fields"); - m_owm->init_waves(level, geom); + m_owm->init_waves(level, geom, m_multiphase_mode); } void OceanWaves::post_init_actions() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::post_init_actions"); - relaxation_zones(); + if (m_multiphase_mode) { + relaxation_zones(); + } } void OceanWaves::post_regrid_actions() @@ -72,12 +80,17 @@ void OceanWaves::post_regrid_actions() void OceanWaves::pre_advance_work() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_advance_work"); + if (!m_multiphase_mode) { + m_owm->update_relax_zones(); + } } void OceanWaves::post_advance_work() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::post_init_actions"); - relaxation_zones(); + if (m_multiphase_mode) { + relaxation_zones(); + } } /** Update ocean waves relaxation zones From dc2870f1ddf2ac35eeed393a3cbcdb1f28d9d956 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 11:15:04 -0700 Subject: [PATCH 10/31] correct conditional --- amr-wind/ocean_waves/OceanWaves.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/ocean_waves/OceanWaves.cpp b/amr-wind/ocean_waves/OceanWaves.cpp index 47778089e9..27961cfda5 100644 --- a/amr-wind/ocean_waves/OceanWaves.cpp +++ b/amr-wind/ocean_waves/OceanWaves.cpp @@ -17,7 +17,7 @@ OceanWaves::OceanWaves(CFDSim& sim) , m_ow_velocity( sim.repo().declare_field("ow_velocity", AMREX_SPACEDIM, 3, 1)) { - if ((!sim.physics_manager().contains("MultiPhase") || + if (!(sim.physics_manager().contains("MultiPhase") || sim.physics_manager().contains("TerrainDrag"))) { amrex::Abort( "OceanWaves requires MultiPhase or TerrainDrag physics to be " From 84a891eee349551dbbbb57d7e2f3d043868cf51a Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 11:58:35 -0700 Subject: [PATCH 11/31] moving toward ocean waves single-phase compatibility --- amr-wind/ocean_waves/OceanWaves.cpp | 6 +- amr-wind/ocean_waves/OceanWavesModel.H | 12 +-- .../relaxation_zones/linear_waves_ops.H | 51 ++++++++----- .../relaxation_zones/stokes_waves_ops.H | 43 +++++++---- .../relaxation_zones/waves2amr_ops.H | 73 +++++++++++-------- 5 files changed, 111 insertions(+), 74 deletions(-) diff --git a/amr-wind/ocean_waves/OceanWaves.cpp b/amr-wind/ocean_waves/OceanWaves.cpp index 27961cfda5..45e176c9e8 100644 --- a/amr-wind/ocean_waves/OceanWaves.cpp +++ b/amr-wind/ocean_waves/OceanWaves.cpp @@ -23,7 +23,7 @@ OceanWaves::OceanWaves(CFDSim& sim) "OceanWaves requires MultiPhase or TerrainDrag physics to be " "active"); } - if (!sim.physics_manager.contains("MultiPhase")) { + if (!sim.physics_manager().contains("MultiPhase")) { m_multiphase_mode = false; } m_ow_levelset.set_default_fillpatch_bc(sim.time()); @@ -81,7 +81,7 @@ void OceanWaves::pre_advance_work() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_advance_work"); if (!m_multiphase_mode) { - m_owm->update_relax_zones(); + m_owm->update_relax_zones(true); } } @@ -99,7 +99,7 @@ void OceanWaves::post_advance_work() void OceanWaves::relaxation_zones() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_relaxation_zones"); - m_owm->update_relax_zones(); + m_owm->update_relax_zones(false); m_owm->apply_relax_zones(); m_owm->reset_regrid_flag(); } diff --git a/amr-wind/ocean_waves/OceanWavesModel.H b/amr-wind/ocean_waves/OceanWavesModel.H index 86fa87bdbf..45e5c4dead 100644 --- a/amr-wind/ocean_waves/OceanWavesModel.H +++ b/amr-wind/ocean_waves/OceanWavesModel.H @@ -32,9 +32,9 @@ public: virtual void read_inputs(const ::amr_wind::utils::MultiParser&) = 0; - virtual void init_waves(int, const amrex::Geometry&) = 0; + virtual void init_waves(int, const amrex::Geometry&, bool) = 0; - virtual void update_relax_zones() = 0; + virtual void update_relax_zones(bool) = 0; virtual void apply_relax_zones() = 0; @@ -91,9 +91,9 @@ public: m_out_op.read_io_options(pp); } - void update_relax_zones() override + void update_relax_zones(bool for_forcing_term) override { - ops::UpdateRelaxZonesOp()(m_data); + ops::UpdateRelaxZonesOp()(m_data, for_forcing_term); } void apply_relax_zones() override @@ -112,9 +112,9 @@ public: void write_outputs() override { m_out_op.write_outputs(); } - void init_waves(int level, const amrex::Geometry& geom) override + void init_waves(int level, const amrex::Geometry& geom, bool multiphase_mode) override { - ops::InitDataOp()(m_data, level, geom); + ops::InitDataOp()(m_data, level, geom, multiphase_mode); } }; diff --git a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H index 8b624d3a5b..41a01d8f43 100644 --- a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H @@ -31,36 +31,42 @@ template <> struct InitDataOp { void operator()( - LinearWaves::DataType& data, int level, const amrex::Geometry& geom) + LinearWaves::DataType& data, + int level, + const amrex::Geometry& geom, + bool multiphase_mode) { const auto& wdata = data.meta(); auto& sim = data.sim(); + Field* levelset{nullptr}; + if (multiphase_mode) { + levelset = &sim.repo().get_field("levelset"); + } // cppcheck-suppress constVariableReference - auto& m_levelset = sim.repo().get_field("levelset"); - // cppcheck-suppress constVariableReference - auto& m_velocity = sim.repo().get_field("velocity"); + auto& velocity = sim.repo().get_field("velocity"); const auto& problo = geom.ProbLoArray(); const auto& probhi = geom.ProbHiArray(); const auto& dx = geom.CellSizeArray(); - const auto& phi = m_levelset(level).arrays(); - const auto& vel = m_velocity(level).arrays(); + const auto& vel = velocity(level).arrays(); + const auto& phi_arrs = multiphase_mode ? (*levelset)(level).arrays() + : amrex::MultiArray4(); const amrex::Real zero_sea_level = wdata.zsl; const amrex::Real gen_length = wdata.gen_length; const amrex::Real beach_length = wdata.beach_length; const amrex::Real g = wdata.g; - const bool has_beach = wdata.has_beach; - const bool init_wave_field = wdata.init_wave_field; + const bool has_beach = wdata.has_beach && multiphase_mode; + const bool init_wave_field = wdata.init_wave_field || !multiphase_mode; const amrex::Real wave_height = wdata.wave_height; const amrex::Real wave_length = wdata.wave_length; const amrex::Real water_depth = wdata.water_depth; amrex::ParallelFor( - m_velocity(level), amrex::IntVect(3), + velocity(level), amrex::IntVect(3), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; @@ -98,32 +104,39 @@ struct InitDataOp x, problo[0], gen_length, probhi[0], beach_length, wave_sol, bulk, outlet); - phi[nbx](i, j, k) = local_profile[3] - z; + const amrex::Real phi = local_profile[3] - z; const amrex::Real cell_length_2D = std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + if (phi + cell_length_2D >= 0) { vel[nbx](i, j, k, 0) = local_profile[0]; vel[nbx](i, j, k, 1) = local_profile[1]; vel[nbx](i, j, k, 2) = local_profile[2]; } + if (multiphase_mode) { + phi_arrs[nbx](i, j, k) = phi; + } }); } -}; +}; // namespace amr_wind::ocean_waves::ops template <> struct UpdateRelaxZonesOp { - void operator()(LinearWaves::DataType& data) + void operator()(LinearWaves::DataType& data, bool for_forcing_term) { const auto& wdata = data.meta(); auto& sim = data.sim(); - const auto& time = sim.time().new_time(); + const amrex::Real time = + sim.time().new_time() - + (for_forcing_term + ? 0.5 * (sim.time().new_time() - sim.time().current_time()) + : 0.0); // cppcheck-suppress constVariableReference - auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); + auto& ow_levelset = sim.repo().get_field("ow_levelset"); // cppcheck-suppress constVariableReference - auto& m_ow_velocity = sim.repo().get_field("ow_velocity"); + auto& ow_velocity = sim.repo().get_field("ow_velocity"); auto nlevels = sim.repo().num_active_levels(); auto geom = sim.mesh().Geom(); @@ -132,8 +145,8 @@ struct UpdateRelaxZonesOp const auto& problo = geom[lev].ProbLoArray(); const auto& dx = geom[lev].CellSizeArray(); - const auto& phi = m_ow_levelset(lev).arrays(); - const auto& vel = m_ow_velocity(lev).arrays(); + const auto& phi = ow_levelset(lev).arrays(); + const auto& vel = ow_velocity(lev).arrays(); const amrex::Real wave_height = wdata.wave_height; const amrex::Real wave_length = wdata.wave_length; @@ -142,7 +155,7 @@ struct UpdateRelaxZonesOp const amrex::Real g = wdata.g; amrex::ParallelFor( - m_ow_velocity(lev), amrex::IntVect(3), + ow_velocity(lev), amrex::IntVect(3), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; diff --git a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H index 0ff477842c..2a5a42ecd0 100644 --- a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H @@ -72,21 +72,27 @@ template <> struct InitDataOp { void operator()( - StokesWaves::DataType& data, int level, const amrex::Geometry& geom) + StokesWaves::DataType& data, + int level, + const amrex::Geometry& geom, + bool multiphase_mode) { const auto& wdata = data.meta(); auto& sim = data.sim(); - // cppcheck-suppress constVariableReference - auto& m_levelset = sim.repo().get_field("levelset"); + Field* levelset{nullptr}; + if (multiphase_mode) { + levelset = &sim.repo().get_field("levelset"); + } // cppcheck-suppress constVariableReference auto& velocity = sim.repo().get_field("velocity"); const auto& problo = geom.ProbLoArray(); const auto& probhi = geom.ProbHiArray(); const auto& dx = geom.CellSizeArray(); - const auto& phi = m_levelset(level).arrays(); const auto& vel = velocity(level).arrays(); + const auto& phi_arrs = multiphase_mode ? (*levelset)(level).arrays() + : amrex::MultiArray4(); const amrex::Real wave_height = wdata.wave_height; const amrex::Real wave_length = wdata.wave_length; @@ -96,8 +102,8 @@ struct InitDataOp const amrex::Real beach_length = wdata.beach_length; const amrex::Real g = wdata.g; const int order = wdata.order; - const bool has_beach = wdata.has_beach; - const bool init_wave_field = wdata.init_wave_field; + const bool has_beach = wdata.has_beach && multiphase_mode; + const bool init_wave_field = wdata.init_wave_field || !multiphase_mode; amrex::ParallelFor( velocity(level), amrex::IntVect(3), @@ -123,14 +129,17 @@ struct InitDataOp x, problo[0], gen_length, probhi[0], beach_length, wave_sol, bulk, outlet); - phi[nbx](i, j, k) = local_profile[3] - z; + const amrex::Real phi = local_profile[3] - z; const amrex::Real cell_length_2D = std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + if (phi + cell_length_2D >= 0) { vel[nbx](i, j, k, 0) = local_profile[0]; vel[nbx](i, j, k, 1) = local_profile[1]; vel[nbx](i, j, k, 2) = local_profile[2]; } + if (multiphase_mode) { + phi_arrs[nbx](i, j, k) = phi; + } }); } }; @@ -138,17 +147,21 @@ struct InitDataOp template <> struct UpdateRelaxZonesOp { - void operator()(StokesWaves::DataType& data) + void operator()(StokesWaves::DataType& data, bool for_forcing_term) { const auto& wdata = data.meta(); auto& sim = data.sim(); - const auto& time = sim.time().new_time(); + const amrex::Real time = + sim.time().new_time() - + (for_forcing_term + ? 0.5 * (sim.time().new_time() - sim.time().current_time()) + : 0.0); // cppcheck-suppress constVariableReference - auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); + auto& ow_levelset = sim.repo().get_field("ow_levelset"); // cppcheck-suppress constVariableReference - auto& m_ow_velocity = sim.repo().get_field("ow_velocity"); + auto& ow_velocity = sim.repo().get_field("ow_velocity"); auto nlevels = sim.repo().num_active_levels(); auto geom = sim.mesh().Geom(); @@ -157,8 +170,8 @@ struct UpdateRelaxZonesOp const auto& problo = geom[lev].ProbLoArray(); const auto& dx = geom[lev].CellSizeArray(); - const auto& phi = m_ow_levelset(lev).arrays(); - const auto& vel = m_ow_velocity(lev).arrays(); + const auto& phi = ow_levelset(lev).arrays(); + const auto& vel = ow_velocity(lev).arrays(); const amrex::Real wave_height = wdata.wave_height; const amrex::Real wave_length = wdata.wave_length; @@ -168,7 +181,7 @@ struct UpdateRelaxZonesOp const int order = wdata.order; amrex::ParallelFor( - m_ow_velocity(lev), amrex::IntVect(3), + ow_velocity(lev), amrex::IntVect(3), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { const amrex::Real x = problo[0] + (i + 0.5) * dx[0]; const amrex::Real z = problo[2] + (k + 0.5) * dx[2]; diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index 259565d842..dedbf36943 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -23,7 +23,7 @@ int evaluate_read_resize( const int new_ntime, const amrex::Real wtinit, const amrex::Real wdt, - const amrex::Real newtime) + const amrex::Real time) { // Flag to indicate that data must be read twice for interpolation int double_data = 0; @@ -41,7 +41,7 @@ int evaluate_read_resize( // Sim time to go with recorded data wtime = new_ntime * wdt + wtinit; // If double reading is deemed necessary, check for convenience - if (double_data == 1 && std::abs(wtime - newtime) <= 1e-10) { + if (double_data == 1 && std::abs(wtime - time) <= 1e-10) { // Reading can be done just once, w2a fields replace ow fields double_data = 2; } @@ -54,7 +54,7 @@ int evaluate_read_resize( resize_flag = true; // Confirm that new time is not coincident with modes time step - if (std::abs(wtime - newtime) > 1e-10) { + if (std::abs(wtime - time) > 1e-10) { // Data must be read before and after wtime double_data = 1; } else { @@ -72,7 +72,7 @@ int evaluate_read_resize( // interpolation is ready to go } // Record latest time as 'last' for next timestep - t_last = newtime; + t_last = time; // Return flag regarding double reading return double_data; } @@ -406,7 +406,7 @@ struct InitDataOp void // cppcheck-suppress constParameterReference operator()( - W2AWaves::DataType & data, int level, const amrex::Geometry & geom) + W2AWaves::DataType & data, int level, const amrex::Geometry & geom, bool multiphase_mode) { #ifdef AMR_WIND_USE_W2A @@ -504,19 +504,23 @@ template <> struct UpdateRelaxZonesOp { // cppcheck-suppress constParameterReference - void operator()(W2AWaves::DataType& data) + void operator()(W2AWaves::DataType& data, bool for_forcing_term) { #ifdef AMR_WIND_USE_W2A auto& wdata = data.meta(); auto& sim = data.sim(); - // Nudge the solution toward where it should be - const amrex::Real newtime = sim.time().new_time(); + // Provide target solution based on intended part of the algorithm + const amrex::Real time = + sim.time().new_time() - + (for_forcing_term + ? 0.5 * (sim.time().new_time() - sim.time().current_time()) + : 0.0); // Update ow fields every time - auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); - auto& m_ow_velocity = sim.repo().get_field("ow_velocity"); + auto& ow_levelset = sim.repo().get_field("ow_levelset"); + auto& ow_velocity = sim.repo().get_field("ow_velocity"); // Update HOS fields when necessary auto& w2a_levelset = sim.repo().get_field("w2a_levelset"); auto& w2a_velocity = sim.repo().get_field("w2a_velocity"); @@ -531,10 +535,10 @@ struct UpdateRelaxZonesOp bool read_flag = false; // Check if time indicates reading must take place int new_ntime = - wdata.rmodes.time2step(newtime + wdata.t_winit, wdata.ntime); + wdata.rmodes.time2step(time + wdata.t_winit, wdata.ntime); int double_data = evaluate_read_resize( wdata.ntime, read_flag, wdata.resize_flag, wdata.t, wdata.t_last, - new_ntime, wdata.t_winit, wdata.dt_modes, newtime); + new_ntime, wdata.t_winit, wdata.dt_modes, time); // Check if stop time is exceeded, introduce offset to ntime if (read_flag) { // Need to only check when reading is occurring @@ -568,7 +572,7 @@ struct UpdateRelaxZonesOp // Get heights for this processor, check overlap in z flag_z = (interp_to_mfab::get_local_height_indices( - wdata.indvec, wdata.hvec, m_ow_velocity.vec_ptrs(), + wdata.indvec, wdata.hvec, ow_velocity.vec_ptrs(), geom) == 1); // No overlap from heights definitely means no interp @@ -576,17 +580,24 @@ struct UpdateRelaxZonesOp const int dir = 0; flag_xlo = (interp_to_mfab::check_lateral_overlap_lo( - wdata.gen_length, dir, m_ow_velocity.vec_ptrs(), + wdata.gen_length, dir, ow_velocity.vec_ptrs(), geom) == 1); // No overlap with gen region means no interp, unless ... if (wdata.has_outprofile) { // ... if overlap exists here, needing interp flag_xhi = (interp_to_mfab::check_lateral_overlap_hi( - wdata.beach_length, dir, m_ow_velocity.vec_ptrs(), + wdata.beach_length, dir, ow_velocity.vec_ptrs(), geom) == 1); } + // Wave generation zones are irrelevant for single-phase mode, + // indicated by the use for a forcing term + if (for_forcing_term) { + flag_xlo = true; + flag_xhi = true; + } + if (flag_z && (flag_xlo || flag_xhi)) { // Interpolation is needed wdata.do_interp = true; @@ -620,30 +631,30 @@ struct UpdateRelaxZonesOp if (double_data == 1) { if (wdata.do_interp) { populate_fields_all_levels( - wdata, geom, m_ow_levelset, m_ow_velocity, -1); + wdata, geom, ow_levelset, ow_velocity, -1); } // Average down to get fine information on coarse grid where // possible (may be unnecessary) for (int lev = nlevels - 1; lev > 0; --lev) { amrex::average_down( - m_ow_velocity(lev), m_ow_velocity(lev - 1), 0, + ow_velocity(lev), ow_velocity(lev - 1), 0, AMREX_SPACEDIM, sim.mesh().refRatio(lev - 1)); amrex::average_down( - m_ow_levelset(lev), m_ow_levelset(lev - 1), 0, 1, + ow_levelset(lev), ow_levelset(lev - 1), 0, 1, sim.mesh().refRatio(lev - 1)); } // Fill patch to get correct ghost cells after average down - m_ow_velocity.fillpatch(sim.time().new_time()); - m_ow_levelset.fillpatch(sim.time().new_time()); + ow_velocity.fillpatch(sim.time().new_time()); + ow_levelset.fillpatch(sim.time().new_time()); // Prior t_last (corresponding to ow fields) t_last = (wdata.ntime - 1) * wdata.dt_modes; } else if (double_data == 2) { // Restarting simulation or taking a big step, new time at ntime // Initialize ow fields to 0 for time interp, will be replaced - m_ow_levelset.setVal(0.0); - m_ow_velocity.setVal(0.0); + ow_levelset.setVal(0.0); + ow_velocity.setVal(0.0); // No modification needed for t_last, leads to interp factor = 1 } @@ -670,28 +681,28 @@ struct UpdateRelaxZonesOp // Temporally interpolate at every timestep to get target solution for (int lev = 0; lev < nlevels; ++lev) { - auto phi = m_ow_levelset(lev).arrays(); - auto vel = m_ow_velocity(lev).arrays(); - auto W2A_phi = w2a_levelset(lev).arrays(); - auto W2A_vel = w2a_velocity(lev).arrays(); + auto phi = ow_levelset(lev).arrays(); + auto vel = ow_velocity(lev).arrays(); + const auto W2A_phi = w2a_levelset(lev).const_arrays(); + const auto W2A_vel = w2a_velocity(lev).const_arrays(); const amrex::Real W2A_t = wdata.t; amrex::ParallelFor( - m_ow_levelset(lev), amrex::IntVect(3), + ow_levelset(lev), amrex::IntVect(3), [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { // Interpolate temporally every time phi[nbx](i, j, k) += (W2A_phi[nbx](i, j, k) - phi[nbx](i, j, k)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); + (time - t_last) / (W2A_t - t_last + 1e-16); vel[nbx](i, j, k, 0) += (W2A_vel[nbx](i, j, k, 0) - vel[nbx](i, j, k, 0)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); + (time - t_last) / (W2A_t - t_last + 1e-16); vel[nbx](i, j, k, 1) += (W2A_vel[nbx](i, j, k, 1) - vel[nbx](i, j, k, 1)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); + (time - t_last) / (W2A_t - t_last + 1e-16); vel[nbx](i, j, k, 2) += (W2A_vel[nbx](i, j, k, 2) - vel[nbx](i, j, k, 2)) * - (newtime - t_last) / (W2A_t - t_last + 1e-16); + (time - t_last) / (W2A_t - t_last + 1e-16); }); } #else From 32783ac3ab55d63e923235bc9904e89fbfbac824 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 12:01:22 -0700 Subject: [PATCH 12/31] renaming update_relax_zones to be more accurate --- amr-wind/ocean_waves/OceanWaves.cpp | 4 ++-- amr-wind/ocean_waves/OceanWavesModel.H | 6 +++--- amr-wind/ocean_waves/OceanWavesOps.H | 2 +- amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H | 2 +- .../ocean_waves/relaxation_zones/relaxation_zones_ops.H | 2 +- amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H | 2 +- amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H | 2 +- 7 files changed, 10 insertions(+), 10 deletions(-) diff --git a/amr-wind/ocean_waves/OceanWaves.cpp b/amr-wind/ocean_waves/OceanWaves.cpp index 45e176c9e8..45268531e1 100644 --- a/amr-wind/ocean_waves/OceanWaves.cpp +++ b/amr-wind/ocean_waves/OceanWaves.cpp @@ -81,7 +81,7 @@ void OceanWaves::pre_advance_work() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_advance_work"); if (!m_multiphase_mode) { - m_owm->update_relax_zones(true); + m_owm->update_target_fields(true); } } @@ -99,7 +99,7 @@ void OceanWaves::post_advance_work() void OceanWaves::relaxation_zones() { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_relaxation_zones"); - m_owm->update_relax_zones(false); + m_owm->update_target_fields(false); m_owm->apply_relax_zones(); m_owm->reset_regrid_flag(); } diff --git a/amr-wind/ocean_waves/OceanWavesModel.H b/amr-wind/ocean_waves/OceanWavesModel.H index 45e5c4dead..c55be252bc 100644 --- a/amr-wind/ocean_waves/OceanWavesModel.H +++ b/amr-wind/ocean_waves/OceanWavesModel.H @@ -34,7 +34,7 @@ public: virtual void init_waves(int, const amrex::Geometry&, bool) = 0; - virtual void update_relax_zones(bool) = 0; + virtual void update_target_fields(bool) = 0; virtual void apply_relax_zones() = 0; @@ -91,9 +91,9 @@ public: m_out_op.read_io_options(pp); } - void update_relax_zones(bool for_forcing_term) override + void update_target_fields(bool for_forcing_term) override { - ops::UpdateRelaxZonesOp()(m_data, for_forcing_term); + ops::UpdateTargetFieldsOp()(m_data, for_forcing_term); } void apply_relax_zones() override diff --git a/amr-wind/ocean_waves/OceanWavesOps.H b/amr-wind/ocean_waves/OceanWavesOps.H index bb48cc8d3f..0c1a1760c5 100644 --- a/amr-wind/ocean_waves/OceanWavesOps.H +++ b/amr-wind/ocean_waves/OceanWavesOps.H @@ -32,7 +32,7 @@ template struct InitDataOp; template -struct UpdateRelaxZonesOp; +struct UpdateTargetFieldsOp; template struct ApplyRelaxZonesOp; diff --git a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H index 41a01d8f43..439e41dd80 100644 --- a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H @@ -120,7 +120,7 @@ struct InitDataOp }; // namespace amr_wind::ocean_waves::ops template <> -struct UpdateRelaxZonesOp +struct UpdateTargetFieldsOp { void operator()(LinearWaves::DataType& data, bool for_forcing_term) { diff --git a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H index e015c61057..a97e200031 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H @@ -51,7 +51,7 @@ struct UseDefaultOp }; template -struct UpdateRelaxZonesOp< +struct UpdateTargetFieldsOp< WaveTheoryTrait, typename std::enable_if_t< std::is_base_of_v>> diff --git a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H index 2a5a42ecd0..0d6e2b4988 100644 --- a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H @@ -145,7 +145,7 @@ struct InitDataOp }; template <> -struct UpdateRelaxZonesOp +struct UpdateTargetFieldsOp { void operator()(StokesWaves::DataType& data, bool for_forcing_term) { diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index dedbf36943..1988fd7497 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -501,7 +501,7 @@ struct InitDataOp }; // namespace ocean_waves template <> -struct UpdateRelaxZonesOp +struct UpdateTargetFieldsOp { // cppcheck-suppress constParameterReference void operator()(W2AWaves::DataType& data, bool for_forcing_term) From fa60d17df4b3f399b0729a4c7077c2ed036f8596 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 14:17:35 -0700 Subject: [PATCH 13/31] init tweaks to w2a --- .../relaxation_zones/waves2amr_ops.H | 28 +++++++++++++------ 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index 1988fd7497..990c476cc4 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -406,7 +406,8 @@ struct InitDataOp void // cppcheck-suppress constParameterReference operator()( - W2AWaves::DataType & data, int level, const amrex::Geometry & geom, bool multiphase_mode) + W2AWaves::DataType & + data, int level, const amrex::Geometry & geom, bool multiphase_mode) { #ifdef AMR_WIND_USE_W2A @@ -417,9 +418,13 @@ struct InitDataOp // Fill ow fields, then populate flow fields according to setup auto& ow_levelset = sim.repo().get_field("ow_levelset"); auto& ow_velocity = sim.repo().get_field("ow_velocity"); - auto& levelset = sim.repo().get_field("levelset"); auto& velocity = sim.repo().get_field("velocity"); + Field* levelset{nullptr}; + if (multiphase_mode) { + levelset = &sim.repo().get_field("levelset"); + } + const auto& problo = geom.ProbLoArray(); const auto& probhi = geom.ProbHiArray(); const auto& dx = geom.CellSizeArray(); @@ -448,15 +453,17 @@ struct InitDataOp // Populate flow fields according to intended forcing and init setup const auto& ow_phi = ow_levelset(level).const_arrays(); const auto& ow_vel = ow_velocity(level).const_arrays(); - const auto& phi = levelset(level).arrays(); const auto& vel = velocity(level).arrays(); + const auto& phi_arrs = multiphase_mode + ? (*levelset)(level).arrays() + : amrex::MultiArray4(); const amrex::Real gen_length = wdata.gen_length; const amrex::Real beach_length = wdata.beach_length; const amrex::Real zero_sea_level = wdata.zsl; - const bool has_beach = wdata.has_beach; - const bool init_wave_field = wdata.init_wave_field; + const bool has_beach = wdata.has_beach && multiphase_mode; + const bool init_wave_field = wdata.init_wave_field || !multiphase_mode; amrex::ParallelFor( velocity(level), amrex::IntVect(3), @@ -479,14 +486,17 @@ struct InitDataOp x, problo[0], gen_length, probhi[0], beach_length, wave_sol, bulk, outlet); - phi[nbx](i, j, k) = local_profile[3] - z; + const amrex::Real phi = local_profile[3] - z; const amrex::Real cell_length_2D = std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + if (phi + cell_length_2D >= 0) { vel[nbx](i, j, k, 0) = local_profile[0]; vel[nbx](i, j, k, 1) = local_profile[1]; vel[nbx](i, j, k, 2) = local_profile[2]; } + if (multiphase_mode) { + phi_arrs[nbx](i, j, k) = phi; + } }); // Start w2a fields at 0, some areas will not be modified @@ -580,8 +590,8 @@ struct UpdateTargetFieldsOp const int dir = 0; flag_xlo = (interp_to_mfab::check_lateral_overlap_lo( - wdata.gen_length, dir, ow_velocity.vec_ptrs(), - geom) == 1); + wdata.gen_length, dir, ow_velocity.vec_ptrs(), geom) == + 1); // No overlap with gen region means no interp, unless ... if (wdata.has_outprofile) { // ... if overlap exists here, needing interp From 8b4090afb505fccef4f8670cb085634bb66a77dd Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 14:43:29 -0700 Subject: [PATCH 14/31] introduce target vel to drag forcing --- .../icns/source_terms/DragForcing.H | 2 ++ .../icns/source_terms/DragForcing.cpp | 29 +++++++++++++++++-- amr-wind/physics/TerrainDrag.H | 2 +- 3 files changed, 29 insertions(+), 4 deletions(-) diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.H b/amr-wind/equation_systems/icns/source_terms/DragForcing.H index 283ccd7bdc..2561640a86 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.H +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.H @@ -34,6 +34,7 @@ private: const CFDSim& m_sim; const amrex::AmrCore& m_mesh; const Field& m_velocity; + const Field* m_target_vel{nullptr}; amrex::Gpu::DeviceVector m_device_vel_ht; amrex::Gpu::DeviceVector m_device_vel_vals; amrex::Real m_drag_coefficient{10.0}; @@ -48,6 +49,7 @@ private: int m_sponge_south{0}; int m_sponge_north{1}; bool m_is_laminar{false}; + bool m_terrain_is_waves{false}; }; } // namespace amr_wind::pde::icns diff --git a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp index 2a6c6654c1..0d81284482 100644 --- a/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp +++ b/amr-wind/equation_systems/icns/source_terms/DragForcing.cpp @@ -4,6 +4,7 @@ #include "AMReX_Gpu.H" #include "AMReX_Random.H" #include "amr-wind/wind_energy/ABL.H" +#include "amr-wind/physics/TerrainDrag.H" #include "amr-wind/utilities/linear_interpolation.H" namespace amr_wind::pde::icns { @@ -42,6 +43,13 @@ DragForcing::DragForcing(const CFDSim& sim) } else { m_sponge_strength = 0.0; } + if (phy_mgr.contains("OceanWaves")) { + const auto terrain_phys = + m_sim.physics_manager().get(); + const auto target_vel_name = terrain_phys.wave_velocity_field_name(); + m_target_vel = &sim.repo().get_field(target_vel_name); + m_terrain_is_waves = true; + } } DragForcing::~DragForcing() = default; @@ -68,6 +76,12 @@ void DragForcing::operator()( const auto& drag = (*m_terrain_drag)(lev).const_array(mfi); auto* const m_terrainz0 = &this->m_sim.repo().get_field("terrainz0"); const auto& terrainz0 = (*m_terrainz0)(lev).const_array(mfi); + + const bool is_waves = m_terrain_is_waves; + const auto& target_vel_arr = is_waves + ? (*m_target_vel)(lev).const_array(mfi) + : amrex::Array4(); + const auto& geom = m_mesh.Geom(lev); const auto& dx = geom.CellSizeArray(); const auto& prob_lo = geom.ProbLoArray(); @@ -158,20 +172,29 @@ void DragForcing::operator()( Dyz = -ustar * ustar * uy1 / (tiny + std::sqrt(ux1 * ux1 + uy1 * uy1)) / dx[2]; } + amrex::Real target_u = 0.; + amrex::Real target_v = 0.; + amrex::Real target_w = 0.; + if (is_waves) { + target_u = target_vel_arr(i, j, k, 0); + target_v = target_vel_arr(i, j, k, 1); + target_w = target_vel_arr(i, j, k, 2); + } + const amrex::Real CdM = std::min(Cd / (m + tiny), cd_max / scale_factor); src_term(i, j, k, 0) -= - (CdM * m * ux1 * blank(i, j, k) + Dxz * drag(i, j, k) + + (CdM * m * (ux1 - target_u) * blank(i, j, k) + Dxz * drag(i, j, k) + bc_forcing_x * drag(i, j, k) + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (ux1 - sponge_density * spongeVelX)); src_term(i, j, k, 1) -= - (CdM * m * uy1 * blank(i, j, k) + Dyz * drag(i, j, k) + + (CdM * m * (uy1 - target_v) * blank(i, j, k) + Dyz * drag(i, j, k) + bc_forcing_y * drag(i, j, k) + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (uy1 - sponge_density * spongeVelY)); src_term(i, j, k, 2) -= - (CdM * m * uz1 * blank(i, j, k) + + (CdM * m * (uz1 - target_w) * blank(i, j, k) + (xstart_damping + xend_damping + ystart_damping + yend_damping) * (uz1 - sponge_density * spongeVelZ)); }); diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index 344334a8fc..68289f2957 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -34,7 +34,7 @@ public: void post_advance_work() override {} - std::string wave_velocity_field_name() { return m_wave_velocity_name; } + std::string wave_velocity_field_name() const { return m_wave_velocity_name; } private: CFDSim& m_sim; From 64fed5e2bdc30a9c1559a0d77d9c48a993a02555 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 17:22:13 -0700 Subject: [PATCH 15/31] wave terrain showing up in forcing at beginning and end of reg test --- amr-wind/ocean_waves/OceanWaves.H | 2 + amr-wind/ocean_waves/OceanWaves.cpp | 25 ++- amr-wind/ocean_waves/OceanWavesModel.H | 7 + amr-wind/ocean_waves/OceanWavesOps.H | 3 + .../relaxation_zones/linear_waves_ops.H | 8 +- .../relaxation_zones/relaxation_zones_ops.H | 19 ++ .../relaxation_zones/relaxation_zones_ops.cpp | 21 ++ .../relaxation_zones/stokes_waves_ops.H | 9 +- amr-wind/physics/TerrainDrag.H | 7 +- amr-wind/physics/TerrainDrag.cpp | 208 +++++++++++------- test/CMakeLists.txt | 1 + .../abl_waves_terrain/abl_waves_terrain.inp | 75 +++++++ .../abl_waves_terrain/static_box.txt | 5 + 13 files changed, 296 insertions(+), 94 deletions(-) create mode 100644 test/test_files/abl_waves_terrain/abl_waves_terrain.inp create mode 100644 test/test_files/abl_waves_terrain/static_box.txt diff --git a/amr-wind/ocean_waves/OceanWaves.H b/amr-wind/ocean_waves/OceanWaves.H index de32157347..35669a75b0 100644 --- a/amr-wind/ocean_waves/OceanWaves.H +++ b/amr-wind/ocean_waves/OceanWaves.H @@ -52,6 +52,8 @@ protected: private: void relaxation_zones(); + void update_target_volume_fraction(); + CFDSim& m_sim; //! Unique pointer to the ocean waves model diff --git a/amr-wind/ocean_waves/OceanWaves.cpp b/amr-wind/ocean_waves/OceanWaves.cpp index 45268531e1..5d6ddf12d4 100644 --- a/amr-wind/ocean_waves/OceanWaves.cpp +++ b/amr-wind/ocean_waves/OceanWaves.cpp @@ -17,12 +17,6 @@ OceanWaves::OceanWaves(CFDSim& sim) , m_ow_velocity( sim.repo().declare_field("ow_velocity", AMREX_SPACEDIM, 3, 1)) { - if (!(sim.physics_manager().contains("MultiPhase") || - sim.physics_manager().contains("TerrainDrag"))) { - amrex::Abort( - "OceanWaves requires MultiPhase or TerrainDrag physics to be " - "active"); - } if (!sim.physics_manager().contains("MultiPhase")) { m_multiphase_mode = false; } @@ -38,6 +32,13 @@ void OceanWaves::pre_init_actions() BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_init_actions"); amrex::ParmParse pp(identifier()); + if (!(m_multiphase_mode || + m_sim.physics_manager().contains("TerrainDrag"))) { + amrex::Abort( + "OceanWaves requires MultiPhase or TerrainDrag physics to be " + "active"); + } + std::string label; pp.query("label", label); const std::string& tname = label; @@ -68,6 +69,9 @@ void OceanWaves::post_init_actions() BL_PROFILE("amr-wind::ocean_waves::OceanWaves::post_init_actions"); if (m_multiphase_mode) { relaxation_zones(); + } else { + m_owm->update_target_fields(true); + m_owm->update_target_volume_fraction(); } } @@ -82,6 +86,7 @@ void OceanWaves::pre_advance_work() BL_PROFILE("amr-wind::ocean_waves::OceanWaves::pre_advance_work"); if (!m_multiphase_mode) { m_owm->update_target_fields(true); + m_owm->update_target_volume_fraction(); } } @@ -98,12 +103,18 @@ void OceanWaves::post_advance_work() */ void OceanWaves::relaxation_zones() { - BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_relaxation_zones"); + BL_PROFILE("amr-wind::ocean_waves::OceanWaves::relaxation_zones"); m_owm->update_target_fields(false); m_owm->apply_relax_zones(); m_owm->reset_regrid_flag(); } +void OceanWaves::update_target_volume_fraction() +{ + BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_target_volume_fraction"); + m_owm->update_target_volume_fraction(); +} + void OceanWaves::prepare_outputs() { const std::string post_dir = m_sim.io_manager().post_processing_directory(); diff --git a/amr-wind/ocean_waves/OceanWavesModel.H b/amr-wind/ocean_waves/OceanWavesModel.H index c55be252bc..077478c578 100644 --- a/amr-wind/ocean_waves/OceanWavesModel.H +++ b/amr-wind/ocean_waves/OceanWavesModel.H @@ -38,6 +38,8 @@ public: virtual void apply_relax_zones() = 0; + virtual void update_target_volume_fraction() = 0; + virtual void prepare_outputs(const std::string&) = 0; virtual void write_outputs() = 0; @@ -101,6 +103,11 @@ public: ops::ApplyRelaxZonesOp()(m_data); } + void update_target_volume_fraction() override + { + ops::UpdateTargetVolumeFractionOp()(m_data); + } + void prepare_outputs(const std::string& out_dir) override { m_out_op.prepare_outputs(out_dir); diff --git a/amr-wind/ocean_waves/OceanWavesOps.H b/amr-wind/ocean_waves/OceanWavesOps.H index 0c1a1760c5..529d508465 100644 --- a/amr-wind/ocean_waves/OceanWavesOps.H +++ b/amr-wind/ocean_waves/OceanWavesOps.H @@ -37,6 +37,9 @@ struct UpdateTargetFieldsOp; template struct ApplyRelaxZonesOp; +template +struct UpdateTargetVolumeFractionOp; + template struct ProcessOutputsOp; diff --git a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H index 439e41dd80..d97306db33 100644 --- a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H @@ -1,7 +1,6 @@ #ifndef LINEAR_WAVES_OPS_H #define LINEAR_WAVES_OPS_H -#include "amr-wind/physics/multiphase/MultiPhase.H" #include "amr-wind/ocean_waves/relaxation_zones/LinearWaves.H" #include "amr-wind/ocean_waves/OceanWavesOps.H" #include "amr-wind/ocean_waves/OceanWaves.H" @@ -19,8 +18,11 @@ struct ReadInputsOp auto& wdata = data.meta(); auto& info = data.info(); relaxation_zones::read_inputs(wdata, info, pp); - // Get gravity from MultiPhase physics, assume negative z - wdata.g = -data.sim().physics_manager().get().gravity()[2]; + // Get gravity, assume negative z + amrex::Vector gravity{0.0, 0.0, -9.81}; + amrex::ParmParse pp_incflo("incflo"); + pp_incflo.queryarr("gravity", gravity); + wdata.g = -gravity[2]; pp.get("wave_length", wdata.wave_length); pp.get("wave_height", wdata.wave_height); diff --git a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H index a97e200031..1cb08c0efd 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H @@ -28,6 +28,8 @@ void init_data_structures(CFDSim&); */ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata); +void update_target_vof(CFDSim& sim, const RelaxZonesBaseData& wdata); + void prepare_netcdf_file( const std::string& /*ncfile*/, const RelaxZonesBaseData& /*meta*/, @@ -76,6 +78,23 @@ struct ApplyRelaxZonesOp< } }; +template +struct UpdateTargetVolumeFractionOp< + WaveTheoryTrait, + typename std::enable_if_t< + std::is_base_of_v>> +{ + void operator()(typename WaveTheoryTrait::DataType& data) + { + BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_target_vof"); + + const auto& wdata = data.meta(); + auto& sim = data.sim(); + + relaxation_zones::update_target_vof(sim, wdata); + } +}; + template struct ProcessOutputsOp< WaveTheoryTrait, 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 4c98d83a0b..e049736cb3 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp @@ -228,6 +228,27 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata) density.fillpatch(time); } +void update_target_vof(CFDSim& sim, const RelaxZonesBaseData& wdata) +{ + const int nlevels = sim.repo().num_active_levels(); + auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); + auto& m_ow_vof = sim.repo().get_field("ow_vof"); + const auto& geom = sim.mesh().Geom(); + + for (int lev = 0; lev < nlevels; ++lev) { + const auto& dx = geom[lev].CellSizeArray(); + const auto target_phi = m_ow_levelset(lev).arrays(); + auto target_volfrac = m_ow_vof(lev).arrays(); + const amrex::Real eps = 2. * std::cbrt(dx[0] * dx[1] * dx[2]); + amrex::ParallelFor( + m_ow_vof(lev), amrex::IntVect(2), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + target_volfrac[nbx](i, j, k) = + multiphase::levelset_to_vof(i, j, k, eps, target_phi[nbx]); + }); + } +} + void prepare_netcdf_file( const std::string& ncfile, const RelaxZonesBaseData& meta, diff --git a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H index 0d6e2b4988..1281c3b237 100644 --- a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H @@ -1,7 +1,6 @@ #ifndef STOKES_WAVES_OPS_H #define STOKES_WAVES_OPS_H -#include "amr-wind/physics/multiphase/MultiPhase.H" #include "amr-wind/ocean_waves/relaxation_zones/StokesWaves.H" #include "amr-wind/ocean_waves/relaxation_zones/stokes_waves_K.H" #include "amr-wind/ocean_waves/OceanWavesOps.H" @@ -21,8 +20,12 @@ struct ReadInputsOp auto& info = data.info(); relaxation_zones::read_inputs(wdata, info, pp); - // Get gravity from MultiPhase physics, assume negative z - wdata.g = -data.sim().physics_manager().get().gravity()[2]; + // Get gravity, assume negative z + amrex::Vector gravity{0.0, 0.0, -9.81}; + amrex::ParmParse pp_incflo("incflo"); + pp_incflo.queryarr("gravity", gravity); + wdata.g = -gravity[2]; + // Get wave attributes pp.get("wave_height", wdata.wave_height); pp.get("order", wdata.order); diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index 68289f2957..e33cb8bc80 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -26,11 +26,11 @@ public: void pre_init_actions() override {} - void post_init_actions() override {} + void post_init_actions() override; void post_regrid_actions() override; - void pre_advance_work() override {} + void pre_advance_work() override; void post_advance_work() override {} @@ -57,13 +57,14 @@ private: amrex::Vector m_yrough; amrex::Vector m_z0rough; - //! Terrain Drag Force Term + //! Terrain Drag for waves bool m_terrain_is_waves{false}; Field* m_wave_volume_fraction{nullptr}; Field* m_wave_negative_elevation{nullptr}; const std::string m_wave_volume_fraction_name{"ow_vof"}; const std::string m_wave_negative_elevation_name{"ow_levelset"}; const std::string m_wave_velocity_name{"ow_velocity"}; + void convert_waves_to_blank_and_drag_flags(); }; } // namespace amr_wind::terraindrag diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 78ff7a6481..04c6e8afd5 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -61,6 +61,10 @@ TerrainDrag::TerrainDrag(CFDSim& sim) void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) { + if (m_terrain_is_waves) { + return; + } + BL_PROFILE("amr-wind::" + this->identifier() + "::initialize_fields"); const auto& dx = geom.CellSizeArray(); const auto& prob_lo = geom.ProbLoArray(); @@ -74,68 +78,125 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) auto levelz0 = terrainz0.arrays(); auto levelheight = terrain_height.arrays(); + const auto xterrain_size = m_xterrain.size(); + const auto yterrain_size = m_yterrain.size(); + const auto zterrain_size = m_zterrain.size(); + amrex::Gpu::DeviceVector device_xterrain(xterrain_size); + amrex::Gpu::DeviceVector device_yterrain(yterrain_size); + amrex::Gpu::DeviceVector device_zterrain(zterrain_size); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), + device_xterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), + device_yterrain.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), + device_zterrain.begin()); + const auto* xterrain_ptr = device_xterrain.data(); + const auto* yterrain_ptr = device_yterrain.data(); + const auto* zterrain_ptr = device_zterrain.data(); + // Copy Roughness to gpu + const auto xrough_size = m_xrough.size(); + const auto yrough_size = m_yrough.size(); + const auto z0rough_size = m_z0rough.size(); + amrex::Gpu::DeviceVector device_xrough(xrough_size); + amrex::Gpu::DeviceVector device_yrough(yrough_size); + amrex::Gpu::DeviceVector device_z0rough(z0rough_size); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_xrough.begin(), m_xrough.end(), + device_xrough.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_yrough.begin(), m_yrough.end(), + device_yrough.begin()); + amrex::Gpu::copy( + amrex::Gpu::hostToDevice, m_z0rough.begin(), m_z0rough.end(), + device_z0rough.begin()); + const auto* xrough_ptr = device_xrough.data(); + const auto* yrough_ptr = device_yrough.data(); + const auto* z0rough_ptr = device_z0rough.data(); + + amrex::ParallelFor( + blanking, m_terrain_blank.num_grow(), + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; + const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; + const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; + const amrex::Real terrainHt = interp::bilinear( + xterrain_ptr, xterrain_ptr + xterrain_size, yterrain_ptr, + yterrain_ptr + yterrain_size, zterrain_ptr, x, y); + levelBlanking[nbx](i, j, k, 0) = + static_cast((z <= terrainHt) && (z > prob_lo[2])); + levelheight[nbx](i, j, k, 0) = + std::max(std::abs(z - terrainHt), 0.5 * dx[2]); + + amrex::Real roughz0 = 0.1; + if (xrough_size > 0) { + roughz0 = interp::bilinear( + xrough_ptr, xrough_ptr + xrough_size, yrough_ptr, + yrough_ptr + yrough_size, z0rough_ptr, x, y); + } + levelz0[nbx](i, j, k, 0) = roughz0; + }); + + amrex::ParallelFor( + blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + if ((levelBlanking[nbx](i, j, k, 0) == 0) && (k > 0) && + (levelBlanking[nbx](i, j, k - 1, 0) == 1)) { + levelDrag[nbx](i, j, k, 0) = 1; + } else { + levelDrag[nbx](i, j, k, 0) = 0; + } + }); +} + +void TerrainDrag::post_init_actions() +{ if (!m_terrain_is_waves) { - const auto xterrain_size = m_xterrain.size(); - const auto yterrain_size = m_yterrain.size(); - const auto zterrain_size = m_zterrain.size(); - amrex::Gpu::DeviceVector device_xterrain(xterrain_size); - amrex::Gpu::DeviceVector device_yterrain(yterrain_size); - amrex::Gpu::DeviceVector device_zterrain(zterrain_size); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_xterrain.begin(), m_xterrain.end(), - device_xterrain.begin()); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_yterrain.begin(), m_yterrain.end(), - device_yterrain.begin()); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_zterrain.begin(), m_zterrain.end(), - device_zterrain.begin()); - const auto* xterrain_ptr = device_xterrain.data(); - const auto* yterrain_ptr = device_yterrain.data(); - const auto* zterrain_ptr = device_zterrain.data(); - // Copy Roughness to gpu - const auto xrough_size = m_xrough.size(); - const auto yrough_size = m_yrough.size(); - const auto z0rough_size = m_z0rough.size(); - amrex::Gpu::DeviceVector device_xrough(xrough_size); - amrex::Gpu::DeviceVector device_yrough(yrough_size); - amrex::Gpu::DeviceVector device_z0rough(z0rough_size); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_xrough.begin(), m_xrough.end(), - device_xrough.begin()); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_yrough.begin(), m_yrough.end(), - device_yrough.begin()); - amrex::Gpu::copy( - amrex::Gpu::hostToDevice, m_z0rough.begin(), m_z0rough.end(), - device_z0rough.begin()); - const auto* xrough_ptr = device_xrough.data(); - const auto* yrough_ptr = device_yrough.data(); - const auto* z0rough_ptr = device_z0rough.data(); + return; + } + BL_PROFILE("amr-wind::" + this->identifier() + "::post_init_actions"); + convert_waves_to_blank_and_drag_flags(); +} - amrex::ParallelFor( - blanking, m_terrain_blank.num_grow(), - [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0]; - const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1]; - const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; - const amrex::Real terrainHt = interp::bilinear( - xterrain_ptr, xterrain_ptr + xterrain_size, yterrain_ptr, - yterrain_ptr + yterrain_size, zterrain_ptr, x, y); - levelBlanking[nbx](i, j, k, 0) = - static_cast((z <= terrainHt) && (z > prob_lo[2])); - levelheight[nbx](i, j, k, 0) = - std::max(std::abs(z - terrainHt), 0.5 * dx[2]); - - amrex::Real roughz0 = 0.1; - if (xrough_size > 0) { - roughz0 = interp::bilinear( - xrough_ptr, xrough_ptr + xrough_size, yrough_ptr, - yrough_ptr + yrough_size, z0rough_ptr, x, y); - } - levelz0[nbx](i, j, k, 0) = roughz0; - }); +void TerrainDrag::pre_advance_work() +{ + if (!m_terrain_is_waves) { + return; + } + BL_PROFILE("amr-wind::" + this->identifier() + "::pre_advance_work"); + convert_waves_to_blank_and_drag_flags(); +} + +void TerrainDrag::post_regrid_actions() +{ + if (m_terrain_is_waves) { + convert_waves_to_blank_and_drag_flags(); } else { + const int nlevels = m_sim.repo().num_active_levels(); + for (int lev = 0; lev < nlevels; ++lev) { + initialize_fields(lev, m_sim.repo().mesh().Geom(lev)); + } + } +} + +void TerrainDrag::convert_waves_to_blank_and_drag_flags() +{ + const int nlevels = m_sim.repo().num_active_levels(); + for (int level = 0; level < nlevels; ++level) { + const auto geom = m_sim.repo().mesh().Geom(level); + const auto& dx = geom.CellSizeArray(); + const auto& prob_lo = geom.ProbLoArray(); + auto& blanking = m_terrain_blank(level); + auto& terrainz0 = m_terrainz0(level); + auto& terrain_height = m_terrain_height(level); + auto& drag = m_terrain_drag(level); + + auto levelBlanking = blanking.arrays(); + auto levelDrag = drag.arrays(); + auto levelz0 = terrainz0.arrays(); + auto levelheight = terrain_height.arrays(); + const auto negative_wave_elevation = (*m_wave_negative_elevation)(level).const_arrays(); const auto wave_vol_frac = @@ -151,26 +212,17 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) levelheight[nbx](i, j, k, 0) = std::max( std::abs(negative_wave_elevation[nbx](i, j, k)), 0.5 * dx[2]); - - // Set uniform roughness (ocean surface) - levelz0[nbx](i, j, k, 0) = 1e-4; }); - } - amrex::ParallelFor( - blanking, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { - if ((levelBlanking[nbx](i, j, k, 0) == 0) && (k > 0) && - (levelBlanking[nbx](i, j, k - 1, 0) == 1)) { - levelDrag[nbx](i, j, k, 0) = 1; - } else { - levelDrag[nbx](i, j, k, 0) = 0; - } - }); -} -void TerrainDrag::post_regrid_actions() -{ - const int nlevels = m_sim.repo().num_active_levels(); - for (int lev = 0; lev < nlevels; ++lev) { - initialize_fields(lev, m_sim.repo().mesh().Geom(lev)); + amrex::ParallelFor( + blanking, + [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept { + if ((levelBlanking[nbx](i, j, k, 0) == 0) && (k > 0) && + (levelBlanking[nbx](i, j, k - 1, 0) == 1)) { + levelDrag[nbx](i, j, k, 0) = 1; + } else { + levelDrag[nbx](i, j, k, 0) = 0; + } + }); } } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8799e0eeee..1dd638918c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -206,6 +206,7 @@ add_test_re(abl_unstable_local_wall_model) add_test_re(abl_unstable_schumann_wall_model) add_test_re(abl_anelastic) add_test_re(abl_multiphase_laminar) +add_test_re(abl_waves_terrain) add_test_re(act_abl_joukowskydisk) add_test_re(act_abl_uniformctdisk) add_test_re(act_fixed_wing) diff --git a/test/test_files/abl_waves_terrain/abl_waves_terrain.inp b/test/test_files/abl_waves_terrain/abl_waves_terrain.inp new file mode 100644 index 0000000000..61633d9414 --- /dev/null +++ b/test/test_files/abl_waves_terrain/abl_waves_terrain.inp @@ -0,0 +1,75 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# 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 # +#.......................................# +time.plot_interval = 10 # Steps between plot files +time.checkpoint_interval = -1 # Steps between checkpoint files +io.outputs = velocity_mueff temperature_mueff ow_vof ow_levelset ow_velocity +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0.0 0.0 -9.81 # Gravitational force (3D) +incflo.density = 1.225 +transport.model = ConstTransport +transport.viscosity = 1e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.333 +turbulence.model = Kosovic +Kosovic.refMOL = -1e30 +incflo.physics = ABL OceanWaves TerrainDrag +ICNS.source_terms = BoussinesqBuoyancy CoriolisForcing GeostrophicForcing NonLinearSGSTerm DragForcing +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 +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 + +OceanWaves.label = Wave1 +OceanWaves.Wave1.type = LinearWaves +OceanWaves.Wave1.wave_height=2.0 +OceanWaves.Wave1.wave_length=40.0 +OceanWaves.Wave1.water_depth=400 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 48 48 48 # 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. -20. # Lo corner coordinates +geometry.prob_hi = 400. 400. 380. # Hi corner coordinates +geometry.is_periodic = 1 1 0 # Periodicity x y z (0/1) +# Boundary conditions +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_waves_terrain/static_box.txt b/test/test_files/abl_waves_terrain/static_box.txt new file mode 100644 index 0000000000..fa81648cb5 --- /dev/null +++ b/test/test_files/abl_waves_terrain/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. From 983a7b7948e0a96caeafa9cefbb3f9614ad255f2 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 17:30:00 -0700 Subject: [PATCH 16/31] formatting --- amr-wind/ocean_waves/OceanWaves.cpp | 3 ++- amr-wind/ocean_waves/OceanWavesModel.H | 6 ++++-- amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H | 5 +++-- amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H | 7 ++++--- 4 files changed, 13 insertions(+), 8 deletions(-) diff --git a/amr-wind/ocean_waves/OceanWaves.cpp b/amr-wind/ocean_waves/OceanWaves.cpp index 5d6ddf12d4..364d67e126 100644 --- a/amr-wind/ocean_waves/OceanWaves.cpp +++ b/amr-wind/ocean_waves/OceanWaves.cpp @@ -111,7 +111,8 @@ void OceanWaves::relaxation_zones() void OceanWaves::update_target_volume_fraction() { - BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_target_volume_fraction"); + BL_PROFILE( + "amr-wind::ocean_waves::OceanWaves::update_target_volume_fraction"); m_owm->update_target_volume_fraction(); } diff --git a/amr-wind/ocean_waves/OceanWavesModel.H b/amr-wind/ocean_waves/OceanWavesModel.H index 077478c578..0a07fce81c 100644 --- a/amr-wind/ocean_waves/OceanWavesModel.H +++ b/amr-wind/ocean_waves/OceanWavesModel.H @@ -119,9 +119,11 @@ public: void write_outputs() override { m_out_op.write_outputs(); } - void init_waves(int level, const amrex::Geometry& geom, bool multiphase_mode) override + void init_waves( + int level, const amrex::Geometry& geom, bool multiphase_mode) override { - ops::InitDataOp()(m_data, level, geom, multiphase_mode); + ops::InitDataOp()( + m_data, level, geom, multiphase_mode); } }; diff --git a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H index d97306db33..afa16dc584 100644 --- a/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H @@ -54,8 +54,9 @@ struct InitDataOp const auto& dx = geom.CellSizeArray(); const auto& vel = velocity(level).arrays(); - const auto& phi_arrs = multiphase_mode ? (*levelset)(level).arrays() - : amrex::MultiArray4(); + const auto& phi_arrs = multiphase_mode + ? (*levelset)(level).arrays() + : amrex::MultiArray4(); const amrex::Real zero_sea_level = wdata.zsl; const amrex::Real gen_length = wdata.gen_length; diff --git a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H index 1281c3b237..706b3f5cdc 100644 --- a/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/stokes_waves_ops.H @@ -25,7 +25,7 @@ struct ReadInputsOp amrex::ParmParse pp_incflo("incflo"); pp_incflo.queryarr("gravity", gravity); wdata.g = -gravity[2]; - + // Get wave attributes pp.get("wave_height", wdata.wave_height); pp.get("order", wdata.order); @@ -94,8 +94,9 @@ struct InitDataOp const auto& dx = geom.CellSizeArray(); const auto& vel = velocity(level).arrays(); - const auto& phi_arrs = multiphase_mode ? (*levelset)(level).arrays() - : amrex::MultiArray4(); + const auto& phi_arrs = multiphase_mode + ? (*levelset)(level).arrays() + : amrex::MultiArray4(); const amrex::Real wave_height = wdata.wave_height; const amrex::Real wave_length = wdata.wave_length; From 86ffb47e41bedaffcfa7afafecaaf5b5ef547481 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 21:31:47 -0700 Subject: [PATCH 17/31] more formatting --- amr-wind/physics/TerrainDrag.H | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index e33cb8bc80..0b56516a7b 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -34,7 +34,10 @@ public: void post_advance_work() override {} - std::string wave_velocity_field_name() const { return m_wave_velocity_name; } + std::string wave_velocity_field_name() const + { + return m_wave_velocity_name; + } private: CFDSim& m_sim; From c2a0b966ab2b91b8462d8fba0cee3cf65c38cd9b Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 21:39:04 -0700 Subject: [PATCH 18/31] correct terrain height --- amr-wind/physics/TerrainDrag.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 04c6e8afd5..ca257d716b 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -209,9 +209,8 @@ void TerrainDrag::convert_waves_to_blank_and_drag_flags() const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; levelBlanking[nbx](i, j, k, 0) = static_cast( (wave_vol_frac[nbx](i, j, k) >= 0.5) && (z > prob_lo[2])); - levelheight[nbx](i, j, k, 0) = std::max( - std::abs(negative_wave_elevation[nbx](i, j, k)), - 0.5 * dx[2]); + levelheight[nbx](i, j, k, 0) = + -negative_wave_elevation[nbx](i, j, k); }); amrex::ParallelFor( blanking, From 3904a337d927513cfbfbf80bc06bf84062898de5 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Fri, 22 Nov 2024 21:40:28 -0700 Subject: [PATCH 19/31] consistent case --- amr-wind/physics/TerrainDrag.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index ca257d716b..d6201ddba1 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -76,7 +76,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) auto levelBlanking = blanking.arrays(); auto levelDrag = drag.arrays(); auto levelz0 = terrainz0.arrays(); - auto levelheight = terrain_height.arrays(); + auto levelHeight = terrain_height.arrays(); const auto xterrain_size = m_xterrain.size(); const auto yterrain_size = m_yterrain.size(); @@ -127,7 +127,7 @@ void TerrainDrag::initialize_fields(int level, const amrex::Geometry& geom) yterrain_ptr + yterrain_size, zterrain_ptr, x, y); levelBlanking[nbx](i, j, k, 0) = static_cast((z <= terrainHt) && (z > prob_lo[2])); - levelheight[nbx](i, j, k, 0) = + levelHeight[nbx](i, j, k, 0) = std::max(std::abs(z - terrainHt), 0.5 * dx[2]); amrex::Real roughz0 = 0.1; @@ -195,7 +195,7 @@ void TerrainDrag::convert_waves_to_blank_and_drag_flags() auto levelBlanking = blanking.arrays(); auto levelDrag = drag.arrays(); auto levelz0 = terrainz0.arrays(); - auto levelheight = terrain_height.arrays(); + auto levelHeight = terrain_height.arrays(); const auto negative_wave_elevation = (*m_wave_negative_elevation)(level).const_arrays(); @@ -209,7 +209,7 @@ void TerrainDrag::convert_waves_to_blank_and_drag_flags() const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2]; levelBlanking[nbx](i, j, k, 0) = static_cast( (wave_vol_frac[nbx](i, j, k) >= 0.5) && (z > prob_lo[2])); - levelheight[nbx](i, j, k, 0) = + levelHeight[nbx](i, j, k, 0) = -negative_wave_elevation[nbx](i, j, k); }); amrex::ParallelFor( From 9dc581d1353d07622470842f902cdd1e0421b57f Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 21:55:58 -0700 Subject: [PATCH 20/31] add dragtempforcing to reg test --- test/test_files/abl_waves_terrain/abl_waves_terrain.inp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_files/abl_waves_terrain/abl_waves_terrain.inp b/test/test_files/abl_waves_terrain/abl_waves_terrain.inp index 61633d9414..f6476b7aae 100644 --- a/test/test_files/abl_waves_terrain/abl_waves_terrain.inp +++ b/test/test_files/abl_waves_terrain/abl_waves_terrain.inp @@ -33,6 +33,9 @@ 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 +temperature.source_terms = DragTempForcing +DragTempForcing.reference_temperature = 265.0 + incflo.velocity = 8.0 0.0 0.0 ABL.reference_temperature = 263.5 ABL.temperature_heights = 0.0 100 400.0 From e309c39dafc6bbd337c5f5e4dc92671474160dad Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 21:57:05 -0700 Subject: [PATCH 21/31] progress on abl + w2a terrain reg test --- test/CMakeLists.txt | 1 + .../abl_w2a_terrain/abl_w2a_terrain.inp | 98 +++++++++++++++++++ .../abl_w2a_terrain/static_box.refine | 5 + 3 files changed, 104 insertions(+) create mode 100644 test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp create mode 100644 test/test_files/abl_w2a_terrain/static_box.refine diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 1dd638918c..3829504144 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -320,6 +320,7 @@ endif() if(AMR_WIND_ENABLE_W2A) add_test_re(ow_w2a) add_test_re(abl_multiphase_w2a) + add_test_red(abl_w2a_terrain abl_bndry_output_native) endif() #============================================================================= diff --git a/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp b/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp new file mode 100644 index 0000000000..1e50d4bef8 --- /dev/null +++ b/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp @@ -0,0 +1,98 @@ +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# SIMULATION STOP # +#.......................................# +time.stop_time = 12.25 # Max (simulated) time to evolve +time.max_step = 10 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# TIME STEP COMPUTATION # +#.......................................# +time.fixed_dt = 0.2 # Use this constant dt if > 0 +time.cfl = 0.95 # CFL factor +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# INPUT AND OUTPUT # +#.......................................# +io.restart_file = "../abl_bndry_output_native/chk00005" +time.plot_interval = 10 # Steps between plot files +time.checkpoint_interval = -1 # Steps between checkpoint files +io.outputs = density velocity p ow_levelset ow_velocity + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# PHYSICS # +#.......................................# +incflo.gravity = 0.0 0.0 -9.81 # Gravitational force (3D) +incflo.density = 1.225 +transport.model = ConstTransport +transport.viscosity = 1e-5 +transport.laminar_prandtl = 0.7 +transport.turbulent_prandtl = 0.333 +turbulence.model = Smagorinsky +Smagorinsky_coeffs.Cs = 0.135 + +incflo.physics = ABL OceanWaves TerrainDrag +ICNS.source_terms = CoriolisForcing DragForcing + +ABL.reference_temperature = 290.0 +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 = 125663.706143592 +GeostrophicForcing.geostrophic_wind = 10.0 0.0 0.0 +DragForcing.sponge_strength = 0 + +temperature.source_terms = DragTempForcing +DragTempForcing.reference_temperature = 265.0 + +incflo.velocity = 10.0 0.0 0.0 +ABL.temperature_heights = 0.0 2000.0 +ABL.temperature_values = 290.0 290.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 = 1.0 +ABL.deltaV = 1.0 +ABL.kappa = .41 +ABL.surface_roughness_z0 = 0.01 +ABL.bndry_file = "../abl_bndry_output_native/bndry_files" +ABL.bndry_io_mode = 1 +ABL.bndry_var_names = velocity temperature +ABL.bndry_output_format = native + +OceanWaves.label = W2A1 +OceanWaves.W2A1.type = W2AWaves +OceanWaves.W2A1.HOS_modes_filename = ../ow_w2a/modes_HOS_SWENSE.dat +# Offset to include wave troughs in domain +OceanWaves.W2A1.zero_sea_level = 0 +# These variables should change with resolution in z +OceanWaves.W2A1.number_interp_points_in_z = 3 +OceanWaves.W2A1.interp_spacing_at_surface = 2. +OceanWaves.W2A1.number_interp_above_surface = 3 + +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# ADAPTIVE MESH REFINEMENT # +#.......................................# +amr.n_cell = 48 48 48 +amr.max_level = 2 +tagging.labels = sr +tagging.sr.type = CartBoxRefinement +tagging.sr.static_refinement_def = static_box.refine +#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# +# GEOMETRY # +#.......................................# +geometry.prob_lo = 0. 0. 0. # Lo corner coordinates +geometry.prob_hi = 1000. 1000. 1000. # Hi corner coordinates +geometry.is_periodic = 0 0 0 # Periodicity x y z (0/1) +# Boundary conditions +xlo.type = "mass_inflow" +xlo.density = 1.0 +xlo.temperature = 0.0 +xhi.type = "pressure_outflow" +ylo.type = "mass_inflow" +ylo.density = 1.0 +ylo.temperature = 0.0 +yhi.type = "pressure_outflow" +zlo.type = "slip_wall" +zhi.type = "slip_wall" diff --git a/test/test_files/abl_w2a_terrain/static_box.refine b/test/test_files/abl_w2a_terrain/static_box.refine new file mode 100644 index 0000000000..a08fe9a30a --- /dev/null +++ b/test/test_files/abl_w2a_terrain/static_box.refine @@ -0,0 +1,5 @@ +2 +1 +150 150 0 1000 1000 25 +1 +170 170 0 1000 1000 20 \ No newline at end of file From ffd6fcbd1d9e1f972b6174a54d6d596aa5e57bf9 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 22:16:39 -0700 Subject: [PATCH 22/31] move function to public --- amr-wind/physics/TerrainDrag.H | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index 0b56516a7b..d5349dcfc5 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -34,6 +34,10 @@ public: void post_advance_work() override {} + //! Terrain Drag for waves + + void convert_waves_to_blank_and_drag_flags(); + std::string wave_velocity_field_name() const { return m_wave_velocity_name; @@ -67,7 +71,6 @@ private: const std::string m_wave_volume_fraction_name{"ow_vof"}; const std::string m_wave_negative_elevation_name{"ow_levelset"}; const std::string m_wave_velocity_name{"ow_velocity"}; - void convert_waves_to_blank_and_drag_flags(); }; } // namespace amr_wind::terraindrag From 21a862ad28bda45f3f14e7d71398a193c22f0552 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 22:17:28 -0700 Subject: [PATCH 23/31] address terrain roughness --- amr-wind/physics/TerrainDrag.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index d6201ddba1..7168f1fdfe 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -183,18 +183,18 @@ void TerrainDrag::post_regrid_actions() void TerrainDrag::convert_waves_to_blank_and_drag_flags() { const int nlevels = m_sim.repo().num_active_levels(); + // Uniform, low roughness for waves + m_terrainz0.setVal(1e-4); for (int level = 0; level < nlevels; ++level) { const auto geom = m_sim.repo().mesh().Geom(level); const auto& dx = geom.CellSizeArray(); const auto& prob_lo = geom.ProbLoArray(); auto& blanking = m_terrain_blank(level); - auto& terrainz0 = m_terrainz0(level); auto& terrain_height = m_terrain_height(level); auto& drag = m_terrain_drag(level); auto levelBlanking = blanking.arrays(); auto levelDrag = drag.arrays(); - auto levelz0 = terrainz0.arrays(); auto levelHeight = terrain_height.arrays(); const auto negative_wave_elevation = From 9b6fc957d7bd1b977f175b42b40fc17fb2e8adb8 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 22:18:55 -0700 Subject: [PATCH 24/31] ignore unused --- amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index f8dc181359..ebdd41e138 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -506,7 +506,7 @@ struct InitDataOp w2a_levelset.setVal(0.0); w2a_velocity.setVal(0.0); #else - amrex::ignore_unused(data, level, geom); + amrex::ignore_unused(data, level, geom, multiphase_mode); #endif } }; // namespace ocean_waves @@ -717,7 +717,7 @@ struct UpdateTargetFieldsOp }); } #else - amrex::ignore_unused(data); + amrex::ignore_unused(data, for_forcing_term); #endif } }; From 29db41bd04316e960a6a5969e7b2d78c74c11494 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 22:20:42 -0700 Subject: [PATCH 25/31] remove unused --- amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H | 3 +-- amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H index 1cb08c0efd..26f102bb99 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H @@ -88,10 +88,9 @@ struct UpdateTargetVolumeFractionOp< { BL_PROFILE("amr-wind::ocean_waves::OceanWaves::update_target_vof"); - const auto& wdata = data.meta(); auto& sim = data.sim(); - relaxation_zones::update_target_vof(sim, wdata); + relaxation_zones::update_target_vof(sim); } }; 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 e049736cb3..82f84a9728 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.cpp @@ -228,7 +228,7 @@ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata) density.fillpatch(time); } -void update_target_vof(CFDSim& sim, const RelaxZonesBaseData& wdata) +void update_target_vof(CFDSim& sim) { const int nlevels = sim.repo().num_active_levels(); auto& m_ow_levelset = sim.repo().get_field("ow_levelset"); From 07e6e420be751a6dd40e3d6e9a24a33a9c034b39 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 22:25:04 -0700 Subject: [PATCH 26/31] more for unused --- amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H index 26f102bb99..c939ff43ab 100644 --- a/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/relaxation_zones_ops.H @@ -28,7 +28,7 @@ void init_data_structures(CFDSim&); */ void apply_relaxation_zones(CFDSim& sim, const RelaxZonesBaseData& wdata); -void update_target_vof(CFDSim& sim, const RelaxZonesBaseData& wdata); +void update_target_vof(CFDSim& sim); void prepare_netcdf_file( const std::string& /*ncfile*/, From 1d1f4123613ba37c4ab158207844d1287395dd12 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 22:25:31 -0700 Subject: [PATCH 27/31] better function name --- amr-wind/physics/TerrainDrag.H | 4 ++-- amr-wind/physics/TerrainDrag.cpp | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/amr-wind/physics/TerrainDrag.H b/amr-wind/physics/TerrainDrag.H index d5349dcfc5..7e9e485cd5 100644 --- a/amr-wind/physics/TerrainDrag.H +++ b/amr-wind/physics/TerrainDrag.H @@ -36,8 +36,8 @@ public: //! Terrain Drag for waves - void convert_waves_to_blank_and_drag_flags(); - + void convert_waves_to_terrain_fields(); + std::string wave_velocity_field_name() const { return m_wave_velocity_name; diff --git a/amr-wind/physics/TerrainDrag.cpp b/amr-wind/physics/TerrainDrag.cpp index 7168f1fdfe..17a0cece59 100644 --- a/amr-wind/physics/TerrainDrag.cpp +++ b/amr-wind/physics/TerrainDrag.cpp @@ -156,7 +156,7 @@ void TerrainDrag::post_init_actions() return; } BL_PROFILE("amr-wind::" + this->identifier() + "::post_init_actions"); - convert_waves_to_blank_and_drag_flags(); + convert_waves_to_terrain_fields(); } void TerrainDrag::pre_advance_work() @@ -165,13 +165,13 @@ void TerrainDrag::pre_advance_work() return; } BL_PROFILE("amr-wind::" + this->identifier() + "::pre_advance_work"); - convert_waves_to_blank_and_drag_flags(); + convert_waves_to_terrain_fields(); } void TerrainDrag::post_regrid_actions() { if (m_terrain_is_waves) { - convert_waves_to_blank_and_drag_flags(); + convert_waves_to_terrain_fields(); } else { const int nlevels = m_sim.repo().num_active_levels(); for (int lev = 0; lev < nlevels; ++lev) { @@ -180,7 +180,7 @@ void TerrainDrag::post_regrid_actions() } } -void TerrainDrag::convert_waves_to_blank_and_drag_flags() +void TerrainDrag::convert_waves_to_terrain_fields() { const int nlevels = m_sim.repo().num_active_levels(); // Uniform, low roughness for waves From a7a411d63ff1a1f7dd8c64a3a7972998509761ed Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 22:45:18 -0700 Subject: [PATCH 28/31] change to warning, filling modes by default at beginning --- .../relaxation_zones/waves2amr_ops.H | 91 +++++++++---------- 1 file changed, 42 insertions(+), 49 deletions(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index ebdd41e138..dfeb31134e 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -294,7 +294,8 @@ struct ReadInputsOp if (std::abs(depth - (wdata.zsl - prob_lo_input[2])) > 1e-3 * depth) { amrex::Print() << "WARNING: Mismatch between water depths from AMR-Wind " - "domain and HOS data interpreted by Waves2AMR"; + "domain and HOS data interpreted by Waves2AMR. \n This " + "warning is not a concern when using waves as terrain.\n"; } // Allocate pointers for FFTW @@ -332,59 +333,51 @@ struct ReadInputsOp std::to_string(flag)); } - // If init_wave_field is activated and initialization will be done, get - // modes on every processor - if (wdata.init_wave_field && data.sim().time().time_index() == 0) { - bool no_EOF = wdata.rmodes.get_data( + // Get modes on every processor (unless requested to skip?) + bool no_EOF = wdata.rmodes.get_data( + wdata.ntime + wdata.n_offset, wdata.mX, wdata.mY, wdata.mZ, + wdata.mFS); + if (!no_EOF) { + // End of file detected, reset reading + wdata.n_offset = update_offset_timestep(wdata.ntime, wdata.n_winit); + // Print warning to screen + amrex::Print() + << "WARNING (waves2amr_ops): end of mode data file " + "detected, resetting to beginning of mode data.\n"; + // Read data again, now from a valid timestep + no_EOF = wdata.rmodes.get_data( wdata.ntime + wdata.n_offset, wdata.mX, wdata.mY, wdata.mZ, wdata.mFS); + // If no valid data is detected at this point, abort if (!no_EOF) { - // End of file detected, reset reading - wdata.n_offset = - update_offset_timestep(wdata.ntime, wdata.n_winit); - // Print warning to screen - amrex::Print() - << "WARNING (waves2amr_ops): end of mode data file " - "detected, resetting to beginning of mode data.\n"; - // Read data again, now from a valid timestep - no_EOF = wdata.rmodes.get_data( - wdata.ntime + wdata.n_offset, wdata.mX, wdata.mY, wdata.mZ, - wdata.mFS); - // If no valid data is detected at this point, abort - if (!no_EOF) { - amrex::Abort( - "waves2amr_ops: end of mode data file detected after " - "resetting to beginning; please evaluate HOS_init_time " - "or HOS_init_timestep and check the length of the mode " - "file."); - } + amrex::Abort( + "waves2amr_ops: end of mode data file detected after " + "resetting to beginning; please evaluate HOS_init_time " + "or HOS_init_timestep and check the length of the mode " + "file."); } + } - // Convert modes to spatial data - modes_hosgrid::copy_complex( - wdata.n0, wdata.n1, wdata.mFS, wdata.eta_mptr); - wdata.sp_eta_vec.resize( - static_cast(wdata.n0) * static_cast(wdata.n1), - 0.0); - modes_hosgrid::populate_hos_eta( - wdata.rmodes, wdata.plan, wdata.eta_mptr, wdata.sp_eta_vec); - // Mesh is not yet created, so get data at every height - const auto n_hts = wdata.hvec.size(); - wdata.sp_u_vec.resize( - static_cast(wdata.n0 * wdata.n1) * n_hts); - wdata.sp_v_vec.resize( - static_cast(wdata.n0 * wdata.n1) * n_hts); - wdata.sp_w_vec.resize( - static_cast(wdata.n0 * wdata.n1) * n_hts); - for (int iht = 0; iht < static_cast(n_hts); ++iht) { - // Get sample height - amrex::Real ht = wdata.hvec[iht]; - // Sample velocity - modes_hosgrid::populate_hos_vel( - wdata.rmodes, ht, wdata.mX, wdata.mY, wdata.mZ, wdata.plan, - wdata.u_mptr, wdata.v_mptr, wdata.w_mptr, wdata.sp_u_vec, - wdata.sp_v_vec, wdata.sp_w_vec, iht * wdata.n0 * wdata.n1); - } + // Convert modes to spatial data + modes_hosgrid::copy_complex( + wdata.n0, wdata.n1, wdata.mFS, wdata.eta_mptr); + wdata.sp_eta_vec.resize( + static_cast(wdata.n0) * static_cast(wdata.n1), 0.0); + modes_hosgrid::populate_hos_eta( + wdata.rmodes, wdata.plan, wdata.eta_mptr, wdata.sp_eta_vec); + // Mesh is not yet created, so get data at every height + const auto n_hts = wdata.hvec.size(); + wdata.sp_u_vec.resize(static_cast(wdata.n0 * wdata.n1) * n_hts); + wdata.sp_v_vec.resize(static_cast(wdata.n0 * wdata.n1) * n_hts); + wdata.sp_w_vec.resize(static_cast(wdata.n0 * wdata.n1) * n_hts); + for (int iht = 0; iht < static_cast(n_hts); ++iht) { + // Get sample height + amrex::Real ht = wdata.hvec[iht]; + // Sample velocity + modes_hosgrid::populate_hos_vel( + wdata.rmodes, ht, wdata.mX, wdata.mY, wdata.mZ, wdata.plan, + wdata.u_mptr, wdata.v_mptr, wdata.w_mptr, wdata.sp_u_vec, + wdata.sp_v_vec, wdata.sp_w_vec, iht * wdata.n0 * wdata.n1); } // Declare fields for HOS From b378762cb12ad2bdd745d8c4fbbcdefbc25fe988 Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 22:46:16 -0700 Subject: [PATCH 29/31] offset for waves in reg test --- test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp b/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp index 1e50d4bef8..9bbc0de2a1 100644 --- a/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp +++ b/test/test_files/abl_w2a_terrain/abl_w2a_terrain.inp @@ -65,9 +65,9 @@ OceanWaves.label = W2A1 OceanWaves.W2A1.type = W2AWaves OceanWaves.W2A1.HOS_modes_filename = ../ow_w2a/modes_HOS_SWENSE.dat # Offset to include wave troughs in domain -OceanWaves.W2A1.zero_sea_level = 0 +OceanWaves.W2A1.zero_sea_level = 10 # These variables should change with resolution in z -OceanWaves.W2A1.number_interp_points_in_z = 3 +OceanWaves.W2A1.number_interp_points_in_z = 6 OceanWaves.W2A1.interp_spacing_at_surface = 2. OceanWaves.W2A1.number_interp_above_surface = 3 From b66f9095646c2ded9da50b601d37864288ffed1d Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 23:28:42 -0700 Subject: [PATCH 30/31] tweak warning message --- amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index dfeb31134e..01e6b3a693 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -294,7 +294,7 @@ struct ReadInputsOp if (std::abs(depth - (wdata.zsl - prob_lo_input[2])) > 1e-3 * depth) { amrex::Print() << "WARNING: Mismatch between water depths from AMR-Wind " - "domain and HOS data interpreted by Waves2AMR. \n This " + "domain and HOS data interpreted by Waves2AMR. \n ^This " "warning is not a concern when using waves as terrain.\n"; } From 2ec52972f26dfd26c4cbb9a634192465a139ac8c Mon Sep 17 00:00:00 2001 From: Michael Kuhn Date: Sun, 24 Nov 2024 23:48:16 -0700 Subject: [PATCH 31/31] correct conditional to get ow_velocity correct --- amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H index 01e6b3a693..ace4308cb9 100644 --- a/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H +++ b/amr-wind/ocean_waves/relaxation_zones/waves2amr_ops.H @@ -90,7 +90,7 @@ void postprocess_velocity_mfab_liquid( // Set velocity to zero if no liquid present const amrex::Real cell_length_2D = std::sqrt(dx[0] * dx[0] + dx[2] * dx[2]); - if (phi[nbx](i, j, k) + cell_length_2D >= 0) { + if (phi[nbx](i, j, k) + cell_length_2D < 0) { vel[nbx](i, j, k, 0) = 0.0; vel[nbx](i, j, k, 1) = 0.0; vel[nbx](i, j, k, 2) = 0.0;