Skip to content

Commit

Permalink
changes to include pressure sharpening and treat dx differently
Browse files Browse the repository at this point in the history
  • Loading branch information
@mbkuhn
mbkuhn committed Jun 11, 2024
1 parent 9f05c8d commit 453ec3c
Show file tree
Hide file tree
Showing 3 changed files with 297 additions and 48 deletions.
72 changes: 52 additions & 20 deletions amr-wind/overset/OversetOps.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -199,27 +199,51 @@ void OversetOps::sharpen_nalu_data()
auto& levelset = repo.get_field("levelset");
auto& rho = repo.get_field("density");
auto& velocity = repo.get_field("velocity");
auto& gp_noghost = repo.get_field("gp");
auto& p = repo.get_field("p");

// Create scratch fields for fluxes
// 5 components are vof, density, and 3 of velocity
auto flux_x = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::XFACE);
auto flux_y = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::YFACE);
auto flux_z = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::ZFACE);
// Create scratch field for approximate signed distance function and grad
// (components 0-2 are gradient, 3 is asdf)
// 9 components are vof, density, 3 of velocity, 3 of gp, and psource flag
auto flux_x = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::XFACE);
auto flux_y = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::YFACE);
auto flux_z = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::ZFACE);
// Create scratch field for pressure source term
auto p_src = repo.create_scratch_field(1, 0, amr_wind::FieldLoc::NODE);
// Create scratch fields for normal vector and target vof field
auto normal_vec = repo.create_scratch_field(3, vof.num_grow()[0] - 1);

auto target_vof = repo.create_scratch_field(1, vof.num_grow()[0]);

// Create levelset field
// Create scratch field for pressure gradient that has ghost cells
// Sharpening fluxes (at faces) have 1 ghost, requiring fields to have >= 2
auto gp_scr = repo.create_scratch_field(3, 2);
auto& gp = *gp_scr;

// Give initial max possible value of pseudo-velocity scale
auto dx_lev0 = (geom[0]).CellSizeArray();
const amrex::Real max_pvscale =
std::min(std::min(dx_lev0[0], dx_lev0[1]), dx_lev0[2]);
amrex::Real pvscale = max_pvscale;

// Prep things that do not change with iterations
for (int lev = 0; lev < nlevels; ++lev) {
// Thickness used here is user parameter, whatever works best
auto dx = (geom[lev]).CellSizeArray();
const amrex::Real i_th = m_rlscale * std::cbrt(dx[0] * dx[1] * dx[2]);

// Populate approximate signed distance function
overset_ops::populate_psi(levelset(lev), vof(lev), i_th, m_asdf_tiny);

// Populate gp scratch field
gp(lev).setVal(0.0); // for external boundaries
amrex::MultiFab::Copy(gp(lev), gp_noghost(lev), 0, 0, 3, 0); // nonghost
gp(lev).FillBoundary(geom[lev].periodicity()); // internal

// Get pseudo-velocity scale, proportional to smallest dx in iblank
const amrex::Real pvscale_lev =
overset_ops::calculate_pseudo_velocity_scale(
iblank_cell(lev), dx, max_pvscale);
pvscale = std::min(pvscale, pvscale_lev);
}
amrex::ParallelDescriptor::ReduceRealMin(pvscale);

// Convert levelset to vof to get target_vof
m_mphase->levelset2vof(iblank_cell, *target_vof);
Expand Down Expand Up @@ -270,13 +294,14 @@ void OversetOps::sharpen_nalu_data()
// Sharpening fluxes for vof, density, and momentum
overset_ops::populate_sharpen_fluxes(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
(*target_vof)(lev), (*normal_vec)(lev), velocity(lev), m_margin,
m_mphase->rho1(), m_mphase->rho2());
(*target_vof)(lev), (*normal_vec)(lev), velocity(lev), gp(lev),
rho(lev), pvscale, m_margin, m_mphase->rho1(),
m_mphase->rho2());

// Process fluxes
overset_ops::process_fluxes(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev),
iblank_cell(lev));
overset_ops::process_fluxes_calc_src(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
(*p_src)(lev), iblank_cell(lev));

// Measure convergence to determine if loop can stop
if (cconv) {
Expand Down Expand Up @@ -312,15 +337,19 @@ void OversetOps::sharpen_nalu_data()

// Apply fluxes
for (int lev = 0; lev < nlevels; ++lev) {
auto dx = (geom[lev]).CellSizeArray();

overset_ops::apply_fluxes(
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
rho(lev), velocity(lev), ptfac, m_vof_tol);
(*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), (*p_src)(lev),
vof(lev), rho(lev), velocity(lev), gp(lev), p(lev), dx, ptfac,
m_vof_tol);

// Update ghost cells
vof(lev).FillBoundary(geom[lev].periodicity());
velocity(lev).FillBoundary(geom[lev].periodicity());
gp(lev).FillBoundary(geom[lev].periodicity());
}

// Fillpatch for ghost cells
vof.fillpatch((*m_sim_ptr).time().current_time());
velocity.fillpatch((*m_sim_ptr).time().current_time());

// Update density (fillpatch built in)
m_mphase->set_density_via_vof();

Expand All @@ -335,6 +364,9 @@ void OversetOps::sharpen_nalu_data()
}
}

// Fillpatch for pressure to make sure pressure stencil has all points
p.fillpatch((*m_sim_ptr).time().current_time());

// Purely for debugging via visualization, should be removed later
// Currently set up to overwrite the levelset field (not used as time
// evolves) with the post-sharpening velocity magnitude
Expand Down
47 changes: 47 additions & 0 deletions amr-wind/overset/overset_ops_K.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,6 +114,53 @@ void AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE velocity_face(
}
}

void AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE gp_rho_face(
int i,
int j,
int k,
int dir,
amrex::Array4<amrex::Real const> const& vof,
amrex::Array4<amrex::Real const> const& gp,
amrex::Array4<amrex::Real const> const& rho,
amrex::Real& uface,
amrex::Real& vface,
amrex::Real& wface)
{
// Set up neighbor indices
int ii = i;
int jj = j;
int kk = k;
ii += (dir == 0) ? -1 : 0;
jj += (dir == 1) ? -1 : 0;
kk += (dir == 2) ? -1 : 0;

// Gradient of phi normal to interface
const amrex::Real gphi = (vof(i, j, k) - vof(ii, jj, kk));

// Get velocities on both sides
const amrex::Real u_ = gp(i, j, k, 0) / rho(i, j, k);
const amrex::Real v_ = gp(i, j, k, 1) / rho(i, j, k);
const amrex::Real w_ = gp(i, j, k, 2) / rho(i, j, k);
const amrex::Real u_nb = gp(ii, jj, kk, 0) / rho(ii, jj, kk);
const amrex::Real v_nb = gp(ii, jj, kk, 1) / rho(ii, jj, kk);
const amrex::Real w_nb = gp(ii, jj, kk, 2) / rho(ii, jj, kk);
// Average value when gphi = 0
uface = 0.5 * (u_ + u_nb);
vface = 0.5 * (v_ + v_nb);
wface = 0.5 * (w_ + w_nb);
// Use simple upwinding
if (gphi > 0.0) {
uface = u_nb;
vface = v_nb;
wface = w_nb;
}
if (gphi < 0.0) {
uface = u_;
vface = v_;
wface = w_;
}
}

} // namespace amr_wind::overset_ops

#endif
Loading

0 comments on commit 453ec3c

Please sign in to comment.