Skip to content

Commit

Permalink
Multiphase Hybrid - FY23 Q4 milestone (reinitialization/sharpening an…
Browse files Browse the repository at this point in the history 
…d pressure approach) (#1065)
  • Loading branch information
@mbkuhn
mbkuhn authored Jul 9, 2024
1 parent b8ab898 commit 41291fe
Show file tree
Hide file tree
Showing 19 changed files with 1,925 additions and 99 deletions.
4 changes: 4 additions & 0 deletions amr-wind/equation_systems/icns/icns_advection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,10 @@ MacProjOp::MacProjOp(
{
amrex::ParmParse pp("incflo");
pp.query("density", m_rho_0);
amrex::ParmParse pp_ovst("Overset");
bool disable_ovst_mac = false;
pp_ovst.query("disable_coupled_mac_proj", disable_ovst_mac);
m_has_overset = m_has_overset && !disable_ovst_mac;
}

void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept
Expand Down
149 changes: 105 additions & 44 deletions amr-wind/equation_systems/vof/volume_fractions.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,27 +16,27 @@ namespace amr_wind::multiphase {
* to zero.
*/

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal(
int i,
int j,
int k,
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_finite_difference_normal(
const int i,
const int j,
const int k,
amrex::Array4<amrex::Real const> const& volfrac,
amrex::Real& mx,
amrex::Real& my,
amrex::Real& mz) noexcept
{
amrex::Real mm1, mm2;

mm1 = volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
volfrac(i - 1, j + 1, k - 1) + volfrac(i - 1, j + 1, k + 1) +
2.0 * (volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1)) +
4.0 * volfrac(i - 1, j, k);
mm2 = volfrac(i + 1, j - 1, k - 1) + volfrac(i + 1, j - 1, k + 1) +
volfrac(i + 1, j + 1, k - 1) + volfrac(i + 1, j + 1, k + 1) +
2.0 * (volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1)) +
4.0 * volfrac(i + 1, j, k);
amrex::Real mm1 =
volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
volfrac(i - 1, j + 1, k - 1) + volfrac(i - 1, j + 1, k + 1) +
2.0 * (volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1)) +
4.0 * volfrac(i - 1, j, k);
amrex::Real mm2 =
volfrac(i + 1, j - 1, k - 1) + volfrac(i + 1, j - 1, k + 1) +
volfrac(i + 1, j + 1, k - 1) + volfrac(i + 1, j + 1, k + 1) +
2.0 * (volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1)) +
4.0 * volfrac(i + 1, j, k);
mx = mm1 - mm2;

mm1 = volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
Expand Down Expand Up @@ -64,25 +64,82 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal(
mz = mm1 - mm2;
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
youngs_finite_difference_normal_neumann(
const int i,
const int j,
const int k,
const int ibdy,
const int jbdy,
const int kbdy,
amrex::Array4<amrex::Real const> const& volfrac,
amrex::Real& mx,
amrex::Real& my,
amrex::Real& mz) noexcept
{
// Do neumann condition via indices
const int im1 = ibdy == -1 ? i : i - 1;
const int jm1 = jbdy == -1 ? j : j - 1;
const int km1 = kbdy == -1 ? k : k - 1;
const int ip1 = ibdy == +1 ? i : i + 1;
const int jp1 = jbdy == +1 ? j : j + 1;
const int kp1 = kbdy == +1 ? k : k + 1;

amrex::Real mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
2.0 * (volfrac(im1, jm1, k) + volfrac(im1, jp1, k) +
volfrac(im1, j, km1) + volfrac(im1, j, kp1)) +
4.0 * volfrac(im1, j, k);
amrex::Real mm2 = volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
2.0 * (volfrac(ip1, jm1, k) + volfrac(ip1, jp1, k) +
volfrac(ip1, j, km1) + volfrac(ip1, j, kp1)) +
4.0 * volfrac(ip1, j, k);
mx = mm1 - mm2;

mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
2.0 * (volfrac(im1, jm1, k) + volfrac(ip1, jm1, k) +
volfrac(i, jm1, km1) + volfrac(i, jm1, kp1)) +
4.0 * volfrac(i, jm1, k);
mm2 = volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
2.0 * (volfrac(im1, jp1, k) + volfrac(ip1, jp1, k) +
volfrac(i, jp1, km1) + volfrac(i, jp1, kp1)) +
4.0 * volfrac(i, jp1, k);
my = mm1 - mm2;

mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jp1, km1) +
volfrac(ip1, jm1, km1) + volfrac(ip1, jp1, km1) +
2.0 * (volfrac(im1, j, km1) + volfrac(ip1, j, km1) +
volfrac(i, jm1, km1) + volfrac(i, jp1, km1)) +
4.0 * volfrac(i, j, km1);
mm2 = volfrac(im1, jm1, kp1) + volfrac(im1, jp1, kp1) +
volfrac(ip1, jm1, kp1) + volfrac(ip1, jp1, kp1) +
2.0 * (volfrac(im1, j, kp1) + volfrac(ip1, j, kp1) +
volfrac(i, jm1, kp1) + volfrac(i, jp1, kp1)) +
4.0 * volfrac(i, j, kp1);
mz = mm1 - mm2;
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
int i,
int j,
int k,
const int i,
const int j,
const int k,
amrex::Array4<amrex::Real const> const& volfrac,
amrex::Real& mx,
amrex::Real& my,
amrex::Real& mz) noexcept
{
amrex::Array2D<amrex::Real, 0, AMREX_SPACEDIM + 1, 0, AMREX_SPACEDIM> m;
amrex::Real m1, m2;
// write the plane as: sgn(mx) X = my Y + mz Z + alpha
// m00 X = m01 Y + m02 Z + alpha
m1 = volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1) +
volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
volfrac(i - 1, j, k);
m2 = volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1) +
volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
volfrac(i + 1, j, k);
amrex::Real m1 = volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1) +
volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
volfrac(i - 1, j, k);
amrex::Real m2 = volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1) +
volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
volfrac(i + 1, j, k);

m(0, 0) = (m1 > m2) ? 1.0 : -1.0;

Expand Down Expand Up @@ -180,7 +237,8 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
}

// Youngs-CIAM scheme */
youngs_fd_normal(i, j, k, volfrac, m(3, 0), m(3, 1), m(3, 2));
youngs_finite_difference_normal(
i, j, k, volfrac, m(3, 0), m(3, 1), m(3, 2));
// normalize the set (mx,my,mz): |mx|+|my|+|mz| = 1
constexpr amrex::Real tiny = 1e-20;
t0 = std::abs(m(3, 0)) + std::abs(m(3, 1)) + std::abs(m(3, 2)) + tiny;
Expand Down Expand Up @@ -214,7 +272,10 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
* given that m1+m2+m3=1 (m1,m2,m3>0) and the volumetric fraction volF
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume_intercept(
amrex::Real b1, amrex::Real b2, amrex::Real b3, amrex::Real volF) noexcept
const amrex::Real b1,
const amrex::Real b2,
const amrex::Real b3,
const amrex::Real volF) noexcept
{
using namespace amrex;

Expand Down Expand Up @@ -296,15 +357,15 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume_intercept(
* (5) calculate volume (NOTE: it is assumed:s0=t0=0; ds0=dt0=1.)
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real cut_volume(
amrex::Real m1,
amrex::Real m2,
amrex::Real m3,
amrex::Real alpha,
amrex::Real r0,
amrex::Real dr0) noexcept
const amrex::Real m1,
const amrex::Real m2,
const amrex::Real m3,
const amrex::Real alpha,
const amrex::Real r0,
const amrex::Real dr0) noexcept
{

amrex::Real const_tiny = std::numeric_limits<amrex::Real>::epsilon();
const amrex::Real const_tiny = std::numeric_limits<amrex::Real>::epsilon();
amrex::Real al;

// move origin to x0
Expand Down Expand Up @@ -428,11 +489,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fit_plane(
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(
int i,
int j,
int k,
const int i,
const int j,
const int k,
amrex::Array4<amrex::Real const> const& volfrac,
int n_band = 1) noexcept
const int n_band = 1) noexcept
{
// n_band must be <= number of vof ghost cells (3)
constexpr amrex::Real tiny = 1e-12;
Expand All @@ -456,14 +517,14 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real levelset_to_vof(
int i,
int j,
int k,
amrex::Real eps,
const int i,
const int j,
const int k,
const amrex::Real eps,
amrex::Array4<amrex::Real const> const& phi) noexcept
{
amrex::Real mx, my, mz;
youngs_fd_normal(i, j, k, phi, mx, my, mz);
youngs_finite_difference_normal(i, j, k, phi, mx, my, mz);
mx = std::abs(mx / 32.);
my = std::abs(my / 32.);
mz = std::abs(mz / 32.);
Expand Down
3 changes: 3 additions & 0 deletions amr-wind/incflo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -319,6 +319,9 @@ void incflo::do_advance()
} else {
advance();
}
if (m_sim.has_overset()) {
m_ovst_ops.post_advance_work();
}
}

// Make a new level from scratch using provided BoxArray and
Expand Down
42 changes: 42 additions & 0 deletions amr-wind/overset/OversetOps.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,10 +13,52 @@ public:

void initialize(CFDSim& sim);
void pre_advance_work();
void post_advance_work();

void update_gradp();

private:
// Functions called within public functions
void parameter_output() const;
void sharpen_nalu_data();
void set_hydrostatic_gradp();
void replace_masked_gradp();

// Check for multiphase sim
bool m_vof_exists{false};
// Check for perturbational pressure
// (will be removed soon)
bool m_perturb_p{false};

// This is the only option for now
const bool m_use_hydrostatic_gradp{true};

// Coupling options
bool m_disable_nodal_proj{false};
bool m_disable_mac_proj{false};
bool m_replace_gradp_postsolve{false};

// Verbosity
int m_verbose{0};

// Reinitialization parameters
int m_n_iterations{10};
int m_calc_convg_interval{1}; // calctolniter, cconv
amrex::Real m_convg_tol = 1e-12;
amrex::Real m_relative_length_scale = 1.5;
amrex::Real m_upw_margin = 0.1;
amrex::Real m_target_cutoff = 0.0; // proc_tgvof_tol

// Tolerance for VOF-related bound checks
const amrex::Real m_vof_tol = 1e-12;
// Small number for approximate signed distance function
const amrex::Real m_asdf_tiny = 1e-12;

// Pointer for pressure gradient copy field
amr_wind::Field* m_gp_copy{nullptr};
// Pointer for MultiPhase physics
amr_wind::MultiPhase* m_mphase{nullptr};

CFDSim* m_sim_ptr;
};

Expand Down
Loading

0 comments on commit 41291fe

Please sign in to comment.