Skip to content

Commit

Permalink
use th0 in dry anelastic buoyancy (erf-model#1981)
Browse files Browse the repository at this point in the history 
* use th0 in dry anelastic buoyancy

* remove unused

* don't average down density when doing anelastic
  • Loading branch information
@asalmgren
asalmgren authored Nov 24, 2024
1 parent a8e96e8 commit d63c6ae
Show file tree
Hide file tree
Showing 5 changed files with 105 additions and 69 deletions.
9 changes: 8 additions & 1 deletion Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -494,7 +494,14 @@ ERF::post_timestep (int nstep, Real time, Real dt_lev0)
// We need to do this before anything else because refluxing changes the
// values of coarse cells underneath fine grids with the assumption they'll
// be over-written by averaging down
AverageDownTo(lev,0,ncomp);
int src_comp;
if (solverChoice.anelastic[lev]) {
src_comp = 1;
} else {
src_comp = 0;
}
int num_comp = ncomp - src_comp;
AverageDownTo(lev,src_comp,num_comp);
}
}

Expand Down
7 changes: 5 additions & 2 deletions Source/IO/ERF_Checkpoint.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -425,13 +425,16 @@ ERF::ReadCheckpointFile ()
}
MultiFab base(grids[lev],dmap[lev],ncomp_base,ng_base);
VisMF::Read(base, MultiFabFileFullPrefix(lev, restart_chkfile, "Level_", "BaseState"));

MultiFab::Copy(base_state[lev],base,0,0,ncomp_base,ng_base);

if (read_old_base_state) {
for (MFIter mfi(base_state[lev],TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.growntilebox(1);
// We only compute theta_0 on valid cells since we will impose domain BC's after restart
const Box& bx = mfi.tilebox();
Array4<Real> const& fab = base_state[lev].array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
fab(i,j,k,BaseState::th0_comp) = getRhoThetagivenP(fab(i,j,k,BaseState::p0_comp))
/ fab(i,j,k,BaseState::r0_comp);
Expand Down
131 changes: 73 additions & 58 deletions Source/SourceTerms/ERF_buoyancy_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,61 +7,32 @@
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_anelastic (int& i,
int& j,
int& k,
buoyancy_dry_anelastic (int& i, int& j, int& k,
amrex::Real const& grav_gpu,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
// Note: this is the same term as the moist anelastic buoyancy when qv = qc = qt = 0
amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1);
amrex::Real theta_d_hi = cell_data(i,j,k ,RhoTheta_comp)/r0_arr(i,j,k);
amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp);
amrex::Real theta_d_hi = cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp);

amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_hi);

amrex::Real theta_d_0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)) / r0_arr(i,j,k-1);
amrex::Real theta_d_0_hi = getRhoThetagivenP(p0_arr(i,j,k )) / r0_arr(i,j,k );
amrex::Real theta_d_0_lo = th0_arr(i,j,k-1);
amrex::Real theta_d_0_hi = th0_arr(i,j,k );

amrex::Real theta_d_0_wface = amrex::Real(0.5) * (theta_d_0_lo + theta_d_0_hi);

return (-grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface * 0.5 * (r0_arr(i,j,k) + r0_arr(i,j,k-1)));
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_anelastic_T (int& i,
int& j,
int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k));
amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp);
amrex::Real t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp);
amrex::Real qplus = (t_hi-t0_hi)/t0_hi;
amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1));

amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1));
amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp);
amrex::Real t_lo = getTgivenPandTh(p0_arr(i,j,k-1), cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1), rd_over_cp);
amrex::Real qminus = (t_lo-t0_lo)/t0_lo;

amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
return (-r0_q_avg * grav_gpu);
return (-rho0_wface * grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface);
}


AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_moist_anelastic (int& i,
int& j,
int& k,
buoyancy_moist_anelastic (int& i, int& j, int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rv_over_rd,
const amrex::Array4<const amrex::Real>& r0_arr,
Expand All @@ -87,15 +58,15 @@ buoyancy_moist_anelastic (int& i,

amrex::Real theta_v_0_wface = amrex::Real(0.5) * (theta_v_0_lo + theta_v_0_hi);

return (-grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface * 0.5 * (r0_arr(i,j,k) + r0_arr(i,j,k-1)));
amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1));

return (-rho0_wface * grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface);
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_default (int& i,
int& j,
int& k,
buoyancy_dry_default (int& i, int& j, int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& r0_arr,
Expand All @@ -118,13 +89,11 @@ buoyancy_dry_default (int& i,
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_type1 (int& i,
int& j,
int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
buoyancy_rhopert (int& i, int& j, int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real rhop_hi = cell_data(i,j,k ,Rho_comp) - r0_arr(i,j,k );
amrex::Real rhop_lo = cell_data(i,j,k-1,Rho_comp) - r0_arr(i,j,k-1);
Expand All @@ -139,9 +108,7 @@ buoyancy_type1 (int& i,
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_type2 (int& i,
int& j,
int& k,
buoyancy_type2 (int& i, int& j, int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
amrex::Real* rho_d_ptr,
Expand Down Expand Up @@ -191,9 +158,7 @@ buoyancy_type2 (int& i,
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_type3 (int& i,
int& j,
int& k,
buoyancy_type3 (int& i, int& j, int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
amrex::Real* rho_d_ptr,
Expand Down Expand Up @@ -235,9 +200,7 @@ buoyancy_type3 (int& i,
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_type4 (int& i,
int& j,
int& k,
buoyancy_type4 (int& i, int& j, int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
amrex::Real* rho_d_ptr,
Expand Down Expand Up @@ -272,4 +235,56 @@ buoyancy_type4 (int& i,

return (-qavg * r0avg * grav_gpu);
}

// **************************************************************************************
// Routines below this line are not currently used
// **************************************************************************************

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_Tpert (int& i, int& j, int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& th0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp);
amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp));
amrex::Real qplus = (t_hi-t0_hi)/t0_hi;

amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp);
amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp));
amrex::Real qminus = (t_lo-t0_lo)/t0_lo;

amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
return (-r0_q_avg * grav_gpu);
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_anelastic_T (int& i, int& j, int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k));
amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp);
amrex::Real t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp);
amrex::Real qplus = (t_hi-t0_hi)/t0_hi;

amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1));
amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp);
amrex::Real t_lo = getTgivenPandTh(p0_arr(i,j,k-1), cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1), rd_over_cp);
amrex::Real qminus = (t_lo-t0_lo)/t0_lo;

amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus);
return (-r0_q_avg * grav_gpu);
}

#endif
7 changes: 3 additions & 4 deletions Source/SourceTerms/ERF_make_buoyancy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,6 @@ void make_buoyancy (Vector<MultiFab>& S_data,

// Base state density and pressure
const Array4<const Real>& r0_arr = r0.const_array(mfi);
const Array4<const Real>& p0_arr = p0.const_array(mfi);
const Array4<const Real>& th0_arr = th0.const_array(mfi);

if (solverChoice.moisture_type == MoistureType::None) {
Expand All @@ -81,7 +80,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
//
buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k,
grav_gpu[2],
r0_arr,p0_arr,cell_data);
r0_arr,th0_arr,cell_data);
});
} else {
// NOTE: For decomposition in the vertical direction, klo may not
Expand Down Expand Up @@ -127,7 +126,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,

ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
buoyancy_fab(i, j, k) = buoyancy_type1(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data);
buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data);
});
} // mfi

Expand Down Expand Up @@ -202,7 +201,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,

ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
buoyancy_fab(i, j, k) = buoyancy_type1(i,j,k,n_qstate,
buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_qstate,
grav_gpu[2],r0_arr,cell_data);
});
} // mfi
Expand Down
20 changes: 16 additions & 4 deletions Source/Utils/ERF_AverageDown.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,10 +16,17 @@ void
ERF::AverageDown ()
{
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay);
int src_comp = 0;
int num_comp = vars_new[0][Vars::cons].nComp();

int src_comp, num_comp;
for (int lev = finest_level-1; lev >= 0; --lev)
{
// If anelastic we don't average down rho because rho == rho0.
if (solverChoice.anelastic[lev]) {
src_comp = 1;
} else {
src_comp = 0;
}
num_comp = vars_new[0][Vars::cons].nComp() - src_comp;
AverageDownTo(lev,src_comp,num_comp);
}
}
Expand All @@ -28,8 +35,13 @@ ERF::AverageDown ()
void
ERF::AverageDownTo (int crse_lev, int scomp, int ncomp) // NOLINT
{
AMREX_ALWAYS_ASSERT(scomp == 0);
AMREX_ALWAYS_ASSERT(ncomp == vars_new[crse_lev][Vars::cons].nComp());
if (solverChoice.anelastic[crse_lev]) {
AMREX_ALWAYS_ASSERT(scomp == 1);
} else {
AMREX_ALWAYS_ASSERT(scomp == 0);
}

AMREX_ALWAYS_ASSERT(ncomp == vars_new[crse_lev][Vars::cons].nComp() - scomp);
AMREX_ALWAYS_ASSERT(solverChoice.coupling_type == CouplingType::TwoWay);

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

0 comments on commit d63c6ae

Please sign in to comment.