Skip to content

Commit

Permalink
Improve initialization of zalesak case (Exawind#1223)
Browse files Browse the repository at this point in the history 
* fix resolution-dependent issues, make levelset field way smoother
* apply the same changes to the ZalesakDiskScalarVel case
* update some things in reg test files because there will be diffs anyway
  • Loading branch information
@mbkuhn @moprak-nrel
mbkuhn authored and moprak-nrel committed Sep 6, 2024
1 parent e92c074 commit ff86d71
Show file tree
Hide file tree
Showing 6 changed files with 80 additions and 60 deletions.
4 changes: 2 additions & 2 deletions amr-wind/physics/multiphase/ZalesakDisk.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,8 +49,8 @@ private:
//! sphere radius value
amrex::Real m_radius{0.16};

//! slot width value
amrex::Real m_width{0.04};
//! slot half width value
amrex::Real m_halfwidth{0.04};

//! slot depth
amrex::Real m_depth{0.2};
Expand Down
59 changes: 37 additions & 22 deletions amr-wind/physics/multiphase/ZalesakDisk.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ void ZalesakDisk::initialize_fields(int level, const amrex::Geometry& geom)
const amrex::Real zc = m_loc[2];
const amrex::Real radius = m_radius;
const amrex::Real TT = m_TT;
const amrex::Real width = m_width;
const amrex::Real hwidth = m_halfwidth;
const amrex::Real depth = m_depth;

for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
Expand All @@ -71,33 +71,48 @@ void ZalesakDisk::initialize_fields(int level, const amrex::Geometry& geom)
vel(i, j, k, 0) = 2.0 * M_PI / TT * (0.5 - y);
vel(i, j, k, 1) = 2.0 * M_PI / TT * (x - 0.5);
vel(i, j, k, 2) = 0.0;
// First define the sphere

phi(i, j, k) =
radius - std::sqrt(
(x - xc) * (x - xc) + (y - yc) * (y - yc) +
(z - zc) * (z - zc));
// then the slot
amrex::Real eps = std::cbrt(dx[0] * dx[1] * dx[2]);
if (y - yc <= radius && y - yc >= radius - depth &&
std::abs(x - xc) <= width &&
std::sqrt(
(x - xc) * (x - xc) + (y - yc) * (y - yc) +
(z - zc) * (z - zc)) < radius + eps) {
amrex::Real d1;
if (x > xc) {
d1 = std::abs(xc + width - x);
// First define the sphere
const amrex::Real r = std::sqrt(
(x - xc) * (x - xc) + (y - yc) * (y - yc) +
(z - zc) * (z - zc));
phi(i, j, k) = radius - r;

// Then the slot
// Signed distances in lateral (x, y) directions
const amrex::Real sd_xr = -hwidth + (x - xc);
const amrex::Real sd_xl = -hwidth - (x - xc);
const amrex::Real sd_x = amrex::max(sd_xr, sd_xl);

const amrex::Real sd_y = radius - depth - (y - yc);
const amrex::Real min_signed_dist = amrex::max(sd_x, sd_y);

// Additional distance if past sphere (distance to corners)
const amrex::Real reduced_radius =
std::sqrt(radius * radius - hwidth * hwidth);
const amrex::Real r_2D =
std::sqrt(std::pow(y - yc, 2) + std::pow(z - zc, 2));
const amrex::Real sd_r = -std::sqrt(
std::pow(r_2D - reduced_radius, 2) + std::pow(sd_x, 2));

const bool in_slot_x_ymin =
y - yc > radius - depth && std::abs(x - xc) < hwidth;
const bool in_slot_r = r_2D < reduced_radius;

if (in_slot_x_ymin) {
// Prescribe slot distances directly (overwrite sphere)
if (in_slot_r) {
phi(i, j, k) = min_signed_dist;
} else {
d1 = std::abs(xc - width - x);
phi(i, j, k) = sd_r;
}
const amrex::Real d2 = std::abs(y - (yc + radius - depth));
const amrex::Real min_dist = amrex::min(d1, d2);

phi(i, j, k) = -min_dist;
} else {
// Select the minimum of the two
phi(i, j, k) = amrex::min(phi(i, j, k), min_signed_dist);
}

amrex::Real smooth_heaviside;
eps = std::cbrt(2. * dx[0] * dx[1] * dx[2]);
const amrex::Real eps = std::cbrt(2. * dx[0] * dx[1] * dx[2]);
if (phi(i, j, k) > eps) {
smooth_heaviside = 1.0;
} else if (phi(i, j, k) < -eps) {
Expand Down
4 changes: 2 additions & 2 deletions amr-wind/physics/multiphase/ZalesakDiskScalarVel.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,8 +68,8 @@ private:
//! sphere radius value
amrex::Real m_radius{0.16};

//! slot width value
amrex::Real m_width{0.04};
//! slot half-width value
amrex::Real m_halfwidth{0.04};

//! slot depth
amrex::Real m_depth{0.2};
Expand Down
59 changes: 37 additions & 22 deletions amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ void ZalesakDiskScalarVel::initialize_fields(
const amrex::Real zc = m_loc[2];
const amrex::Real radius = m_radius;
const amrex::Real TT = m_TT;
const amrex::Real width = m_width;
const amrex::Real hwidth = m_halfwidth;
const amrex::Real depth = m_depth;

for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
Expand All @@ -103,33 +103,48 @@ void ZalesakDiskScalarVel::initialize_fields(

vel(i, j, k, 1) = 0.0;
vel(i, j, k, 2) = 0.0;
// First define the sphere

phi(i, j, k) =
radius - std::sqrt(
(x - xc) * (x - xc) + (y - yc) * (y - yc) +
(z - zc) * (z - zc));
// then the slot
amrex::Real eps = std::cbrt(dx[0] * dx[1] * dx[2]);
if (y - yc <= radius && y - yc >= radius - depth &&
std::abs(x - xc) <= width &&
std::sqrt(
(x - xc) * (x - xc) + (y - yc) * (y - yc) +
(z - zc) * (z - zc)) < radius + eps) {
amrex::Real d1;
if (x > xc) {
d1 = std::abs(xc + width - x);
// First define the sphere
const amrex::Real r = std::sqrt(
(x - xc) * (x - xc) + (y - yc) * (y - yc) +
(z - zc) * (z - zc));
phi(i, j, k) = radius - r;

// Then the slot
// Signed distances in lateral (x, y) directions
const amrex::Real sd_xr = -hwidth + (x - xc);
const amrex::Real sd_xl = -hwidth - (x - xc);
const amrex::Real sd_x = amrex::max(sd_xr, sd_xl);

const amrex::Real sd_y = radius - depth - (y - yc);
const amrex::Real min_signed_dist = amrex::max(sd_x, sd_y);

// Additional distance if past sphere (distance to corners)
const amrex::Real reduced_radius =
std::sqrt(radius * radius - hwidth * hwidth);
const amrex::Real r_2D =
std::sqrt(std::pow(y - yc, 2) + std::pow(z - zc, 2));
const amrex::Real sd_r = -std::sqrt(
std::pow(r_2D - reduced_radius, 2) + std::pow(sd_x, 2));

const bool in_slot_x_ymin =
y - yc > radius - depth && std::abs(x - xc) < hwidth;
const bool in_slot_r = r_2D < reduced_radius;

if (in_slot_x_ymin) {
// Prescribe slot distances directly (overwrite sphere)
if (in_slot_r) {
phi(i, j, k) = min_signed_dist;
} else {
d1 = std::abs(xc - width - x);
phi(i, j, k) = sd_r;
}
const amrex::Real d2 = std::abs(y - (yc + radius - depth));
const amrex::Real min_dist = amrex::min(d1, d2);

phi(i, j, k) = -min_dist;
} else {
// Select the minimum of the two
phi(i, j, k) = amrex::min(phi(i, j, k), min_signed_dist);
}

amrex::Real smooth_heaviside;
eps = std::cbrt(2. * dx[0] * dx[1] * dx[2]);
const amrex::Real eps = std::cbrt(2. * dx[0] * dx[1] * dx[2]);
if (phi(i, j, k) > eps) {
smooth_heaviside = 1.0;
} else if (phi(i, j, k) < -eps) {
Expand Down
11 changes: 1 addition & 10 deletions test/test_files/zalesak_disk_godunov/zalesak_disk_godunov.inp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,22 +15,13 @@ time.cfl = 0.95 # CFL factor
#.......................................#
time.plot_interval = 10 # Steps between plot files
time.checkpoint_interval = -100 # Steps between checkpoint files
io.output_default_fields = 0
io.outputs = density velocity velocity_mueff vof
io.outputs = density velocity levelset vof

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# PHYSICS #
#.......................................#
incflo.initial_iterations = 0
incflo.do_initial_proj = 0
incflo.use_godunov = 1
incflo.godunov_type = "weno"
transport.model = TwoPhaseTransport
transport.viscosity_fluid1=0.0
transport.viscosity_fluid2=0.0
transport.laminar_prandtl = 0.7
transport.turbulent_prandtl = 0.3333
turbulence.model = Laminar

incflo.physics = MultiPhase ZalesakDisk
#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
Expand Down
3 changes: 1 addition & 2 deletions test/test_files/zalesak_disk_scalar_vel/zalesak_disk_scalar_vel.inp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,7 @@ io.outputs = density velocity vof
#.......................................#
incflo.initial_iterations = 0
incflo.do_initial_proj = 0
incflo.use_godunov = 1
incflo.godunov_type = "weno"
incflo.godunov_type = "weno_z"
incflo.mflux_type = "minmod"

MultiPhase.density_fluid1 = 1e3
Expand Down

0 comments on commit ff86d71

Please sign in to comment.