Skip to content

Commit

Permalink
Modify buoyancy and pert press. This touches uninitialized data; need…
Browse files Browse the repository at this point in the history 
… to fix.
  • Loading branch information
@AMLattanzi
AMLattanzi committed Dec 19, 2023
1 parent a078cad commit 5bced24
Show file tree
Hide file tree
Showing 11 changed files with 84 additions and 150 deletions.
1 change: 1 addition & 0 deletions Exec/RegTests/Bubble/inputs_squall2d_x
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ erf.use_rayleigh_damping = false # TODO: enable
erf.les_type = "Deardorff"
erf.pbl_type = "None"
erf.moisture_model = "SAM"
erf.buoyancy_type = 1

erf.molec_diff_type = "Constant"
erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities
Expand Down
10 changes: 7 additions & 3 deletions Source/BoundaryConditions/ERF_FillPatch.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -100,9 +100,11 @@ ERF::FillPatch (int lev, Real time, const Vector<MultiFab*>& mfs, bool fillset)
// We call this even if init_type == real because this routine will fill the vertical bcs
(*physbcs[lev])(mfs,icomp_cons,ncomp_cons,ngvect_cons,ngvect_vels,init_type,cons_only,BCVars::cons_bc,time);

/*
// Update vars in the micro model
if (solverChoice.moisture_type != MoistureType::None)
micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
*/
}

/*
Expand Down Expand Up @@ -262,9 +264,11 @@ ERF::FillIntermediatePatch (int lev, Real time,
if (!(cons_only && ncomp_cons == 1) && m_most && allow_most_bcs)
m_most->impose_most_bcs(lev,mfs,eddyDiffs_lev[lev].get(),z_phys_nd[lev].get());

// Update vars in the micro model
if (solverChoice.moisture_type != MoistureType::None)
micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
/*
// Update vars in the micro model
if (solverChoice.moisture_type != MoistureType::None && (!cons_only && ncomp_cons > 2))
micro.Update_Micro_Vars_Lev(lev, *mfs[Vars::cons]);
*/
}

/*
Expand Down
2 changes: 2 additions & 0 deletions Source/DataStructs/DataStruct.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,7 @@ struct SolverChoice {
amrex::Abort("buoyancy_type must be 1, 2, 3 or 4");
}

/*
// Set a different default for moist vs dry
if (moisture_type != MoistureType::None) {
if (moisture_type == MoistureType::FastEddy) {
Expand All @@ -93,6 +94,7 @@ struct SolverChoice {
buoyancy_type = 2; // uses Tprime
}
}
*/

// Is the terrain static or moving?
static std::string terrain_type_string = "Static";
Expand Down
2 changes: 1 addition & 1 deletion Source/ERF_make_new_level.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba,
init_stuff(lev, ba, dm);

int n_qstate = micro.Get_Qstate_Size();
int ncomp_cons = NVAR_max - (3 - n_qstate);
int ncomp_cons = NVAR_max - (NMOIST_max - n_qstate);

lev_new[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state);
lev_old[Vars::cons].define(ba, dm, ncomp_cons, ngrow_state);
Expand Down
2 changes: 1 addition & 1 deletion Source/IO/Plotfile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
Vector<std::string> tmp_plot_names;

int n_qstate = micro.Get_Qstate_Size();
int ncomp_cons = NVAR_max - (3 - n_qstate);
int ncomp_cons = NVAR_max - (NMOIST_max - n_qstate);

for (int i = 0; i < ncomp_cons; ++i) {
if ( containerHasElement(plot_var_names, cons_names[i]) ) {
Expand Down
3 changes: 3 additions & 0 deletions Source/IndexDefines.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,9 @@
// but not to allocate actual solution data
#define NVAR_max 8

// This defines the maximum number of moisture vars
#define NMOIST_max 3

// Cell-centered primitive variables
#define PrimTheta_comp (RhoTheta_comp -1)
#define PrimKE_comp (RhoKE_comp -1)
Expand Down
41 changes: 0 additions & 41 deletions Source/TimeIntegration/ERF_advance_dycore.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,47 +250,6 @@ void ERF::advance_dycore(int level,
cons_to_prim(state_old[IntVar::cons], state_old[IntVar::cons].nGrow());
} // profile

int q_size = micro.Get_Qmoist_Size();
int ncell = fine_geom.Domain().length(2);

Gpu::DeviceVector<Real> qv_d(ncell), qc_d(ncell), qi_d(ncell);

if (solverChoice.moisture_type != MoistureType::None) {
if (q_size >= 1) {
MultiFab qvapor (*(qmoist[level][0]), make_alias, 0, 1);
PlaneAverage qv_ave(&qvapor, geom[level], solverChoice.ave_plane);
qv_ave.compute_averages(ZDir(), qv_ave.field());

Gpu::HostVector <Real> qv_h(ncell);

qv_ave.line_average(0, qv_h);
Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin());
}

if (q_size >= 2) {
MultiFab qcloud (*(qmoist[level][1]), make_alias, 0, 1);
PlaneAverage qc_ave(&qcloud, geom[level], solverChoice.ave_plane);
qc_ave.compute_averages(ZDir(), qc_ave.field());

Gpu::HostVector <Real> qc_h(ncell);

qc_ave.line_average(0, qc_h);
Gpu::copyAsync(Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin());
}

if (q_size >= 3) {
MultiFab qice (*(qmoist[level][2]), make_alias, 0, 1);
PlaneAverage qi_ave(&qice , geom[level], solverChoice.ave_plane);
qi_ave.compute_averages(ZDir(), qi_ave.field());

Gpu::HostVector <Real> qi_h(ncell);

qi_ave.line_average(0, qi_h);
Gpu::copyAsync(Gpu::hostToDevice, qi_h.begin(), qi_h.end(), qi_d.begin());
Gpu::streamSynchronize();
}
}

#include "TI_no_substep_fun.H"
#include "TI_slow_rhs_fun.H"
#include "TI_fast_rhs_fun.H"
Expand Down
155 changes: 68 additions & 87 deletions Source/TimeIntegration/ERF_make_buoyancy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,10 +36,6 @@ using namespace amrex;
void make_buoyancy (Vector<MultiFab>& S_data,
const MultiFab& S_prim,
MultiFab& buoyancy,
Vector<MultiFab*> qmoist,
Gpu::DeviceVector<Real> qv_d,
Gpu::DeviceVector<Real> qc_d,
Gpu::DeviceVector<Real> qi_d,
const amrex::Geometry geom,
const SolverChoice& solverChoice,
const MultiFab* r0)
Expand Down Expand Up @@ -138,35 +134,10 @@ void make_buoyancy (Vector<MultiFab>& S_data,
// ******************************************************************************************
// Moist versions of buoyancy expressions
// ******************************************************************************************
if (solverChoice.moisture_type == MoistureType::FastEddy) {
if (solverChoice.moisture_type != MoistureType::None) {

AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1);

for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
Box tbz = mfi.tilebox();

// We don't compute a source term for z-momentum on the bottom or top domain boundary
if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1);
if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1);

const Array4<const Real> & cell_data = S_data[IntVar::cons].array(mfi);
const Array4< Real> & buoyancy_fab = buoyancy.array(mfi);

// Base state density
const Array4<const Real>& r0_arr = r0->const_array(mfi);

amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Real rhop_hi = cell_data(i,j,k ,Rho_comp) + cell_data(i,j,k ,RhoQ1_comp)
+ cell_data(i,j,k ,RhoQ2_comp) - r0_arr(i,j,k );
Real rhop_lo = cell_data(i,j,k-1,Rho_comp) + cell_data(i,j,k-1,RhoQ1_comp)
+ cell_data(i,j,k-1,RhoQ2_comp) - r0_arr(i,j,k-1);
buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo );
});
} // mfi

} else if (solverChoice.moisture_type != MoistureType::None) {
if (solverChoice.moisture_type == MoistureType::FastEddy)
AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1);

if (solverChoice.buoyancy_type == 1) {

Expand All @@ -184,16 +155,13 @@ void make_buoyancy (Vector<MultiFab>& S_data,
// Base state density
const Array4<const Real>& r0_arr = r0->const_array(mfi);

const Array4<const Real> & qv_data = qmoist[0]->array(mfi);
const Array4<const Real> & qc_data = qmoist[1]->array(mfi);
const Array4<const Real> & qi_data = qmoist[2]->array(mfi);
// TODO: A microphysics model may have more than q1 & q2 components for the
// non-precipitating phase.

amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Real rhop_hi = cell_data(i,j,k ,Rho_comp) * (1.0 + qv_data(i,j,k ) + qc_data(i,j,k )
+ qi_data(i,j,k )) + cell_data(i,j,k,RhoQ2_comp) - r0_arr(i,j,k );
Real rhop_lo = cell_data(i,j,k-1,Rho_comp) * (1.0 + qv_data(i,j,k-1) + qc_data(i,j,k-1)
+ qi_data(i,j,k-1)) + cell_data(i,j,k-1,RhoQ2_comp) - r0_arr(i,j,k-1);
Real rhop_hi = cell_data(i,j,k ,Rho_comp) + cell_data(i,j,k ,RhoQ1_comp) + cell_data(i,j,k ,RhoQ2_comp) - r0_arr(i,j,k );
Real rhop_lo = cell_data(i,j,k-1,Rho_comp) + cell_data(i,j,k-1,RhoQ2_comp) + cell_data(i,j,k-1,RhoQ2_comp) - r0_arr(i,j,k-1);
buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( rhop_hi + rhop_lo );
});
} // mfi
Expand All @@ -209,7 +177,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,

int ncell = state_ave.ncell_line();

Gpu::HostVector <Real> rho_h(ncell), theta_h(ncell), qp_h(ncell);
Gpu::HostVector <Real> rho_h(ncell), theta_h(ncell);
Gpu::DeviceVector<Real> rho_d(ncell), theta_d(ncell);

state_ave.line_average(Rho_comp, rho_h);
Expand All @@ -221,21 +189,28 @@ void make_buoyancy (Vector<MultiFab>& S_data,
Real* rho_d_ptr = rho_d.data();
Real* theta_d_ptr = theta_d.data();

if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) {

prim_ave.line_average(PrimQ2_comp, qp_h);

Gpu::DeviceVector<Real> qp_d(ncell);

// Copy data to device
Gpu::copyAsync(Gpu::hostToDevice, qp_h.begin(), qp_h.end(), qp_d.begin());
Gpu::streamSynchronize();

Real* qp_d_ptr = qp_d.data();
Real* qv_d_ptr = qv_d.data();
Real* qc_d_ptr = qc_d.data();
Real* qi_d_ptr = qi_d.data();
// Average valid moisture vars
int n_prim_max = NVAR_max - 1;
int n_moist_var = NMOIST_max - (S_prim.nComp() - n_prim_max);
Gpu::HostVector <Real> qv_h(ncell) , qc_h(ncell) , qp_h(ncell);
Gpu::DeviceVector<Real> qv_d(ncell,0.0), qc_d(ncell,0.0), qp_d(ncell,0.0);
if (n_moist_var >=1) {
prim_ave.line_average(PrimQ1_comp, qv_h);
Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin());
}
if (n_moist_var >=2) {
prim_ave.line_average(PrimQ2_comp, qc_h);
Gpu::copyAsync(Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin());
}
if (n_moist_var >=3) {
prim_ave.line_average(PrimQ3_comp, qp_h);
Gpu::copyAsync(Gpu::hostToDevice, qp_h.begin(), qp_h.end(), qp_d.begin());
}
Real* qv_d_ptr = qv_d.data();
Real* qc_d_ptr = qc_d.data();
Real* qp_d_ptr = qp_d.data();

if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) {
#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
Expand All @@ -252,10 +227,7 @@ void make_buoyancy (Vector<MultiFab>& S_data,
const Array4<const Real> & cell_data = S_data[IntVar::cons].array(mfi);
const Array4<const Real> & cell_prim = S_prim.array(mfi);

const Array4<const Real> & qv_data = qmoist[0]->const_array(mfi);
const Array4<const Real> & qc_data = qmoist[1]->const_array(mfi);
const Array4<const Real> & qi_data = qmoist[2]->const_array(mfi);

// TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp)
amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ]);
Expand All @@ -266,32 +238,36 @@ void make_buoyancy (Vector<MultiFab>& S_data,

Real qplus, qminus;

Real qv_plus = (n_moist_var >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0;
Real qv_minus = (n_moist_var >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;

Real qc_plus = (n_moist_var >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0;
Real qc_minus = (n_moist_var >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0;

Real qp_plus = (n_moist_var >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0;
Real qp_minus = (n_moist_var >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0;

if (solverChoice.buoyancy_type == 2) {
qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) -
(qc_data(i,j,k)-qc_d_ptr[k]+
qi_data(i,j,k)-qi_d_ptr[k]+
cell_prim(i,j,k,PrimQ2_comp)-qp_d_ptr[k])
+ (tempp3d-tempp1d)/tempp1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qi_d_ptr[k]-qp_d_ptr[k]);

qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) -
(qc_data(i,j,k-1)-qc_d_ptr[k-1]+
qi_data(i,j,k-1)-qi_d_ptr[k-1]+
cell_prim(i,j,k-1,PrimQ2_comp)-qp_d_ptr[k-1])
+ (tempm3d-tempm1d)/tempm1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k-1]-qi_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]);
qplus = 0.61 * ( qv_plus - qv_d_ptr[k] ) -
( qc_plus - qc_d_ptr[k] +
qp_plus - qp_d_ptr[k] )
+ (tempp3d-tempp1d)/tempp1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qp_d_ptr[k]);

} else if (solverChoice.buoyancy_type == 4) {
qminus = 0.61 * ( qv_minus - qv_d_ptr[k-1] ) -
( qc_minus - qc_d_ptr[k-1] +
qp_minus - qp_d_ptr[k-1] )
+ (tempm3d-tempm1d)/tempm1d*(Real(1.0) + Real(0.61)*qv_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]);

qplus = 0.61* ( qv_data(i,j,k)-qv_d_ptr[k]) -
(qc_data(i,j,k)-qc_d_ptr[k]+
qi_data(i,j,k)-qi_d_ptr[k]+
cell_prim(i,j,k,PrimQ2_comp)-qp_d_ptr[k])
} else if (solverChoice.buoyancy_type == 4) {
qplus = 0.61 * ( qv_plus - qv_d_ptr[k] ) -
( qc_plus - qc_d_ptr[k] +
qp_plus - qp_d_ptr[k] )
+ (cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - theta_d_ptr[k ])/theta_d_ptr[k ];

qminus = 0.61 *( qv_data(i,j,k-1)-qv_d_ptr[k-1]) -
(qc_data(i,j,k-1)-qc_d_ptr[k-1]+
qi_data(i,j,k-1)-qi_d_ptr[k-1]+
cell_prim(i,j,k-1,PrimQ2_comp)-qp_d_ptr[k-1])
+ (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1];
qminus = 0.61 * ( qv_minus - qv_d_ptr[k-1] ) -
( qc_minus - qc_d_ptr[k-1] +
qp_minus - qp_d_ptr[k-1] )
+ (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1];
}

Real qavg = Real(0.5) * (qplus + qminus);
Expand All @@ -302,7 +278,6 @@ void make_buoyancy (Vector<MultiFab>& S_data,
} // mfi

} else if (solverChoice.buoyancy_type == 3) {

#ifdef _OPENMP
#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
#endif
Expand All @@ -318,9 +293,8 @@ void make_buoyancy (Vector<MultiFab>& S_data,

const Array4<const Real> & cell_data = S_data[IntVar::cons].array(mfi);
const Array4<const Real> & cell_prim = S_prim.array(mfi);
const Array4<const Real> & qv_data = qmoist[0]->const_array(mfi);
const Array4<const Real> & qc_data = qmoist[1]->const_array(mfi);
const Array4<const Real> & qi_data = qmoist[2]->const_array(mfi);

// TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp)

amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k)
{
Expand All @@ -330,11 +304,18 @@ void make_buoyancy (Vector<MultiFab>& S_data,
Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp));
Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp));

Real qplus = 0.61 * qv_data(i,j,k) - (qc_data(i,j,k)+ qi_data(i,j,k)+ cell_prim(i,j,k,PrimQ2_comp))
+ (tempp3d-tempp1d)/tempp1d;
Real qv_plus = (n_moist_var >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0;
Real qv_minus = (n_moist_var >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;

Real qc_plus = (n_moist_var >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0;
Real qc_minus = (n_moist_var >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0;

Real qp_plus = (n_moist_var >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0;
Real qp_minus = (n_moist_var >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0;

Real qplus = 0.61 * qv_plus - (qc_plus + qp_plus) + (tempp3d-tempp1d)/tempp1d;

Real qminus = 0.61 *qv_data(i,j,k-1) - (qc_data(i,j,k-1)+ qi_data(i,j,k-1)+ cell_prim(i,j,k-1,PrimQ2_comp))
+ (tempm3d-tempm1d)/tempm1d;
Real qminus = 0.61 * qv_minus - (qc_minus + qp_minus) + (tempm3d-tempm1d)/tempm1d;

Real qavg = Real(0.5) * (qplus + qminus);
Real r0avg = Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]);
Expand Down
7 changes: 1 addition & 6 deletions Source/TimeIntegration/ERF_slow_rhs_pre.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,6 @@ void erf_slow_rhs_pre (int level, int finest_level,
const MultiFab& xvel,
const MultiFab& yvel,
const MultiFab& zvel,
Vector<MultiFab*> qmoist,
std::unique_ptr<MultiFab>& z_t_mf,
MultiFab& Omega,
const MultiFab& source,
Expand Down Expand Up @@ -543,18 +542,14 @@ void erf_slow_rhs_pre (int level, int finest_level,
FArrayBox pprime; pprime.resize(gbx,1);
Elixir pp_eli = pprime.elixir();
const Array4<Real > & pp_arr = pprime.array();
const Array4<Real const> & qv_arr = (use_moisture) ? qmoist[0]->const_array(mfi) : Array4<Real> {};
{
BL_PROFILE("slow_rhs_pre_pprime");
ParallelFor(gbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
//if (cell_data(i,j,k,RhoTheta_comp) < 0.) printf("BAD THETA AT %d %d %d %e %e \n",
// i,j,k,cell_data(i,j,k,RhoTheta_comp),cell_data(i,j,k+1,RhoTheta_comp));
AMREX_ASSERT(cell_data(i,j,k,RhoTheta_comp) > 0.);
Real qv_for_p = 0.0;
if (use_moisture) {
qv_for_p = qv_arr(i,j,k);
}
Real qv_for_p = (use_moisture) ? cell_data(i,j,k,RhoTheta_comp) : 0.0;
pp_arr(i,j,k) = getPgivenRTh(cell_data(i,j,k,RhoTheta_comp),qv_for_p) - p0_arr(i,j,k);
});
} // end profile
Expand Down
Loading

0 comments on commit 5bced24

Please sign in to comment.