From ff86d719c378be218670128a3261f25a796973be Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Tue, 3 Sep 2024 10:30:12 -0600 Subject: [PATCH] Improve initialization of zalesak case (#1223) * 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 --- amr-wind/physics/multiphase/ZalesakDisk.H | 4 +- amr-wind/physics/multiphase/ZalesakDisk.cpp | 59 ++++++++++++------- .../physics/multiphase/ZalesakDiskScalarVel.H | 4 +- .../multiphase/ZalesakDiskScalarVel.cpp | 59 ++++++++++++------- .../zalesak_disk_godunov.inp | 11 +--- .../zalesak_disk_scalar_vel.inp | 3 +- 6 files changed, 80 insertions(+), 60 deletions(-) diff --git a/amr-wind/physics/multiphase/ZalesakDisk.H b/amr-wind/physics/multiphase/ZalesakDisk.H index 05cdc04b8e..48a8f72baa 100644 --- a/amr-wind/physics/multiphase/ZalesakDisk.H +++ b/amr-wind/physics/multiphase/ZalesakDisk.H @@ -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}; diff --git a/amr-wind/physics/multiphase/ZalesakDisk.cpp b/amr-wind/physics/multiphase/ZalesakDisk.cpp index fee7728c23..d3377c0b7f 100644 --- a/amr-wind/physics/multiphase/ZalesakDisk.cpp +++ b/amr-wind/physics/multiphase/ZalesakDisk.cpp @@ -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) { @@ -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) { diff --git a/amr-wind/physics/multiphase/ZalesakDiskScalarVel.H b/amr-wind/physics/multiphase/ZalesakDiskScalarVel.H index aca98cacd2..3f756dfdcf 100644 --- a/amr-wind/physics/multiphase/ZalesakDiskScalarVel.H +++ b/amr-wind/physics/multiphase/ZalesakDiskScalarVel.H @@ -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}; diff --git a/amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp b/amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp index 83593438b2..b4c5e2b570 100644 --- a/amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp +++ b/amr-wind/physics/multiphase/ZalesakDiskScalarVel.cpp @@ -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) { @@ -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) { diff --git a/test/test_files/zalesak_disk_godunov/zalesak_disk_godunov.inp b/test/test_files/zalesak_disk_godunov/zalesak_disk_godunov.inp index 13a05c750d..0d195e7f0d 100644 --- a/test/test_files/zalesak_disk_godunov/zalesak_disk_godunov.inp +++ b/test/test_files/zalesak_disk_godunov/zalesak_disk_godunov.inp @@ -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 #¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨# diff --git a/test/test_files/zalesak_disk_scalar_vel/zalesak_disk_scalar_vel.inp b/test/test_files/zalesak_disk_scalar_vel/zalesak_disk_scalar_vel.inp index 1a96641c50..923a35081c 100644 --- a/test/test_files/zalesak_disk_scalar_vel/zalesak_disk_scalar_vel.inp +++ b/test/test_files/zalesak_disk_scalar_vel/zalesak_disk_scalar_vel.inp @@ -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