Skip to content

Commit

Permalink
Implement sponge_use_basestate for lateral boundaries
Browse files Browse the repository at this point in the history
  • Loading branch information
@ewquon
ewquon committed Sep 3, 2024
1 parent 83e7a36 commit 8ab2f86
Show file tree
Hide file tree
Showing 6 changed files with 105 additions and 48 deletions.
9 changes: 8 additions & 1 deletion Source/DataStructs/SpongeStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,10 +36,13 @@ struct SpongeChoice {
pp.query("zlo_sponge_end" , zlo_sponge_end);
pp.query("zhi_sponge_start", zhi_sponge_start);

pp.query("sponge_use_basestate", use_basestate); // alternative to specifying a constant sponge_density (and sponge_temperature)
pp.query("sponge_density" , sponge_density);
pp.query("sponge_x_velocity" , sponge_x_velocity);
pp.query("sponge_y_velocity" , sponge_y_velocity);
pp.query("sponge_z_velocity" , sponge_z_velocity);
pp.query("sponge_temperature", sponge_temperature);
pp.query("sponge_damp_temperature", damp_temperature);
}

void display()
Expand All @@ -63,6 +66,10 @@ struct SpongeChoice {
amrex::Real xlo_sponge_end, xhi_sponge_start;
amrex::Real ylo_sponge_end, yhi_sponge_start;
amrex::Real zlo_sponge_end, zhi_sponge_start;
amrex::Real sponge_density, sponge_x_velocity, sponge_y_velocity, sponge_z_velocity;

bool use_basestate = false;
bool damp_temperature = false;
amrex::Real sponge_density, sponge_temperature;
amrex::Real sponge_x_velocity, sponge_y_velocity, sponge_z_velocity;
};
#endif
88 changes: 59 additions & 29 deletions Source/SourceTerms/ERF_ApplySpongeZoneBCs.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,9 @@ ApplySpongeZoneBCsForCC (
const Geometry geom,
const Box& bx,
const Array4<Real>& cell_rhs,
const Array4<const Real>& cell_data)
const Array4<const Real>& cell_data,
const Array4<const Real>& rho0,
const Array4<const Real>& p0)
{
// Domain cell size and real bounds
auto dx = geom.CellSizeArray();
Expand All @@ -35,8 +37,12 @@ ApplySpongeZoneBCsForCC (
const Real zlo_sponge_end = spongeChoice.zlo_sponge_end;
const Real zhi_sponge_start = spongeChoice.zhi_sponge_start;

const bool damp_to_rhse = spongeChoice.use_basestate;
const Real sponge_density = spongeChoice.sponge_density;

const bool damp_temperature = spongeChoice.damp_temperature;
const Real sponge_temperature = spongeChoice.sponge_temperature;

// Domain valid box
const Box& domain = geom.Domain();
int domlo_x = domain.smallEnd(0);
Expand All @@ -63,48 +69,64 @@ ApplySpongeZoneBCsForCC (
Real y = ProbLoArr[1] + (jj+0.5) * dx[1];
Real z = ProbLoArr[2] + (kk+0.5) * dx[2];

Real rho_ref = (damp_to_rhse) ? rho0(i,j,k) : sponge_density;
Real rhotheta_ref = (damp_to_rhse) ? getRhoThetagivenP(p0(i,j,k))
: sponge_density*sponge_temperature;

// x left sponge
if(use_xlo_sponge_damping){
if (x < xlo_sponge_end) {
Real xi = (xlo_sponge_end - x) / (xlo_sponge_end - ProbLoArr[0]);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - sponge_density);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - rho_ref);
if (damp_temperature)
cell_rhs(i, j, k, 1) -= sponge_strength * xi * xi * (cell_data(i, j, k, 1) - rhotheta_ref);
}
}
// x right sponge
if(use_xhi_sponge_damping){
if (x > xhi_sponge_start) {
Real xi = (x - xhi_sponge_start) / (ProbHiArr[0] - xhi_sponge_start);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - sponge_density);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - rho_ref);
if (damp_temperature)
cell_rhs(i, j, k, 1) -= sponge_strength * xi * xi * (cell_data(i, j, k, 1) - rhotheta_ref);
}
}

// y left sponge
if(use_ylo_sponge_damping){
if (y < ylo_sponge_end) {
Real xi = (ylo_sponge_end - y) / (ylo_sponge_end - ProbLoArr[1]);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - sponge_density);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - rho_ref);
if (damp_temperature)
cell_rhs(i, j, k, 1) -= sponge_strength * xi * xi * (cell_data(i, j, k, 1) - rhotheta_ref);
}
}
// x right sponge
// y right sponge
if(use_yhi_sponge_damping){
if (y > yhi_sponge_start) {
Real xi = (y - yhi_sponge_start) / (ProbHiArr[1] - yhi_sponge_start);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - sponge_density);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - rho_ref);
if (damp_temperature)
cell_rhs(i, j, k, 1) -= sponge_strength * xi * xi * (cell_data(i, j, k, 1) - rhotheta_ref);
}
}

// x left sponge
// z left sponge
if(use_zlo_sponge_damping){
if (z < zlo_sponge_end) {
Real xi = (zlo_sponge_end - z) / (zlo_sponge_end - ProbLoArr[2]);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - sponge_density);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - rho_ref);
if (damp_temperature)
cell_rhs(i, j, k, 1) -= sponge_strength * xi * xi * (cell_data(i, j, k, 1) - rhotheta_ref);
}
}
// x right sponge
// z right sponge
if(use_zhi_sponge_damping){
if (z > zhi_sponge_start) {
Real xi = (z - zhi_sponge_start) / (ProbHiArr[2] - zhi_sponge_start);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - sponge_density);
cell_rhs(i, j, k, 0) -= sponge_strength * xi * xi * (cell_data(i, j, k, 0) - rho_ref);
if (damp_temperature)
cell_rhs(i, j, k, 1) -= sponge_strength * xi * xi * (cell_data(i, j, k, 1) - rhotheta_ref);
}
}
});
Expand All @@ -122,7 +144,8 @@ ApplySpongeZoneBCsForMom (
const Array4<Real>& rho_w_rhs,
const Array4<const Real>& rho_u,
const Array4<const Real>& rho_v,
const Array4<const Real>& rho_w)
const Array4<const Real>& rho_w,
const Array4<const Real>& rho0)
{
// Domain cell size and real bounds
auto dx = geom.CellSizeArray();
Expand All @@ -144,6 +167,7 @@ ApplySpongeZoneBCsForMom (
const Real zlo_sponge_end = spongeChoice.zlo_sponge_end;
const Real zhi_sponge_start = spongeChoice.zhi_sponge_start;

const bool damp_to_rhse = spongeChoice.use_basestate;
const Real sponge_density = spongeChoice.sponge_density;
const Real sponge_x_velocity = spongeChoice.sponge_x_velocity;
const Real sponge_y_velocity = spongeChoice.sponge_y_velocity;
Expand Down Expand Up @@ -175,41 +199,43 @@ ApplySpongeZoneBCsForMom (
Real y = ProbLoArr[1] + (jj+0.5) * dx[1];
Real z = ProbLoArr[2] + (kk+0.5) * dx[2];

Real rho_ref = (damp_to_rhse) ? rho0(i,j,k) : sponge_density;

// x lo sponge
if(use_xlo_sponge_damping){
if (x < xlo_sponge_end) {
Real xi = (xlo_sponge_end - x) / (xlo_sponge_end - ProbLoArr[0]);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - sponge_density*sponge_x_velocity);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - rho_ref*sponge_x_velocity);
}
}
// x hi sponge
if(use_xhi_sponge_damping){
if (x > xhi_sponge_start) {
Real xi = (x - xhi_sponge_start) / (ProbHiArr[0] - xhi_sponge_start);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - sponge_density*sponge_x_velocity);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - rho_ref*sponge_x_velocity);
}
}

// y lo sponge
if(use_ylo_sponge_damping){
if (y < ylo_sponge_end) {
Real xi = (ylo_sponge_end - y) / (ylo_sponge_end - ProbLoArr[1]);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - sponge_density*sponge_x_velocity);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - rho_ref*sponge_x_velocity);
}
}
// x right sponge
if(use_yhi_sponge_damping){
if (y > yhi_sponge_start) {
Real xi = (y - yhi_sponge_start) / (ProbHiArr[1] - yhi_sponge_start);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - sponge_density*sponge_x_velocity);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - rho_ref*sponge_x_velocity);
}
}

// z lo sponge
if(use_zlo_sponge_damping){
if (z < zlo_sponge_end) {
Real xi = (zlo_sponge_end - z) / (zlo_sponge_end - ProbLoArr[2]);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - sponge_density*sponge_x_velocity);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - rho_ref*sponge_x_velocity);
}
}

Expand All @@ -218,7 +244,7 @@ ApplySpongeZoneBCsForMom (
if(use_zhi_sponge_damping){
if (z > zhi_sponge_start) {
Real xi = (z - zhi_sponge_start) / (ProbHiArr[2] - zhi_sponge_start);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - sponge_density*sponge_x_velocity);
rho_u_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_u(i, j, k) - rho_ref*sponge_x_velocity);
}
}
});
Expand All @@ -234,41 +260,43 @@ ApplySpongeZoneBCsForMom (
Real y = ProbLoArr[1] + jj * dx[1];
Real z = ProbLoArr[2] + (kk+0.5) * dx[2];

Real rho_ref = (damp_to_rhse) ? rho0(i,j,k) : sponge_density;

// x lo sponge
if(use_xlo_sponge_damping){
if (x < xlo_sponge_end) {
Real xi = (xlo_sponge_end - x) / (xlo_sponge_end - ProbLoArr[0]);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - sponge_density*sponge_y_velocity);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - rho_ref*sponge_y_velocity);
}
}
// x hi sponge
if(use_xhi_sponge_damping){
if (x > xhi_sponge_start) {
Real xi = (x - xhi_sponge_start) / (ProbHiArr[0] - xhi_sponge_start);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - sponge_density*sponge_y_velocity);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - rho_ref*sponge_y_velocity);
}
}

// y lo sponge
if(use_ylo_sponge_damping){
if (y < ylo_sponge_end) {
Real xi = (ylo_sponge_end - y) / (ylo_sponge_end - ProbLoArr[1]);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - sponge_density*sponge_y_velocity);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - rho_ref*sponge_y_velocity);
}
}
// x right sponge
if(use_yhi_sponge_damping){
if (y > yhi_sponge_start) {
Real xi = (y - yhi_sponge_start) / (ProbHiArr[1] - yhi_sponge_start);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - sponge_density*sponge_y_velocity);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - rho_ref*sponge_y_velocity);
}
}

// z lo sponge
if(use_zlo_sponge_damping){
if (z < zlo_sponge_end) {
Real xi = (zlo_sponge_end - z) / (zlo_sponge_end - ProbLoArr[2]);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - sponge_density*sponge_y_velocity);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - rho_ref*sponge_y_velocity);
}
}

Expand All @@ -277,7 +305,7 @@ ApplySpongeZoneBCsForMom (
if(use_zhi_sponge_damping){
if (z > zhi_sponge_start) {
Real xi = (z - zhi_sponge_start) / (ProbHiArr[2] - zhi_sponge_start);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - sponge_density*sponge_y_velocity);
rho_v_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_v(i, j, k) - rho_ref*sponge_y_velocity);
}
}
});
Expand All @@ -293,41 +321,43 @@ ApplySpongeZoneBCsForMom (
Real y = ProbLoArr[1] + (jj+0.5) * dx[1];
Real z = ProbLoArr[2] + kk * dx[2];

Real rho_ref = (damp_to_rhse) ? rho0(i,j,k) : sponge_density;

// x left sponge
if(use_xlo_sponge_damping){
if (x < xlo_sponge_end) {
Real xi = (xlo_sponge_end - x) / (xlo_sponge_end - ProbLoArr[0]);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - sponge_density*sponge_z_velocity);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - rho_ref*sponge_z_velocity);
}
}
// x right sponge
if(use_xhi_sponge_damping){
if (x > xhi_sponge_start) {
Real xi = (x - xhi_sponge_start) / (ProbHiArr[0] - xhi_sponge_start);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - sponge_density*sponge_z_velocity);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - rho_ref*sponge_z_velocity);
}
}

// y lo sponge
if(use_ylo_sponge_damping){
if (y < ylo_sponge_end) {
Real xi = (ylo_sponge_end - y) / (ylo_sponge_end - ProbLoArr[1]);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - sponge_density*sponge_z_velocity);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - rho_ref*sponge_z_velocity);
}
}
// x right sponge
if(use_yhi_sponge_damping){
if (y > yhi_sponge_start) {
Real xi = (y - yhi_sponge_start) / (ProbHiArr[1] - yhi_sponge_start);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - sponge_density*sponge_z_velocity);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - rho_ref*sponge_z_velocity);
}
}

// z lo sponge
if(use_zlo_sponge_damping){
if (z < zlo_sponge_end) {
Real xi = (zlo_sponge_end - z) / (zlo_sponge_end - ProbLoArr[2]);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - sponge_density*sponge_z_velocity);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - rho_ref*sponge_z_velocity);
}
}

Expand All @@ -336,7 +366,7 @@ ApplySpongeZoneBCsForMom (
if(use_zhi_sponge_damping){
if (z > zhi_sponge_start) {
Real xi = (z - zhi_sponge_start) / (ProbHiArr[2] - zhi_sponge_start);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - sponge_density*sponge_z_velocity);
rho_w_rhs(i, j, k) -= sponge_strength * xi * xi * (rho_w(i, j, k) - rho_ref*sponge_z_velocity);
}
}
});
Expand Down
12 changes: 9 additions & 3 deletions Source/SourceTerms/ERF_make_mom_sources.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,6 +41,7 @@ void make_mom_sources (int /*level*/,
int /*nrk*/, Real dt, Real time,
Vector<MultiFab>& S_data,
const MultiFab & S_prim,
const MultiFab & base_state,
const MultiFab & /*xvel*/,
const MultiFab & /*yvel*/,
MultiFab & xmom_src,
Expand Down Expand Up @@ -205,6 +206,9 @@ void make_mom_sources (int /*level*/,
const Array4<const Real>& rho_v = S_data[IntVars::ymom].array(mfi);
const Array4<const Real>& rho_w = S_data[IntVars::zmom].array(mfi);

MultiFab r_hse(base_state, make_alias, 0, 1); // r_0 is first component
const Array4<const Real>& rho0_arr = r_hse.const_array(mfi);

const Array4< Real>& xmom_src_arr = xmom_src.array(mfi);
const Array4< Real>& ymom_src_arr = ymom_src.array(mfi);
const Array4< Real>& zmom_src_arr = zmom_src.array(mfi);
Expand Down Expand Up @@ -410,13 +414,15 @@ void make_mom_sources (int /*level*/,
// *****************************************************************************
if(solverChoice.spongeChoice.sponge_type == "input_sponge")
{
ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby, cell_data,
xmom_src_arr, ymom_src_arr, rho_u, rho_v, d_sponge_ptrs_at_lev);
ApplySpongeZoneBCsForMom_ReadFromFile(solverChoice.spongeChoice, geom, tbx, tby,
cell_data, xmom_src_arr, ymom_src_arr,
rho_u, rho_v, d_sponge_ptrs_at_lev);
}
else
{
ApplySpongeZoneBCsForMom(solverChoice.spongeChoice, geom, tbx, tby, tbz,
xmom_src_arr, ymom_src_arr, zmom_src_arr, rho_u, rho_v, rho_w);
xmom_src_arr, ymom_src_arr, zmom_src_arr,
rho_u, rho_v, rho_w, rho0_arr);
}

} // mfi
Expand Down
9 changes: 8 additions & 1 deletion Source/SourceTerms/ERF_make_sources.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ void make_sources (int level,
int /*nrk*/, Real dt, Real time,
Vector<MultiFab>& S_data,
const MultiFab & S_prim,
const MultiFab & base_state,
MultiFab & source,
#ifdef ERF_USE_RRTMGP
const MultiFab* qheating_rates,
Expand Down Expand Up @@ -182,6 +183,11 @@ void make_sources (int level,
const Array4<const Real> & cell_prim = S_prim.array(mfi);
const Array4<Real> & cell_src = source.array(mfi);

MultiFab r_hse(base_state, make_alias, 0, 1); // r_0 is first component
MultiFab p_hse(base_state, make_alias, 1, 1); // p_0 is second componentt
const Array4<const Real> & rho0_arr = r_hse.const_array(mfi);
const Array4<const Real> & p0_arr = p_hse.const_array(mfi);

#ifdef ERF_USE_RRTMGP
// *************************************************************************************
// 2. Add radiation source terms to (rho theta)
Expand Down Expand Up @@ -348,7 +354,8 @@ void make_sources (int level,
// 7. Add sponging
// *************************************************************************************
if(!(solverChoice.spongeChoice.sponge_type == "input_sponge")){
ApplySpongeZoneBCsForCC(solverChoice.spongeChoice, geom, bx, cell_src, cell_data);
ApplySpongeZoneBCsForCC(solverChoice.spongeChoice, geom, bx, cell_src, cell_data,
rho0_arr, p0_arr);
}

// *************************************************************************************
Expand Down
Loading

0 comments on commit 8ab2f86

Please sign in to comment.