Skip to content

Commit

Permalink
Align waves2amr initialization with other wave profiles; add +abl reg…
Browse files Browse the repository at this point in the history 
… test (#1368)

* documentation (for arguments now in reg test)
* replaced MFIters in ocean waves
  • Loading branch information
@mbkuhn
mbkuhn authored Nov 22, 2024
1 parent 3191948 commit d665147
Show file tree
Hide file tree
Showing 8 changed files with 522 additions and 407 deletions.
230 changes: 111 additions & 119 deletions amr-wind/ocean_waves/relaxation_zones/linear_waves_ops.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,71 +45,68 @@ struct InitDataOp<LinearWaves>
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];
}
});
}
};

Expand All @@ -135,61 +132,56 @@ struct UpdateRelaxZonesOp<LinearWaves>
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.;
}
});
}
}
};
Expand Down
Loading

0 comments on commit d665147

Please sign in to comment.