Skip to content

Commit

Permalink
get this working for fully general adaptive gridding -- this includes… (
Browse files Browse the repository at this point in the history 
erf-model#1750)

* get this working for fully general adaptive gridding -- this includes base state, terrain, moisture with two levels and factor 2 or 3 refinement

* update amrex

* fix box size
  • Loading branch information
@asalmgren
asalmgren authored Aug 21, 2024
1 parent 0535948 commit f10ae3d
Show file tree
Hide file tree
Showing 9 changed files with 124 additions and 76 deletions.
6 changes: 4 additions & 2 deletions Source/Diffusion/ComputeStrain_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,8 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
tbxxy.growLo(0,-1);
ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i, j-1, k)/mf_u(i,j-1,0))*dxInv[1] +
(-(8./3.) * v(i-1,j,k)/mf_v(i-1,j,0) + 3. * v(i,j,k)/mf_v(i,j,0) - (1./3.) * v(i+1,j,k)/mf_v(i+1,j,0))*dxInv[0] ) * mf_u(i,j,0)*mf_u(i,j,0);
(-(8./3.) * v(i-1,j,k)/mf_v(i-1,j,0) + 3. * v(i,j,k)/mf_v(i,j,0) -
(1./3.) * v(i+1,j,k)/mf_v(i+1,j,0))*dxInv[0] ) * mf_u(i,j,0)*mf_u(i,j,0);
});
}
if (xh_v_dir) {
Expand All @@ -104,7 +105,8 @@ ComputeStrain_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
tbxxy.growHi(0,-1);
ParallelFor(planexy,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i, j-1, k)/mf_u(i,j-1,0))*dxInv[1] +
-(-(8./3.) * v(i,j,k)/mf_v(i,j,0) + 3. * v(i-1,j,k)/mf_v(i-1,j,0) - (1./3.) * v(i-2,j,k)/mf_v(i-2,j,0))*dxInv[0] ) * mf_u(i,j,0)*mf_u(i,j,0);
-(-(8./3.) * v(i,j,k)/mf_v(i,j,0) + 3. * v(i-1,j,k)/mf_v(i-1,j,0) -
(1./3.) * v(i-2,j,k)/mf_v(i-2,j,0))*dxInv[0] ) * mf_u(i,j,0)*mf_u(i,j,0);
});
}

Expand Down
19 changes: 12 additions & 7 deletions Source/Diffusion/ComputeStrain_T.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,7 +98,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
(bc_ptr[BCVars::yvel_bc].hi(2) == ERFBCType::ext_dir_ingested) );
zh_v_dir = ( zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2)) );


//***********************************************************************************
// X-Dirichlet
//***********************************************************************************
if (xl_v_dir) {
Expand Down Expand Up @@ -181,6 +181,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
});
}

//***********************************************************************************
// Y-Dirichlet
//***********************************************************************************
if (yl_u_dir) {
Expand Down Expand Up @@ -263,6 +264,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
});
}

//***********************************************************************************
// Z-Dirichlet
//***********************************************************************************
if (zl_u_dir) {
Expand Down Expand Up @@ -312,7 +314,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
tau32(i,j,k) = tau23(i,j,k);
});
}
if (zh_v_dir) {
if (zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
tbxyz.growHi(2,-1);
ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Expand Down Expand Up @@ -368,10 +370,11 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
});
}

//***********************************************************************************
// Z-lo w/out Z-Dirichlet (GradWz extrapolation)
//***********************************************************************************
if (!zl_u_dir) {
Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) );
if (!zl_u_dir && (tbxxz.smallEnd(2) == domain_xz.smallEnd(2)) ) {
Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2));
tbxxz.growLo(2,-1);
ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Real GradWz = 0.5 * dxInv[2] * ( w(i ,j ,k+1) + w(i-1,j ,k+1)
Expand All @@ -387,7 +390,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
tau31(i,j,k) = tau13(i,j,k);
});
}
if (!zl_v_dir) {
if (!zl_v_dir && (tbxyz.smallEnd(2) == domain_yz.smallEnd(2))) {
Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) );
tbxyz.growLo(2,-1);
ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Expand All @@ -405,9 +408,10 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
});
}

//***********************************************************************************
// Z-hi w/out Z-Dirichlet (h_xi = h_eta = 0)
//***********************************************************************************
if (!zh_u_dir) {
if (!zh_u_dir && (tbxxz.bigEnd(2) == domain_xz.bigEnd(2))) {
Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) );
tbxxz.growHi(2,-1);
ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Expand All @@ -419,7 +423,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
tau31(i,j,k) = tau13(i,j,k);
});
}
if (!zh_v_dir) {
if (!zh_v_dir && (tbxyz.bigEnd(2) == domain_yz.bigEnd(2))) {
Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) );
tbxyz.growHi(2,-1);
ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept {
Expand All @@ -432,6 +436,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain,
});
}

//***********************************************************************************
// Fill the interior cells
//***********************************************************************************
// Cell centered strains
Expand Down
4 changes: 2 additions & 2 deletions Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22,
{
// NOTE: This gets us the lateral ghost cells for lev>0; which
// have been filled from FP Two Levels.
Box bxcc = mfi.growntilebox() & domain;
Box bxcc = mfi.growntilebox(1) & domain;

const Array4<Real>& mu_turb = eddyViscosity.array(mfi);
const Array4<Real>& hfx_x = Hfx1.array(mfi);
Expand Down Expand Up @@ -201,7 +201,7 @@ void ComputeTurbulentViscosityLES (const MultiFab& Tau11, const MultiFab& Tau22,
dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp)
- cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dzInv;
}
Real E = cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp);
Real E = amrex::max(cell_data(i,j,k,RhoKE_comp)/cell_data(i,j,k,Rho_comp),Real(0.0));
Real stratification = l_abs_g * dtheta_dz * l_inv_theta0; // stratification
Real length;
if (stratification <= eps) {
Expand Down
7 changes: 5 additions & 2 deletions Source/ERF.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -427,7 +427,9 @@ private:

void update_diffusive_arrays (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm);

void update_terrain_arrays (int lev, amrex::Real time);
void init_zphys (int lev, amrex::Real time);
void remake_zphys (int lev, amrex::Real time, std::unique_ptr<amrex::MultiFab>& temp_zphys_nd);
void update_terrain_arrays (int lev);

void Construct_ERFFillPatchers (int lev);

Expand All @@ -443,7 +445,8 @@ private:

void init_stuff (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
amrex::Vector<amrex::MultiFab>& lev_new, amrex::Vector<amrex::MultiFab>& lev_old,
amrex::MultiFab& tmp_base_state);
amrex::MultiFab& tmp_base_state,
std::unique_ptr<amrex::MultiFab>& tmp_zphys_nd);

// Initialize the Turbulent perturbation
void turbPert_update (const int lev, const amrex::Real dt);
Expand Down
47 changes: 39 additions & 8 deletions Source/ERF_make_new_arrays.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,8 @@ using namespace amrex;
void
ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
Vector<MultiFab>& lev_new, Vector<MultiFab>& lev_old,
MultiFab& tmp_base_state)
MultiFab& tmp_base_state,
std::unique_ptr<MultiFab>& tmp_zphys_nd)
{
if (lev == 0) {
min_k_at_level[lev] = 0;
Expand Down Expand Up @@ -80,10 +81,11 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,

// We need this to be one greater than the ghost cells to handle levels > 0
int ngrow = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 2;
z_phys_nd[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,1));
tmp_zphys_nd = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));

if (solverChoice.terrain_type != TerrainType::Static) {
z_phys_nd_new[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,1));
z_phys_nd_src[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,1));
z_phys_nd_new[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
z_phys_nd_src[lev] = std::make_unique<MultiFab>(ba_nd,dm,1,IntVect(ngrow,ngrow,ngrow));
}

} else {
Expand Down Expand Up @@ -437,7 +439,7 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap
}

if (l_use_kturb) {
eddyDiffs_lev[lev] = std::make_unique<MultiFab>( ba, dm, EddyDiff::NumDiffs, 1 );
eddyDiffs_lev[lev] = std::make_unique<MultiFab>(ba, dm, EddyDiff::NumDiffs, 2);
eddyDiffs_lev[lev]->setVal(0.0);
if(l_use_ddorf) {
SmnSmn_lev[lev] = std::make_unique<MultiFab>( ba, dm, 1, 0 );
Expand All @@ -451,10 +453,9 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap
}

void
ERF::update_terrain_arrays (int lev, Real time)
ERF::init_zphys (int lev, Real time)
{
if (solverChoice.use_terrain) {

//
// First interpolate from coarser level if there is one
//
Expand All @@ -477,7 +478,7 @@ ERF::update_terrain_arrays (int lev, Real time)
}

//
// NOTE: we are not currently doing this as described -- we will simply use
// NOTE: we are NOT currently doing this as described -- we will simply use
// the interpolation from coarse done above
// Then, if not using real/metgrid data,
// 1) redefine the terrain at k=0 for every fine box which includes k=0
Expand All @@ -500,7 +501,37 @@ ERF::update_terrain_arrays (int lev, Real time)
init_terrain_grid(lev,geom[lev],*z_phys_nd[lev],zlevels_stag,phys_bc_type);
} // init_type
} // lev == 0
}
}

void
ERF::remake_zphys (int lev, Real time, std::unique_ptr<MultiFab>& temp_zphys_nd)
{
if (solverChoice.use_terrain && lev > 0) {

Vector<MultiFab*> fmf = {z_phys_nd[lev].get(), z_phys_nd[lev].get()};
Vector<MultiFab*> cmf = {z_phys_nd[lev-1].get(), z_phys_nd[lev-1].get()};
Vector<Real> ftime = {time, time};
Vector<Real> ctime = {time, time};

PhysBCFunctNoOp null_bc;
Interpolater* mapper = &node_bilinear_interp;

FillPatchTwoLevels(*temp_zphys_nd, time,
cmf, ctime, fmf, ftime,
0, 0, 1, geom[lev-1], geom[lev],
null_bc, 0, null_bc, 0, refRatio(lev-1),
mapper, domain_bcs_type, 0);

std::swap(temp_zphys_nd, z_phys_nd[lev]);

} // use_terrain && lev > 0
}

void
ERF::update_terrain_arrays (int lev)
{
if (solverChoice.use_terrain) {
make_J(geom[lev],*z_phys_nd[lev],*detJ_cc[lev]);
make_areas(geom[lev],*z_phys_nd[lev],*ax[lev],*ay[lev],*az[lev]);
make_zcc(geom[lev],*z_phys_nd[lev],*z_phys_cc[lev]);
Expand Down
Loading

0 comments on commit f10ae3d

Please sign in to comment.