Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Substep moisture advection #1932

Merged
merged 4 commits into from
Nov 6, 2024
Merged

Substep moisture advection #1932

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
cc6436e
Ran with kessler and sam successfully.
AMLattanzi Nov 6, 2024
82c0ab3
Merge branch 'development' into SubstepMoist
AMLattanzi Nov 6, 2024
48e7860
Merge branch 'development' into SubstepMoist
AMLattanzi Nov 6, 2024
a0aa615
Use numerical limits.
Nov 6, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions Source/Microphysics/Kessler/ERF_Kessler.H
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading