Skip to content

Commit

Permalink
BOMEX Subsidence mod (#1674)
Browse files Browse the repository at this point in the history 
* Match FE source code.

* Use FE qsat and balance with moisture.

* Move wbar to z-faces and average.

* Correct z location for wbar on z-face.

* Revert qsat.
  • Loading branch information
@AMLattanzi
AMLattanzi authored Jul 9, 2024
1 parent 60c1245 commit bc278b3
Show file tree
Hide file tree
Showing 5 changed files with 185 additions and 69 deletions.
8 changes: 4 additions & 4 deletions Exec/DevTests/Bomex/prob.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -325,7 +325,7 @@ Problem::update_w_subsidence (const Real& /*time*/,
{
if (wbar.empty()) return;

const int khi = geom.Domain().bigEnd()[2];
const int khi = geom.Domain().bigEnd()[2] + 1; // lives on z-faces
const Real* prob_lo = geom.ProbLo();
const auto dx = geom.CellSize();

Expand All @@ -343,9 +343,9 @@ Problem::update_w_subsidence (const Real& /*time*/,
Real z_0 = 0.0; // (z_phys_cc) ? zlevels[0] : prob_lo[2] + 0.5 * dx[2];
Real slope1 = parms.wbar_sub_max / (parms.wbar_cutoff_max - z_0);
Real slope2 = -parms.wbar_sub_max / (parms.wbar_cutoff_min - parms.wbar_cutoff_max);
// wbar[0] = 0.0;
for (int k = 0; k <= khi; k++) {
const Real z_cc = (z_phys_cc) ? zlevels[k] : prob_lo[2] + (k+0.5)* dx[2];
wbar[0] = 0.0;
for (int k = 1; k <= khi; k++) {
const Real z_cc = (z_phys_cc) ? zlevels[k] : prob_lo[2] + k*dx[2];
if (z_cc <= parms.wbar_cutoff_max) {
wbar[k] = slope1 * (z_cc - z_0);
} else if (z_cc <= parms.wbar_cutoff_min) {
Expand Down
2 changes: 1 addition & 1 deletion Source/ERF.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -774,7 +774,7 @@ ERF::InitData ()
h_w_subsid.resize(max_level+1, Vector<Real>(0));
d_w_subsid.resize(max_level+1, Gpu::DeviceVector<Real>(0));
for (int lev = 0; lev <= finest_level; lev++) {
const int domlen = geom[lev].Domain().length(2);
const int domlen = geom[lev].Domain().length(2) + 1; // lives on z-faces
h_w_subsid[lev].resize(domlen, 0.0_rt);
d_w_subsid[lev].resize(domlen, 0.0_rt);
prob->update_w_subsidence(t_new[0],
Expand Down
90 changes: 72 additions & 18 deletions Source/SourceTerms/ERF_make_mom_sources.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,12 +110,38 @@ void make_mom_sources (int /*level*/,
// *****************************************************************************
// Planar averages for subsidence terms
// *****************************************************************************
Table1D<Real> dptr_u_plane, dptr_v_plane;
TableData<Real, 1> u_plane_tab, v_plane_tab;
Table1D<Real> dptr_r_plane, dptr_u_plane, dptr_v_plane;
TableData<Real, 1> r_plane_tab, u_plane_tab, v_plane_tab;

if (dptr_wbar_sub) {
PlaneAverage u_ave(&xvel , geom, solverChoice.ave_plane, true);
PlaneAverage v_ave(&yvel , geom, solverChoice.ave_plane, true);
// Rho
PlaneAverage r_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane, true);
r_ave.compute_averages(ZDir(), r_ave.field());

int ncell = r_ave.ncell_line();
Gpu::HostVector< Real> r_plane_h(ncell);
Gpu::DeviceVector< Real> r_plane_d(ncell);

r_ave.line_average(Rho_comp, r_plane_h);

Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());

Real* dptr_r = r_plane_d.data();

IntVect ng_c = S_data[IntVars::cons].nGrowVect();
Box tdomain = domain; tdomain.grow(2,ng_c[2]);
r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});

int offset = ng_c[2];
dptr_r_plane = r_plane_tab.table();
ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
{
dptr_r_plane(k-offset) = dptr_r[k];
});

// U and V momentum
PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, true);
PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, true);

u_ave.compute_averages(ZDir(), u_ave.field());
v_ave.compute_averages(ZDir(), v_ave.field());
Expand All @@ -125,17 +151,17 @@ void make_mom_sources (int /*level*/,
Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell);
Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell);

u_ave.line_average(0 , u_plane_h);
v_ave.line_average(0 , v_plane_h);
u_ave.line_average(0 , u_plane_h);
v_ave.line_average(0 , v_plane_h);

Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin());
Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin());

Real* dptr_u = u_plane_d.data();
Real* dptr_v = v_plane_d.data();

IntVect ng_u = xvel.nGrowVect();
IntVect ng_v = yvel.nGrowVect();
IntVect ng_u = S_data[IntVars::xmom].nGrowVect();
IntVect ng_v = S_data[IntVars::ymom].nGrowVect();
Box udomain = domain; udomain.grow(2,ng_u[2]);
Box vdomain = domain; vdomain.grow(2,ng_v[2]);
u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)});
Expand Down Expand Up @@ -280,16 +306,44 @@ void make_mom_sources (int /*level*/,
// Add custom SUBSIDENCE terms
// *****************************************************************************
if (solverChoice.custom_w_subsidence) {
ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
xmom_src_arr(i, j, k) -= dptr_wbar_sub[k] *
0.5 * (dptr_u_plane(k+1) - dptr_u_plane(k-1)) * dxInv[2];
});
ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
ymom_src_arr(i, j, k) -= dptr_wbar_sub[k] *
0.5 * (dptr_v_plane(k+1) - dptr_v_plane(k-1)) * dxInv[2];
});
if (solverChoice.custom_forcing_prim_vars) {
const int nr = Rho_comp;
ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf *
0.5 * (U_hi - U_lo) * dxInv[2];
});
ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf *
0.5 * (V_hi - V_lo) * dxInv[2];
});
} else {
ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
xmom_src_arr(i, j, k) -= wbar_xf *
0.5 * (U_hi - U_lo) * dxInv[2];
});
ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
ymom_src_arr(i, j, k) -= wbar_yf *
0.5 * (V_hi - V_lo) * dxInv[2];
});
}
}

// *****************************************************************************
Expand Down
118 changes: 72 additions & 46 deletions Source/SourceTerms/ERF_make_sources.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,28 +70,50 @@ void make_sources (int level,
// *****************************************************************************
// Planar averages for subsidence terms
// *****************************************************************************
Table1D<Real> dptr_t_plane, dptr_qv_plane, dptr_qc_plane;
TableData<Real, 1> t_plane_tab, qv_plane_tab, qc_plane_tab;
Table1D<Real> dptr_r_plane, dptr_t_plane, dptr_qv_plane, dptr_qc_plane;
TableData<Real, 1> r_plane_tab, t_plane_tab, qv_plane_tab, qc_plane_tab;
if (dptr_wbar_sub)
{
PlaneAverage t_ave(&S_prim, geom, solverChoice.ave_plane, true);
// Rho
PlaneAverage r_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane, true);
r_ave.compute_averages(ZDir(), r_ave.field());

int ncell = r_ave.ncell_line();
Gpu::HostVector< Real> r_plane_h(ncell);
Gpu::DeviceVector< Real> r_plane_d(ncell);

r_ave.line_average(Rho_comp, r_plane_h);

Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin());

Real* dptr_r = r_plane_d.data();

IntVect ng_c = S_data[IntVars::cons].nGrowVect();
Box tdomain = domain; tdomain.grow(2,ng_c[2]);
r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});

int offset = ng_c[2];
dptr_r_plane = r_plane_tab.table();
ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
{
dptr_r_plane(k-offset) = dptr_r[k];
});

// Rho * Theta
PlaneAverage t_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane, true);
t_ave.compute_averages(ZDir(), t_ave.field());

int ncell = t_ave.ncell_line();
Gpu::HostVector< Real> t_plane_h(ncell);
Gpu::DeviceVector< Real> t_plane_d(ncell);

t_ave.line_average(PrimTheta_comp, t_plane_h);
t_ave.line_average(RhoTheta_comp, t_plane_h);

Gpu::copyAsync(Gpu::hostToDevice, t_plane_h.begin(), t_plane_h.end(), t_plane_d.begin());

Real* dptr_t = t_plane_d.data();

IntVect ng_c = S_prim.nGrowVect();
Box tdomain = domain; tdomain.grow(2,ng_c[2]);
t_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});

int offset = ng_c[2];
dptr_t_plane = t_plane_tab.table();
ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
{
Expand All @@ -101,13 +123,13 @@ void make_sources (int level,
if (solverChoice.moisture_type != MoistureType::None)
{
// Water vapor
PlaneAverage qv_ave(&S_prim, geom, solverChoice.ave_plane, true);
PlaneAverage qv_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane, true);
qv_ave.compute_averages(ZDir(), qv_ave.field());

Gpu::HostVector< Real> qv_plane_h(ncell);
Gpu::DeviceVector<Real> qv_plane_d(ncell);

qv_ave.line_average(PrimQ1_comp, qv_plane_h);
qv_ave.line_average(RhoQ1_comp, qv_plane_h);
Gpu::copyAsync(Gpu::hostToDevice, qv_plane_h.begin(), qv_plane_h.end(), qv_plane_d.begin());

Real* dptr_qv = qv_plane_d.data();
Expand All @@ -119,26 +141,6 @@ void make_sources (int level,
{
dptr_qv_plane(k-offset) = dptr_qv[k];
});

// Cloud water
PlaneAverage qc_ave(&S_prim, geom, solverChoice.ave_plane, true);
qc_ave.compute_averages(ZDir(), qc_ave.field());

Gpu::HostVector< Real> qc_plane_h(ncell);
Gpu::DeviceVector<Real> qc_plane_d(ncell);

qc_ave.line_average(PrimQ2_comp, qc_plane_h);
Gpu::copyAsync(Gpu::hostToDevice, qc_plane_h.begin(), qc_plane_h.end(), qc_plane_d.begin());

Real* dptr_qc = qc_plane_d.data();

qc_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)});

dptr_qc_plane = qc_plane_tab.table();
ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept
{
dptr_qc_plane(k-offset) = dptr_qc[k];
});
}
}

Expand Down Expand Up @@ -239,12 +241,26 @@ void make_sources (int level,
// *************************************************************************************
if (solverChoice.custom_w_subsidence) {
const int n = RhoTheta_comp;
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
cell_src(i, j, k, n) -= dptr_wbar_sub[k] *
0.5 * (dptr_t_plane(k+1) - dptr_t_plane(k-1)) * dxInv[2];

});
if (solverChoice.custom_forcing_prim_vars) {
const int nr = Rho_comp;
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real T_hi = dptr_t_plane(k+1) / dptr_r_plane(k+1);
Real T_lo = dptr_t_plane(k-1) / dptr_r_plane(k-1);
Real wbar_cc = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
cell_src(i, j, k, n) -= cell_data(i,j,k,nr) * wbar_cc *
0.5 * (T_hi - T_lo) * dxInv[2];
});
} else {
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real T_hi = dptr_t_plane(k+1) / dptr_r_plane(k+1);
Real T_lo = dptr_t_plane(k-1) / dptr_r_plane(k-1);
Real wbar_cc = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
cell_src(i, j, k, n) -= wbar_cc *
0.5 * (T_hi - T_lo) * dxInv[2];
});
}
}

// *************************************************************************************
Expand All @@ -253,16 +269,26 @@ void make_sources (int level,
if (solverChoice.custom_w_subsidence && (solverChoice.moisture_type != MoistureType::None)) {
const int nv = RhoQ1_comp;
const int nc = RhoQ2_comp;
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
cell_src(i,j,k,nv) -= dptr_wbar_sub[k] *
0.5 * (dptr_qv_plane(k+1) - dptr_qv_plane(k-1)) * dxInv[2];
});
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
cell_src(i,j,k,nc) -= dptr_wbar_sub[k] *
0.5 * (dptr_qc_plane(k+1) - dptr_qc_plane(k-1)) * dxInv[2];
});
if (solverChoice.custom_forcing_prim_vars) {
const int nr = Rho_comp;
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real Qv_hi = dptr_qv_plane(k+1) / dptr_r_plane(k+1);
Real Qv_lo = dptr_qv_plane(k-1) / dptr_r_plane(k-1);
Real wbar_cc = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
cell_src(i, j, k, nv) -= cell_data(i,j,k,nr) * wbar_cc *
0.5 * (Qv_hi - Qv_lo) * dxInv[2];
});
} else {
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real Qv_hi = dptr_qv_plane(k+1) / dptr_r_plane(k+1);
Real Qv_lo = dptr_qv_plane(k-1) / dptr_r_plane(k-1);
Real wbar_cc = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
cell_src(i, j, k, nv) -= wbar_cc *
0.5 * (Qv_hi - Qv_lo) * dxInv[2];
});
}
}

// *************************************************************************************
Expand Down
36 changes: 36 additions & 0 deletions Source/Utils/Microphysics_Utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ amrex::Real erf_esati (amrex::Real t) {

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real erf_esatw (amrex::Real t) {
#if 1
amrex::Real const a0 = 6.105851;
amrex::Real const a1 = 0.4440316;
amrex::Real const a2 = 0.1430341e-1;
Expand All @@ -62,6 +63,37 @@ amrex::Real erf_esatw (amrex::Real t) {
esatw = 2.0*0.01*std::exp(9.550426 - 5723.265/t + 3.53068*std::log(t) - 0.00728332*t);
}
return esatw;
#else
amrex::Real a0 = 6.11239921;
amrex::Real a1 = 0.443987641;
amrex::Real a2 = 0.142986287e-1;
amrex::Real a3 = 0.264847430e-3;
amrex::Real a4 = 0.302950461e-5;
amrex::Real a5 = 0.206739458e-7;
amrex::Real a6 = 0.640689451e-10;
amrex::Real a7 = -0.952447341e-13;
amrex::Real a8 = -0.976195544e-15;
amrex::Real a0i = 6.11147274;
amrex::Real a1i = 0.503160820;
amrex::Real a2i = 0.188439774e-1;
amrex::Real a3i = 0.420895665e-3;
amrex::Real a4i = 0.615021634e-5;
amrex::Real a5i = 0.602588177e-7;
amrex::Real a6i = 0.385852041e-9;
amrex::Real a7i = 0.146898966e-11;
amrex::Real a8i = 0.252751365e-14;
amrex::Real Tc = std::max(-80.0, t-273.16);
amrex::Real esatw;
if (t>=273.15){
esatw = (a0 + Tc*(a1 +Tc*(a2 +Tc*(a3 +Tc*(a4 +Tc*(a5 +Tc*(a6 +Tc*(a7 +a8 *Tc))))))))*100.;
}else{
esatw = (a0i + Tc*(a1i+Tc*(a2i+Tc*(a3i+Tc*(a4i+Tc*(a5i+Tc*(a6i+Tc*(a7i+a8i*Tc))))))))*100.;
}
return esatw;
#endif
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand Down Expand Up @@ -122,7 +154,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void erf_qsatw (amrex::Real t, amrex::Real p, amrex::Real &qsatw) {
amrex::Real esatw;
esatw = erf_esatw(t);
#if 1
qsatw = Rd_on_Rv*esatw/std::max(esatw,p-esatw);
#else
qsatw = esatw/(R_v*t);
#endif
}

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand Down

0 comments on commit bc278b3

Please sign in to comment.