Skip to content

Commit

Permalink
Substep moisture advection (erf-model#1932)
Browse files Browse the repository at this point in the history 
* Ran with kessler and sam successfully.

* Use numerical limits.

---------

Co-authored-by: Aaron Lattanzi <[email protected]>
  • Loading branch information
@AMLattanzi
AMLattanzi and Aaron Lattanzi authored Nov 6, 2024
1 parent 18c66d1 commit 1b67c57
Show file tree
Hide file tree
Showing 5 changed files with 336 additions and 151 deletions.
3 changes: 3 additions & 0 deletions Source/Microphysics/Kessler/ERF_Kessler.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,6 +143,9 @@ private:
// Number of qstate variables
int m_qstate_size = 3;

// CFL MAX for vertical advection
static constexpr amrex::Real CFL_MAX = 0.5;

// MicVar map (Qmoist indices -> MicVar enum)
amrex::Vector<int> MicVarMap;

Expand Down
185 changes: 124 additions & 61 deletions Source/Microphysics/Kessler/ERF_Kessler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,53 +22,10 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
auto dm = tabs->DistributionMap();
fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells

Real dtn = dt;

for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){
auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);
auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);
auto rain_accum_array = mic_fab_vars[MicVar_Kess::rain_accum]->array(mfi);

auto fz_array = fz.array(mfi);
const Box& tbz = mfi.tilebox();

ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real rho_avg, qp_avg;

if (k==k_lo) {
rho_avg = rho_array(i,j,k);
qp_avg = qp_array(i,j,k);
} else if (k==k_hi+1) {
rho_avg = rho_array(i,j,k-1);
qp_avg = qp_array(i,j,k-1);
} else {
rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3
qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
}

qp_avg = std::max(0.0, qp_avg);

Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s

// NOTE: Fz is the sedimentation flux from the advective operator.
// In the terrain-following coordinate system, the z-deriv in
// the divergence uses the normal velocity (Omega). However,
// there are no u/v components to the sedimentation velocity.
// Therefore, we simply end up with a division by detJ when
// evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;

if(k==k_lo){
rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*qp_avg*V_terminal*dtn/1000.0*1000.0; // Divide by rho_water and convert to mm
}

/*if(k==0){
fz_array(i,j,k) = 0;
}*/
});
}
Real dtn = dt;
Real coef = dtn/dz;

// Saturation and evaporation calculations go first
for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi);
auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi);
Expand All @@ -79,20 +36,13 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi);
auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);

const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4<const Real>{};

const auto& box3d = mfi.tilebox();

auto fz_array = fz.array(mfi);
const auto& tbx = mfi.tilebox();

// Expose for GPU
Real d_fac_cond = m_fac_cond;

ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Jacobian determinant
Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0;

qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k));
qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k));
qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
Expand Down Expand Up @@ -163,14 +113,9 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
dq_clwater_to_rain = std::min(dq_clwater_to_rain, qc_array(i,j,k));
}

if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0;
if(std::fabs(fz_array(i,j,k )) < 1e-14) fz_array(i,j,k ) = 0.0;
Real dq_sed = dtn * dJinv * (1.0/rho_array(i,j,k)) * (fz_array(i,j,k+1) - fz_array(i,j,k))/dz;
if(std::fabs(dq_sed) < 1e-14) dq_sed = 0.0;

qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor + dq_rain_to_vapor;
qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;
qp_array(i,j,k) += dq_sed + dq_clwater_to_rain - dq_rain_to_vapor;
qp_array(i,j,k) += dq_clwater_to_rain - dq_rain_to_vapor;

Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k);
theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor);
Expand All @@ -182,8 +127,126 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice)
qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k);
});
}

// Precompute terminal velocity for substepping
for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){
auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);
auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);

auto fz_array = fz.array(mfi);
const Box& tbz = mfi.tilebox();

ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real rho_avg, qp_avg;

if (k==k_lo) {
rho_avg = rho_array(i,j,k);
qp_avg = qp_array(i,j,k);
} else if (k==k_hi+1) {
rho_avg = rho_array(i,j,k-1);
qp_avg = qp_array(i,j,k-1);
} else {
rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3
qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
}

qp_avg = std::max(0.0, qp_avg);

Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s

// NOTE: Fz is the sedimentation flux from the advective operator.
// In the terrain-following coordinate system, the z-deriv in
// the divergence uses the normal velocity (Omega). However,
// there are no u/v components to the sedimentation velocity.
// Therefore, we simply end up with a division by detJ when
// evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;
});
} // mfi

// Compute number of substeps from maximum terminal velocity
Real wt_max;
int n_substep;
auto const& ma_fz_arr = fz.const_arrays();
GpuTuple<Real> max = ParReduce(TypeList<ReduceOpMax>{},
TypeList<Real>{},
fz, IntVect(0),
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-> GpuTuple<Real>
{
return { ma_fz_arr[box_no](i,j,k) };
});
wt_max = get<0>(max) + std::numeric_limits<Real>::epsilon();
n_substep = int( std::ceil(wt_max * coef / CFL_MAX) );
AMREX_ALWAYS_ASSERT(n_substep >= 1);
coef /= Real(n_substep);
dtn /= Real(n_substep);

// Substep the vertical advection
for (int nsub(0); nsub<n_substep; ++nsub) {
for ( MFIter mfi(*tabs, TilingIfNotGPU()); mfi.isValid(); ++mfi ){
auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi);
auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi);
auto rain_accum_array = mic_fab_vars[MicVar_Kess::rain_accum]->array(mfi);
auto fz_array = fz.array(mfi);

const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4<const Real>{};

const Box& tbx = mfi.tilebox();
const Box& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0));

// Update vertical flux every substep
ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real rho_avg, qp_avg;
if (k==k_lo) {
rho_avg = rho_array(i,j,k);
qp_avg = qp_array(i,j,k);
} else if (k==k_hi+1) {
rho_avg = rho_array(i,j,k-1);
qp_avg = qp_array(i,j,k-1);
} else {
rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3
qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k));
}

qp_avg = std::max(0.0, qp_avg);

Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s

// NOTE: Fz is the sedimentation flux from the advective operator.
// In the terrain-following coordinate system, the z-deriv in
// the divergence uses the normal velocity (Omega). However,
// there are no u/v components to the sedimentation velocity.
// Therefore, we simply end up with a division by detJ when
// evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
fz_array(i,j,k) = rho_avg*V_terminal*qp_avg;

if(k==k_lo){
rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*qp_avg*V_terminal*dtn/1000.0*1000.0; // Divide by rho_water and convert to mm
}
});

// Update precip every substep
ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Jacobian determinant
Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0;

if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0;
if(std::fabs(fz_array(i,j,k )) < 1e-14) fz_array(i,j,k ) = 0.0;
Real dq_sed = dJinv * (1.0/rho_array(i,j,k)) * (fz_array(i,j,k+1) - fz_array(i,j,k)) * coef;
if(std::fabs(dq_sed) < 1e-14) dq_sed = 0.0;

qp_array(i,j,k) += dq_sed;
qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k));
});
} // mfi
} // nsub
}


if (solverChoice.moisture_type == MoistureType::Kessler_NoRain){

// get the temperature, dentisy, theta, qt and qc from input
Expand Down
115 changes: 82 additions & 33 deletions Source/Microphysics/SAM/ERF_IceFall.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#include <AMReX_ParReduce.H>
#include "ERF_SAM.H"
#include "ERF_TileNoZ.H"
#include <cmath>

using namespace amrex;

Expand Down Expand Up @@ -55,7 +56,7 @@ void SAM::IceFall (const SolverChoice& sc) {
rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k));
qci_avg = 0.5*(qci_array(i,j,k-1) + qci_array(i,j,k));
}
Real vt_ice = min( 0.4 , 8.66 * pow( (max(0.,qci_avg)+1.e-10) , 0.24) );
Real vt_ice = min( 0.4 , 8.66 * pow( (std::max(0.,qci_avg)+1.e-10) , 0.24) );

// NOTE: Fz is the sedimentation flux from the advective operator.
// In the terrain-following coordinate system, the z-deriv in
Expand All @@ -67,37 +68,85 @@ void SAM::IceFall (const SolverChoice& sc) {
});
}

for (MFIter mfi(*qci, TileNoZ()); mfi.isValid(); ++mfi) {
auto qci_array = qci->array(mfi);
auto qn_array = qn->array(mfi);
auto qt_array = qt->array(mfi);
auto rho_array = rho->array(mfi);
auto fz_array = fz.array(mfi);

const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4<const Real>{};

const auto& box3d = mfi.tilebox();

ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Jacobian determinant
Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0;

//==================================================
// Cloud ice sedimentation (A32)
//==================================================
Real dqi = dJinv * (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef;
dqi = std::max(-qci_array(i,j,k), dqi);

// Add this increment to both non-precipitating and total water.
qci_array(i,j,k) += dqi;
qn_array(i,j,k) += dqi;
qt_array(i,j,k) += dqi;

// NOTE: Sedimentation does not affect the potential temperature,
// but it does affect the liquid/ice static energy.
// No source to Theta occurs here.
});
}
// Compute number of substeps from maximum terminal velocity
Real wt_max;
int n_substep;
auto const& ma_fz_arr = fz.const_arrays();
GpuTuple<Real> max = ParReduce(TypeList<ReduceOpMax>{},
TypeList<Real>{},
fz, IntVect(0),
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-> GpuTuple<Real>
{
return { ma_fz_arr[box_no](i,j,k) };
});
wt_max = get<0>(max) + std::numeric_limits<Real>::epsilon();
n_substep = int( std::ceil(wt_max * coef / CFL_MAX) );
AMREX_ALWAYS_ASSERT(n_substep >= 1);
coef /= Real(n_substep);
dtn /= Real(n_substep);

// Substep the vertical advection
for (int nsub(0); nsub<n_substep; ++nsub) {
for (MFIter mfi(*qci, TileNoZ()); mfi.isValid(); ++mfi) {
auto qci_array = qci->array(mfi);
auto qn_array = qn->array(mfi);
auto qt_array = qt->array(mfi);
auto rho_array = rho->array(mfi);
auto fz_array = fz.array(mfi);

const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4<const Real>{};

const auto& tbx = mfi.tilebox();
const auto& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0));

// Update vertical flux every substep
ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
Real rho_avg, qci_avg;
if (k==k_lo) {
rho_avg = rho_array(i,j,k);
qci_avg = qci_array(i,j,k);
} else if (k==k_hi+1) {
rho_avg = rho_array(i,j,k-1);
qci_avg = qci_array(i,j,k-1);
} else {
rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k));
qci_avg = 0.5*(qci_array(i,j,k-1) + qci_array(i,j,k));
}
Real vt_ice = min( 0.4 , 8.66 * pow( (std::max(0.,qci_avg)+1.e-10) , 0.24) );

// NOTE: Fz is the sedimentation flux from the advective operator.
// In the terrain-following coordinate system, the z-deriv in
// the divergence uses the normal velocity (Omega). However,
// there are no u/v components to the sedimentation velocity.
// Therefore, we simply end up with a division by detJ when
// evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
fz_array(i,j,k) = rho_avg*vt_ice*qci_avg;
});

// Update precip every substep
ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
{
// Jacobian determinant
Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0;

//==================================================
// Cloud ice sedimentation (A32)
//==================================================
Real dqi = dJinv * (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef;
dqi = std::max(-qci_array(i,j,k), dqi);

// Add this increment to both non-precipitating and total water.
qci_array(i,j,k) += dqi;
qn_array(i,j,k) += dqi;
qt_array(i,j,k) += dqi;

// NOTE: Sedimentation does not affect the potential temperature,
// but it does affect the liquid/ice static energy.
// No source to Theta occurs here.
});
} // mfi
} // nsub
}

Loading

0 comments on commit 1b67c57

Please sign in to comment.