Skip to content

Commit

Permalink
use cc project and face-centered SF
Browse files Browse the repository at this point in the history
  • Loading branch information
Hua Tan committed Oct 1, 2024
1 parent 1cbede4 commit c8e1cd3
Show file tree
Hide file tree
Showing 15 changed files with 88,898 additions and 54 deletions.
5 changes: 3 additions & 2 deletions src/convection/incflo_compute_MAC_projected_velocities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,8 +178,9 @@ incflo::compute_MAC_projected_velocities (
allow_inflow_on_outflow, BC_MF.get());

//add surface tension
//if(m_vof_advect_tracer)
// get_volume_of_fluid ()->velocity_face_source(lev,l_dt, AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]));
if(m_vof_advect_tracer && m_use_cc_proj)
get_volume_of_fluid ()->velocity_face_source(lev,0.5*l_dt, AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]),
AMREX_D_DECL(nullptr, nullptr, nullptr));


if(0){
Expand Down
11 changes: 7 additions & 4 deletions src/incflo.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,14 +103,16 @@ public:
amrex::Vector<amrex::MultiFab const*> const& density,
amrex::Vector<amrex::MultiFab const*> const& tracer_old,
amrex::Vector<amrex::MultiFab const*> const& tracer_new,
bool include_pressure_gradient = true);
bool include_pressure_gradient = true,
bool include_SF = false);
void compute_vel_forces_on_level ( int lev,
amrex::MultiFab& vel_forces,
const amrex::MultiFab& velocity,
const amrex::MultiFab& density,
const amrex::MultiFab& tracer_old,
const amrex::MultiFab& tracer_new,
bool include_pressure_gradient = true);
bool include_pressure_gradient = true,
bool include_SF = false);


///////////////////////////////////////////////////////////////////////////
Expand Down Expand Up @@ -225,8 +227,8 @@ public:
AMREX_D_DECL(amrex::Vector<amrex::MultiFab const*> const& u_mac,
amrex::Vector<amrex::MultiFab const*> const& v_mac,
amrex::Vector<amrex::MultiFab const*> const& w_mac));
void update_vof_density (amrex::Vector<amrex::MultiFab*> const& density,
amrex::Vector<amrex::MultiFab*> const& tracer);
void update_vof_density (int lev, amrex::Vector<amrex::MultiFab*> const& density,
amrex::Vector<amrex::MultiFab*> const& tracer);


[[nodiscard]] amrex::Array<amrex::MultiFab,AMREX_SPACEDIM>
Expand Down Expand Up @@ -662,6 +664,7 @@ private:

// cell-centered pressure gradient
amrex::MultiFab gp;
amrex::MultiFab gp_mac;

amrex::MultiFab conv_velocity;
amrex::MultiFab conv_velocity_o;
Expand Down
3 changes: 2 additions & 1 deletion src/incflo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,8 @@ void incflo::InitData ()
}

InitialIterations();

//get_volume_of_fluid()->WriteTecPlotFile (m_cur_time,m_nstep);
//amrex::Abort("finish initial projection");
// Set m_nstep to 0 before entering time loop
m_nstep = 0;

Expand Down
3 changes: 2 additions & 1 deletion src/incflo_apply_predictor.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,7 +128,8 @@ void incflo::ApplyPredictor (bool incremental_projection)
//when VOF method is used to advect the tracer, density and viscosity of each cell will
// depend the VOF field value of the cell.
if (m_vof_advect_tracer)
update_vof_density (get_density_old(),get_tracer_old());
for (int lev = 0; lev <= finest_level; ++lev)
update_vof_density (lev, get_density_old(),get_tracer_old());

// *************************************************************************************
// Compute explicit viscous term
Expand Down
14 changes: 8 additions & 6 deletions src/incflo_compute_forces.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,11 +40,12 @@ void incflo::compute_vel_forces (Vector<MultiFab*> const& vel_forces,
Vector<MultiFab const*> const& density,
Vector<MultiFab const*> const& tracer_old,
Vector<MultiFab const*> const& tracer_new,
bool include_pressure_gradient)
bool include_pressure_gradient,
bool include_SF)
{
for (int lev = 0; lev <= finest_level; ++lev)
compute_vel_forces_on_level (lev, *vel_forces[lev], *velocity[lev], *density[lev],
*tracer_old[lev], *tracer_new[lev], include_pressure_gradient);
*tracer_old[lev], *tracer_new[lev], include_pressure_gradient, include_SF);
}

void incflo::compute_vel_forces_on_level (int lev,
Expand All @@ -53,7 +54,8 @@ void incflo::compute_vel_forces_on_level (int lev,
const MultiFab& density,
const MultiFab& tracer_old,
const MultiFab& tracer_new,
bool include_pressure_gradient)
bool include_pressure_gradient,
bool include_SF)
{
GpuArray<Real,3> l_gravity{m_gravity[0],m_gravity[1],m_gravity[2]};
GpuArray<Real,3> l_gp0{m_gp0[0], m_gp0[1], m_gp0[2]};
Expand Down Expand Up @@ -157,7 +159,7 @@ void incflo::compute_vel_forces_on_level (int lev,
// rho: density

//fixme: we just consider the surface tension for first tracer
if (m_vof_advect_tracer && m_sigma[0]!=0.){
if (m_vof_advect_tracer && m_sigma[0]!=0.&&!m_use_cc_proj&&include_SF){

VolumeOfFluid* vof_p = get_volume_of_fluid ();

Expand Down Expand Up @@ -528,8 +530,8 @@ static int oct[3][2] = { { 1, 2 }, { 0, 2 }, { 0, 1 } };

}
}
else if(choice==4) {
//cell-centered surface tension force
else if (choice ==4) {
//cell-centered surface tension force
#ifdef _OPENMP
#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
Expand Down
7 changes: 7 additions & 0 deletions src/prob/prob_init_fluid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -207,6 +207,13 @@ void incflo::prob_init_fluid (int lev)
if (1109 == m_probtype) {
get_volume_of_fluid ()->tracer_vof_init_fraction(lev, ld.tracer);
MultiFab::Copy(ld.tracer_o, ld.tracer, 0, 0, 1, ld.tracer.nGrow());
ld.tracer_o.FillBoundary(geom[lev].periodicity());
if (m_vof_advect_tracer){
update_vof_density (lev, get_density_new(),get_tracer_new());
MultiFab::Copy(ld.density_o, ld.density, 0, 0, 1, ld.density.nGrow());
fillpatch_density(lev, m_t_new[lev], ld.density_o, 3);
}

}
}

Expand Down
66 changes: 55 additions & 11 deletions src/projection/incflo_apply_cc_projection.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -297,7 +297,7 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,
LPInfo lp_info;
lp_info.setMaxCoarseningLevel(m_mac_mg_max_coarsening_level);
#ifndef AMREX_USE_EB
if (m_constant_density) {
if (m_constant_density&&!m_vof_advect_tracer) {
Vector<BoxArray> ba;
Vector<DistributionMapping> dm;
for (auto const& ir : inv_rho) {
Expand Down Expand Up @@ -340,13 +340,14 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,
mac_vec[lev][2] = w_mac[lev];);
}

Vector<MultiFab> sfu_mac(finest_level+1), sfv_mac(finest_level+1), sfw_mac(finest_level+1);
// Compute velocity on faces
for (int lev = 0; lev <= finest_level; ++lev)
{
// Predict normal velocity to faces -- note that the {u_mac, v_mac, w_mac}
// returned from this call are on face CENTROIDS
vel[lev]->FillBoundary(geom[lev].periodicity());
#if 1
#if 0
MOL::ExtrapVelToFaces(*vel[lev],
AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]),
geom[lev],
Expand All @@ -355,10 +356,27 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,
#else
average_ccvel_to_mac( mac_vec[lev], *vel[lev]);
#endif


//add surface tension
AMREX_D_TERM(sfu_mac[lev].define(u_mac[lev]->boxArray(), dmap[lev], 1, u_mac[lev]->nGrow(), MFInfo(), Factory(lev));,
sfv_mac[lev].define(v_mac[lev]->boxArray(), dmap[lev], 1, v_mac[lev]->nGrow(), MFInfo(), Factory(lev));,
sfw_mac[lev].define(w_mac[lev]->boxArray(), dmap[lev], 1, w_mac[lev]->nGrow(), MFInfo(), Factory(lev)););

AMREX_D_TERM(sfu_mac[lev].setVal(0.0);,
sfv_mac[lev].setVal(0.0);,
sfw_mac[lev].setVal(0.0););


if(m_vof_advect_tracer)
get_volume_of_fluid ()->velocity_face_source(lev, scaling_factor, AMREX_D_DECL(*u_mac[lev], *v_mac[lev], *w_mac[lev]),
AMREX_D_DECL(&sfu_mac[lev], &sfv_mac[lev], &sfw_mac[lev]));

}

macproj->setUMAC(mac_vec);


if (m_verbose > 2) amrex::Print() << "CC Projection:\n";
//
// Perform MAC projection: - del dot (dt/rho) grad phi = div(U)
Expand Down Expand Up @@ -391,12 +409,22 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,

for (int lev=0; lev <= finest_level; ++lev)
{
#ifdef AMREX_USE_EB
amrex::Abort("Haven't written mac_to_ccvel for EB");
#else
//#ifdef AMREX_USE_EB
// amrex::Abort("Haven't written mac_to_ccvel for EB");
//#else
average_mac_to_ccvel(GetArrOfPtrs(m_fluxes[lev]),*cc_gphi[lev]);
#endif
//#endif
}
// computer the cell-centered surface tension term (see note in VolumeOfFluid:: velocity_face_source)
VolumeOfFluid* vof_p = get_volume_of_fluid ();
if(m_vof_advect_tracer)
for (int lev=0; lev <= finest_level; ++lev)
{
AMREX_D_TERM(Copy(m_fluxes[lev][0],sfu_mac[lev], 0, 0, 1, 0);,
Copy(m_fluxes[lev][1],sfv_mac[lev], 0, 0, 1, 0);,
Copy(m_fluxes[lev][2],sfw_mac[lev], 0, 0, 1, 0););
average_mac_to_ccvel(GetArrOfPtrs(m_fluxes[lev]),vof_p->force[lev]);
}

for(int lev = 0; lev <= finest_level; lev++)
{
Expand All @@ -413,55 +441,71 @@ void incflo::ApplyCCProjection (Vector<MultiFab const*> density,

Array4<Real> const& u = ld.velocity.array(mfi);
Array4<Real const> const& rho = density[lev]->const_array(mfi);
Array4<Real const> const& gsf = vof_p->force[lev].const_array(mfi);

Real r0 = m_ro_0;

amrex::ParallelFor(tbx, [u,gphi,p_cc,phi,incremental] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
AMREX_D_TERM(u(i,j,k,0) += gphi(i,j,k,0);,
u(i,j,k,1) += gphi(i,j,k,1);,
u(i,j,k,2) += gphi(i,j,k,2););
//we need to add the surface-tension effect to the cell-centered velocity
if(m_vof_advect_tracer){
AMREX_D_TERM(u(i,j,k,0) -= gsf(i,j,k,0)*scaling_factor;,
u(i,j,k,1) -= gsf(i,j,k,1)*scaling_factor;,
u(i,j,k,2) -= gsf(i,j,k,2)*scaling_factor;);
}
if (incremental)
p_cc (i,j,k) += phi(i,j,k);
else
p_cc (i,j,k) = phi(i,j,k);
});

if (incremental && m_constant_density) {
if (incremental && m_constant_density&&!m_vof_advect_tracer) {
amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
AMREX_D_TERM(gp_cc(i,j,k,0) -= gphi(i,j,k,0) * r0 / scaling_factor;,
gp_cc(i,j,k,1) -= gphi(i,j,k,1) * r0 / scaling_factor;,
gp_cc(i,j,k,2) -= gphi(i,j,k,2) * r0 / scaling_factor;);

});
} else if (incremental && !m_constant_density) {
} else if (incremental && (!m_constant_density||m_vof_advect_tracer)) {
amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
AMREX_D_TERM(gp_cc(i,j,k,0) -= gphi(i,j,k,0) * rho(i,j,k) / scaling_factor;,
gp_cc(i,j,k,1) -= gphi(i,j,k,1) * rho(i,j,k) / scaling_factor;,
gp_cc(i,j,k,2) -= gphi(i,j,k,2) * rho(i,j,k) / scaling_factor;);
});
} else if (!incremental && m_constant_density) {
} else if (!incremental && m_constant_density&&!m_vof_advect_tracer) {
amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
AMREX_D_TERM(gp_cc(i,j,k,0) = -gphi(i,j,k,0) * r0 / scaling_factor;,
gp_cc(i,j,k,1) = -gphi(i,j,k,1) * r0 / scaling_factor;,
gp_cc(i,j,k,2) = -gphi(i,j,k,2) * r0 / scaling_factor;);
});
} else if (!incremental && !m_constant_density) {
} else if (!incremental && (!m_constant_density||m_vof_advect_tracer)) {
amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
AMREX_D_TERM(gp_cc(i,j,k,0) = -gphi(i,j,k,0) * rho(i,j,k) / scaling_factor;,
gp_cc(i,j,k,1) = -gphi(i,j,k,1) * rho(i,j,k) / scaling_factor;,
gp_cc(i,j,k,2) = -gphi(i,j,k,2) * rho(i,j,k) / scaling_factor;);
if(m_vof_advect_tracer){
AMREX_D_TERM(gp_cc(i,j,k,0) += gsf(i,j,k,0) * rho(i,j,k);,
gp_cc(i,j,k,1) += gsf(i,j,k,1) * rho(i,j,k);,
gp_cc(i,j,k,2) += gsf(i,j,k,2) * rho(i,j,k););
}

});
}
}
ld.gp.FillBoundary(geom[lev].periodicity());
ld.p_cc.FillBoundary(geom[lev].periodicity());
}

//get_volume_of_fluid()->WriteTecPlotFile (m_cur_time,m_nstep);
//amrex::Abort("finish initial projection");

// ***************************************************************************************
// END OF MAC STUFF
// ***************************************************************************************
Expand Down
5 changes: 4 additions & 1 deletion src/vof/VolumeOfFluid.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,9 +20,12 @@ public:
void write_tecplot_surface(amrex::Real time, int nstep);
void WriteTecPlotFile (amrex::Real time, int nstep);
void output_droplet (amrex::Real time, int nstep);
int domain_tag_droplets (int finest_level, amrex::Vector<amrex::BoxArray > const &grids, amrex::Vector<amrex::Geometry > const& geom,
amrex::Vector<amrex::MultiFab*> const& vof,amrex::Vector<amrex::MultiFab*> const& tag);
void apply_velocity_field(amrex::Real time, int nstep);
void velocity_face_source(int lev,amrex::Real dt, AMREX_D_DECL(amrex::MultiFab& u_mac, amrex::MultiFab& v_mac,
amrex::MultiFab& w_mac));
amrex::MultiFab& w_mac),
AMREX_D_DECL(amrex::MultiFab* sfu_mac, amrex::MultiFab* sfv_mac, amrex::MultiFab* sfw_mac));

// normal vector of interface
amrex::Vector<amrex::MultiFab> normal;
Expand Down
Loading

0 comments on commit c8e1cd3

Please sign in to comment.