Skip to content

Commit

Permalink
Add psuedo-shock past partially refined sphere problem in 2D and 3D
Browse files Browse the repository at this point in the history
  • Loading branch information
cgilet committed Jul 12, 2023
1 parent 0582b05 commit 686f97d
Show file tree
Hide file tree
Showing 5 changed files with 309 additions and 7 deletions.
114 changes: 114 additions & 0 deletions Exec/eb_run2d/regtest.2d.shock_past_cylinder
Original file line number Diff line number Diff line change
@@ -0,0 +1,114 @@
#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# SIMULATION STOP #
#.......................................#
max_step = 10 # Max number of time steps
#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# SOLVER SETTINGS #
#.......................................#
ns.init_shrink = 1.0
#ns.init_dt = 0.01
ns.init_iter = 0
ns.do_init_proj = 0

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# Algorithm options
#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
ns.do_mom_diff=0
ns.do_cons_trac=0
ns.do_trac2 = 1
ns.do_cons_trac2=1

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# TIME STEP COMPUTATION #
#.......................................#
ns.cfl = 0.45

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# INPUT AND OUTPUT #
#.......................................#
amr.plot_int = 50 # Steps between plot files
amr.check_int = 1000 # Steps between checkpoint files

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# PHYSICS #
#.......................................#
ns.gravity = 0.

ns.vel_visc_coef = 0.
ns.scal_diff_coefs = 0. 0. #1 #0.03

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# ADAPTIVE MESH REFINEMENT #
#.......................................#
amr.n_cell = 64 32 # Grid cells at coarsest AMRlevel
amr.max_level = 1 # Max AMR level in hierarchy
# Refinement around the EB by default
#amr.n_error_buf = 8

# AMR tagging criteria
amr.refinement_indicators = box
amr.box.in_box_lo = 0.05 0.01 0.
amr.box.in_box_hi = 0.15 0.05 0.

ns.refine_cutcells = 0

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# GEOMETRY #
#.......................................#
geometry.prob_lo = -0.4 -0.2 -0.05 # Lo corner coordinates
geometry.prob_hi = 0.4 0.2 0.05 # Hi corner coordinates
geometry.is_periodic = 0 1 1 # Periodicity x y z (0/1)
geometry.coord_sys = 0

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# BOUNDARY CONDITIONS #
#.......................................#
# 0 = Interior/Periodic 3 = Symmetry
# 1 = Inflow 4 = SlipWall
# 2 = Outflow 5 = NoSlipWal
# Boundary conditions on the low end of the domain.
ns.lo_bc = 1 0 0
# Boundary conditions on the high end of the domain.
ns.hi_bc = 2 0 0

xlo.velocity = 1. 0. 0.
xlo.density = 10.
xlo.tracer = 1.
xlo.tracer2 = 0.

# Add cylinder
eb2.geom_type = sphere
eb2.sphere_radius = 0.08
eb2.sphere_center = 0.1 0.0 0.0
eb2.sphere_has_fluid_inside =0

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# PROBLEM PARAMETERS #
#.......................................#

prob.probtype = 3
prob.rho_1 = 10.
prob.rho_2 = 1.
prob.tra_1 = 1.
prob.tra_2 = 0.
prob.interface_width = 0.01
prob.blob_center = -0.25 0.0 0.
prob.blob_radius = 0.04
# Constant density initial condition
prob.density_ic = 1.0
# Set up a flow, defaults to zero
prob.velocity_ic = 1. 0. 0.

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# VERBOSITY #
#.......................................#
ns.v = 0 # NS solver
mac_proj.verbose = 0 # MacProjector
#diffusion.v = 1 # Diffusion operator
nodal_proj.verbose = 0 # Nodal projection

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# for testing #
#.......................................#
# turn tiling on for testing
fabarray.mfiter_tile_size = 8 8 8
114 changes: 114 additions & 0 deletions Exec/eb_run3d/regtest.3d.shock_past_sphere
Original file line number Diff line number Diff line change
@@ -0,0 +1,114 @@
#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# SIMULATION STOP #
#.......................................#
max_step = 10 # Max number of time steps
#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# SOLVER SETTINGS #
#.......................................#
ns.init_shrink = 1.0
#ns.init_dt = 0.01
ns.init_iter = 0
ns.do_init_proj = 0

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# Algorithm options
#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
ns.do_mom_diff=0
ns.do_cons_trac=0
ns.do_trac2 = 1
ns.do_cons_trac2=1

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# TIME STEP COMPUTATION #
#.......................................#
ns.cfl = 0.45

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# INPUT AND OUTPUT #
#.......................................#
amr.plot_int = 50 # Steps between plot files
amr.check_int = 1000 # Steps between checkpoint files

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# PHYSICS #
#.......................................#
ns.gravity = 0.

ns.vel_visc_coef = 0.
ns.scal_diff_coefs = 0. 0. #1 #0.03

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# ADAPTIVE MESH REFINEMENT #
#.......................................#
amr.n_cell = 64 32 32 # Grid cells at coarsest AMRlevel
amr.max_level = 1 # Max AMR level in hierarchy
# Refinement around the EB by default
#amr.n_error_buf = 8

# AMR tagging criteria
amr.refinement_indicators = box
amr.box.in_box_lo = 0.05 0.01 -0.05
amr.box.in_box_hi = 0.15 0.05 0.05

ns.refine_cutcells = 0

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# GEOMETRY #
#.......................................#
geometry.prob_lo = -0.4 -0.2 -0.2 # Lo corner coordinates
geometry.prob_hi = 0.4 0.2 0.2 # Hi corner coordinates
geometry.is_periodic = 0 1 1 # Periodicity x y z (0/1)
geometry.coord_sys = 0

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# BOUNDARY CONDITIONS #
#.......................................#
# 0 = Interior/Periodic 3 = Symmetry
# 1 = Inflow 4 = SlipWall
# 2 = Outflow 5 = NoSlipWal
# Boundary conditions on the low end of the domain.
ns.lo_bc = 1 0 0
# Boundary conditions on the high end of the domain.
ns.hi_bc = 2 0 0

xlo.velocity = 1. 0. 0.
xlo.density = 10.
xlo.tracer = 1.
xlo.tracer2 = 0.

# Add cylinder
eb2.geom_type = sphere
eb2.sphere_radius = 0.08
eb2.sphere_center = 0.1 0.0 0.0
eb2.sphere_has_fluid_inside =0

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# PROBLEM PARAMETERS #
#.......................................#

prob.probtype = 3
prob.rho_1 = 10.
prob.rho_2 = 1.
prob.tra_1 = 1.
prob.tra_2 = 0.
prob.interface_width = 0.01
prob.blob_center = -0.25 0.0 0.
prob.blob_radius = 0.04
# Constant density initial condition
prob.density_ic = 1.0
# Set up a flow, defaults to zero
prob.velocity_ic = 1. 0. 0.

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# VERBOSITY #
#.......................................#
ns.v = 0 # NS solver
mac_proj.verbose = 0 # MacProjector
#diffusion.v = 1 # Diffusion operator
nodal_proj.verbose = 0 # Nodal projection

#¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨#
# for testing #
#.......................................#
# turn tiling on for testing
fabarray.mfiter_tile_size = 8 8 8
8 changes: 1 addition & 7 deletions Source/NavierStokesBase.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1163,12 +1163,6 @@ NavierStokesBase::create_umac_grown (int nGrow,
// NOTE that this does not fill grid edges or corners.
//

#ifdef AMREX_USE_EB
if (!refine_cutcells) {
amrex::Abort("NSB::create_umac_grown not implemented for EB crossing the coarse-fine boundary.");
}
#endif

// Need 2 ghost cells here so we can safely check the status of all faces of a
// u_mac ghost cell
if (coarse_fine_mask == nullptr) {
Expand Down Expand Up @@ -1335,7 +1329,7 @@ NavierStokesBase::errorEst (TagBoxArray& tb,
if ( !ebfactory.isAllRegular() )
{
if (!refine_cutcells) {
amrex::Abort("For now, cutcells must always exist at finest level.");
amrex::Warning("Partially refined EB is still under development. This is not garanteed to work! Please use ns.refine_cutcells=1 for now.");
}

// Refine on cut cells
Expand Down
11 changes: 11 additions & 0 deletions Source/prob/prob_init.H
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,17 @@ static void init_constant_vel_rho (amrex::Box const& vbx,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi,
InitialConditions IC);

static void init_jump (amrex::Box const& vbx,
amrex::Array4<amrex::Real> const& press,
amrex::Array4<amrex::Real> const& vel,
amrex::Array4<amrex::Real> const& scal,
int nscal,
amrex::Box const& domain,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& dx,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& problo,
amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> const& probhi,
InitialConditions IC);

static void init_DoubleShearLayer (amrex::Box const& vbx,
amrex::Array4<amrex::Real> const& press,
amrex::Array4<amrex::Real> const& vel,
Expand Down
69 changes: 69 additions & 0 deletions Source/prob/prob_init.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,12 @@ void NavierStokes::prob_initData ()
S_new.array(mfi, Density), nscal,
domain, dx, problo, probhi, IC);
}
else if ( 3 == probtype )
{
init_jump(vbx, P_new.array(mfi), S_new.array(mfi, Xvel),
S_new.array(mfi, Density), nscal,
domain, dx, problo, probhi, IC);
}
else if ( 4 == probtype )
{
init_constant_vel_rho(vbx, P_new.array(mfi), S_new.array(mfi, Xvel),
Expand Down Expand Up @@ -274,6 +280,69 @@ void NavierStokes::init_constant_vel_rho (Box const& vbx,
});
}

void NavierStokes::init_jump (Box const& vbx,
Array4<Real> const& /*press*/,
Array4<Real> const& vel,
Array4<Real> const& scal,
const int nscal,
Box const& domain,
GpuArray<Real, AMREX_SPACEDIM> const& dx,
GpuArray<Real, AMREX_SPACEDIM> const& problo,
GpuArray<Real, AMREX_SPACEDIM> const& /*probhi*/,
InitialConditions IC)
{
const auto domlo = amrex::lbound(domain);

amrex::ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real x = problo[0] + (i - domlo.x + 0.5)*dx[0];
Real y = problo[1] + (j - domlo.y + 0.5)*dx[1];

//
// Fill Velocity
//
vel(i,j,k,0) = IC.v_x;
vel(i,j,k,1) = IC.v_y;

#if (AMREX_SPACEDIM == 3)
Real z = problo[2] + (k - domlo.z + 0.5)*dx[2];

vel(i,j,k,2) = IC.v_z;
#endif

Real dist = std::sqrt( (x-IC.blob_x)*(x-IC.blob_x)
+ (y-IC.blob_y)*(y-IC.blob_y)
#if (AMREX_SPACEDIM == 3)
+ (z-IC.blob_z)*(z-IC.blob_z)
#endif
);
//
// Scalars, ordered as Density, Tracer(s), Temp (if using)
//

Real x_jump = -0.1;

// Smooth the density interface just a bit
scal(i,j,k,0) = //dens1 * (2. + std::tanh(100.*(-.2-x)/IC.interface_width) ) / 2.;
IC.rho_1 + ((IC.rho_2-IC.rho_1)/2.0)*(1.0+std::tanh((x_jump-x)/IC.interface_width));

if (x <= x_jump ) {
//scal(i,j,k,0) = IC.rho_1;
scal(i,j,k,1) = IC.tra_1;
} else {
//scal(i,j,k,0) = IC.rho_2;
scal(i,j,k,1) = IC.tra_2;
}

for ( int nt=2; nt<nscal; nt++)
{
scal(i,j,k,nt) = 0.5*(1.0-std::tanh(25.*(dist-IC.blob_radius)/IC.interface_width));
// scal(i,j,k,nt) = dist < IC.blob_radius ? 1.0 : 0.0;
// scal(i,j,k,nt) = 0.0;
}
});
}

void NavierStokes::init_DoubleShearLayer (Box const& vbx,
Array4<Real> const& /*press*/,
Array4<Real> const& vel,
Expand Down

0 comments on commit 686f97d

Please sign in to comment.