Skip to content

Commit

Permalink
Refluxing2 (erf-model#1284)
Browse files Browse the repository at this point in the history 
* update particle functionality

* modify for refluxing plus some other cleanup such as terrain_type and coupling_type becoming enum structs

* clean up and add breaks in switch statement
  • Loading branch information
@asalmgren
asalmgren authored Nov 3, 2023
1 parent e8b228e commit 194cd60
Show file tree
Hide file tree
Showing 37 changed files with 498 additions and 823 deletions.
2 changes: 1 addition & 1 deletion Docs/sphinx_doc/Inputs.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -860,7 +860,7 @@ List of Parameters
| **erf.use_terrain** | use terrain-fitted | true / false | false |
| | coordinates? | | |
+-----------------------------+--------------------+--------------------+------------+
| **erf.terrain_type** | static or moving? | 0 / 1 | 0 |
| **erf.terrain_type** | static or moving? | Static / Moving | Static |
+-----------------------------+--------------------+--------------------+------------+
| **erf.terrain_smoothing** | specify terrain | 0, | 0 |
| | following | 1, | |
Expand Down
2 changes: 1 addition & 1 deletion Docs/sphinx_doc/RegressionTests.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ The following problems are currently tested in the CI:
| DensityCurrent_detJ2_MT | 256 4 64 | Symmetry | Periodic | SlipWall | None | use_terrain = true |
| | | Outflow | | SlipWall | | uses zlevels |
| | | Outflow | | SlipWall | | detJ = 2 everywhere |
| | | | | | | terrain_type = 1 |
| | | | | | | terrain_type = Moving |
+-------------------------------+----------+----------+----------+------------+-------+-----------------------+
| EkmanSpiral | 4 4 400 | Periodic | Periodic | NoSlipWall | Geo | +Coriolis |
| | | | | SlipWall | | +gravity |
Expand Down
2 changes: 1 addition & 1 deletion Docs/sphinx_doc/theory/WetEquations.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ The governing equations for this model are
\frac{\partial \rho_d}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u})
\frac{\partial (\rho_d \mathbf{u})}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \mathbf{u}) -
\frac{1}{1 + q_v + q_c} \nabla p^\prime - \nabla \cdot \tau + \mathbf{F} + \delta_{i,3}\mathbf{B}
\frac{1}{1 + q_v + q_c} ( \nabla p^\prime + \delta_{i,3}\mathbf{B} ) - \nabla \cdot \tau + \mathbf{F}
\frac{\partial (\rho_d \theta_d)}{\partial t} &= - \nabla \cdot (\rho_d \mathbf{u} \theta_d)
+ \nabla \cdot ( \rho_d \alpha_{T}\ \nabla \theta_d) + \frac{\theta_d L_v}{T_d C_{pd}} f_{cond}
Expand Down
4 changes: 2 additions & 2 deletions Exec/DevTests/MovingTerrain/inputs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@ zlo.type = "SlipWall"
zhi.type = "SlipWall"

# TERRRAIN GRID TYPE
erf.use_terrain = 1 # enable terrain stencils
erf.terrain_type = 1 # moving terrain
erf.use_terrain = true # enable terrain stencils
erf.terrain_type = Moving # moving terrain
erf.terrain_smoothing = 2 # Sullivan 2004 approach

# FOR NO SUBSTEPPING
Expand Down
5 changes: 3 additions & 2 deletions Exec/RegTests/ScalarAdvDiff/prob.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,8 +118,9 @@ Problem::init_custom_pert(
state(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz, r0) / r0));
} else if (parms.prob_type == 13) {
const Real r0_z = parms.rad_0 * (prob_hi[2] - prob_lo[2]);
const Real r2d_xz = std::sqrt((x-xc)*(x-xc)/(r0*r0) + (z-zc)*(z-zc)/(r0_z*r0_z)); //ellipse for mapfac shear validation
state(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz, r0_z) / r0_z));
//ellipse for mapfac shear validation
const Real r2d_xz_ell = std::sqrt((x-xc)*(x-xc)/(r0*r0) + (z-zc)*(z-zc)/(r0_z*r0_z));
state(i, j, k, RhoScalar_comp) = parms.A_0 * 0.25 * (1.0 + std::cos(PI * std::min(r2d_xz_ell, r0_z) / r0_z));
} else if (parms.prob_type == 14) {
state(i, j, k, RhoScalar_comp) = std::cos(PI*x);
} else {
Expand Down
11 changes: 7 additions & 4 deletions Source/Advection/Advection.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,21 +28,24 @@ void AdvectionSrcForRhoAndTheta (const amrex::Box& bx, const amrex::Box& valid_b
const amrex::Array4<const amrex::Real>& mf_u,
const amrex::Array4<const amrex::Real>& mf_v,
const AdvType horiz_adv_type, const AdvType vert_adv_type,
const int use_terrain);
const bool use_terrain);

/** Compute advection tendency for all scalars other than density and potential temperature */
void AdvectionSrcForScalars (const amrex::Box& bx,
void AdvectionSrcForScalars (int level, int finest_level,
const amrex::MFIter& mfi,
const amrex::Box& bx,
const int icomp, const int ncomp,
const amrex::Array4<const amrex::Real>& avg_xmom,
const amrex::Array4<const amrex::Real>& avg_ymom,
const amrex::Array4<const amrex::Real>& avg_zmom,
const amrex::Array4<const amrex::Real>& cell_prim,
const amrex::Array4<amrex::Real>& src,
const amrex::Array4<const amrex::Real>& detJ,
const amrex::Real dt,
const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM>& cellSize,
const amrex::Array4<const amrex::Real>& mf_m,
const AdvType horiz_adv_type, const AdvType vert_adv_type,
const int use_terrain);
const bool use_terrain, const bool is_two_way_coupling);

/** Compute advection tendencies for all components of momentum */
void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amrex::Box& bxz,
Expand All @@ -58,7 +61,7 @@ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amr
const amrex::Array4<const amrex::Real>& mf_u,
const amrex::Array4<const amrex::Real>& mf_v,
const AdvType horiz_adv_type, const AdvType vert_adv_type,
const int use_terrain, const int domhi_z);
const bool use_terrain, const int domhi_z);

AMREX_GPU_HOST_DEVICE
AMREX_FORCE_INLINE
Expand Down
2 changes: 1 addition & 1 deletion Source/Advection/AdvectionSrcForMom.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz,
const Array4<const Real>& mf_v,
const AdvType horiz_adv_type,
const AdvType vert_adv_type,
const int use_terrain,
const bool use_terrain,
const int domhi_z)
{
BL_PROFILE_VAR("AdvectionSrcForMom", AdvectionSrcForMom);
Expand Down
66 changes: 40 additions & 26 deletions Source/Advection/AdvectionSrcForMom_N.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,22 +42,27 @@ AdvectionSrcForXMom_N (int i, int j, int k,
amrex::Real interp_hi(0.), interp_lo(0.);

rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0));
rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0));
interp_u_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo);
interp_u_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
xflux_hi = rho_u_avg_hi * interp_hi;

rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0));
interp_u_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
xflux_lo = rho_u_avg_lo * interp_lo;

rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) * mf_v_inv(i,j+1,0) + rho_v(i-1, j+1, k) * mf_v_inv(i-1,j+1,0));
rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_v_inv(i,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0));
interp_u_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo);
interp_u_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
yflux_hi = rho_v_avg_hi * interp_hi;

rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_v_inv(i,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0));
interp_u_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
yflux_lo = rho_v_avg_lo * interp_lo;

rho_w_avg_hi = 0.5 * (rho_w(i, j, k+1) + rho_w(i-1, j, k+1));
rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i-1, j, k ));
interp_u_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi);
interp_u_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo);
interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
zflux_hi = rho_w_avg_hi * interp_hi;

rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i-1, j, k ));
interp_u_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo);
zflux_lo = rho_w_avg_lo * interp_lo;

amrex::Real mfsq = 1 / (mf_u_inv(i,j,0) * mf_u_inv(i,j,0));
Expand Down Expand Up @@ -110,22 +115,27 @@ AdvectionSrcForYMom_N (int i, int j, int k,
amrex::Real interp_hi(0.), interp_lo(0.);

rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j ,0) + rho_u(i+1, j-1, k) * mf_u_inv(i+1,j-1,0));
rho_u_avg_lo = 0.5 * (rho_u(i , j, k) * mf_u_inv(i ,j ,0) + rho_u(i , j-1, k) * mf_u_inv(i ,j-1,0));
interp_v_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo);
interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
xflux_hi = rho_u_avg_hi * interp_hi;

rho_u_avg_lo = 0.5 * (rho_u(i , j, k) * mf_u_inv(i ,j ,0) + rho_u(i , j-1, k) * mf_u_inv(i ,j-1,0));
interp_v_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
xflux_lo = rho_u_avg_lo * interp_lo;

rho_v_avg_hi = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j+1, k) * mf_v_inv(i ,j+1,0));
rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j-1, k) * mf_v_inv(i ,j-1,0));
interp_v_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo);
interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
yflux_hi = rho_v_avg_hi * interp_hi;

rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j-1, k) * mf_v_inv(i ,j-1,0));
interp_v_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
yflux_lo = rho_v_avg_lo * interp_lo;

rho_w_avg_hi = 0.5 * (rho_w(i, j, k+1) + rho_w(i, j-1, k+1));
rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i, j-1, k ));
interp_v_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi);
interp_v_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo);
interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
zflux_hi = rho_w_avg_hi * interp_hi;

rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i, j-1, k ));
interp_v_v.InterpolateInZ(i,j,k ,0,interp_lo,rho_w_avg_lo);
zflux_lo = rho_w_avg_lo * interp_lo;

amrex::Real mfsq = 1 / (mf_v_inv(i,j,0) * mf_v_inv(i,j,0));
Expand Down Expand Up @@ -186,15 +196,19 @@ AdvectionSrcForZMom_N (int i, int j, int k,
amrex::Real interp_hi(0.), interp_lo(0.);

rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) + rho_u(i+1, j, k-1)) * mf_u_inv(i+1,j ,0);
rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0);
interp_w_h.InterpolateInX(i,j,k,0,interp_hi,interp_lo,rho_u_avg_hi,rho_u_avg_lo);
interp_w_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
xflux_hi = rho_u_avg_hi * interp_hi;

rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0);
interp_w_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
xflux_lo = rho_u_avg_lo * interp_lo;

rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) + rho_v(i, j+1, k-1)) * mf_v_inv(i ,j+1,0);
rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0);
interp_w_h.InterpolateInY(i,j,k,0,interp_hi,interp_lo,rho_v_avg_hi,rho_v_avg_lo);
interp_w_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
yflux_hi = rho_v_avg_hi * interp_hi;

rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0);
interp_w_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
yflux_lo = rho_v_avg_lo * interp_lo;

// int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(domhi_z+1-k)), 2*(k+1));
Expand All @@ -208,15 +222,15 @@ AdvectionSrcForZMom_N (int i, int j, int k,
rho_w_avg_hi = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k+1));
if (k == domhi_z)
{
interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,AdvType::Centered_2nd);
interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_2nd);
} else if (k == domhi_z-1 || k == 1) {
if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,AdvType::Centered_4th);
interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_4th);
} else {
interp_w_wall.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi,vert_adv_type);
interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,vert_adv_type);
}
} else {
interp_w_v.InterpolateInZ_hi(i,j,k,0,interp_hi,rho_w_avg_hi);
interp_w_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
}
zflux_hi = rho_w_avg_hi * interp_hi;
}
Expand All @@ -231,15 +245,15 @@ AdvectionSrcForZMom_N (int i, int j, int k,
} else {
rho_w_avg_lo = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k-1));
if (k == 1) {
interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_2nd);
interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_2nd);
} else if (k == 2 || k == domhi_z) {
if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_4th);
interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_4th);
} else {
interp_w_wall.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo,vert_adv_type);
interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,vert_adv_type);
}
} else {
interp_w_v.InterpolateInZ_lo(i,j,k,0,interp_lo,rho_w_avg_lo);
interp_w_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo);
}
zflux_lo = rho_w_avg_lo * interp_lo;
}
Expand Down
Loading

0 comments on commit 194cd60

Please sign in to comment.