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

Refactor buoyancy and add moist anelastic option. #1914

Merged
merged 7 commits into from
Oct 28, 2024
Merged

Refactor buoyancy and add moist anelastic option. #1914

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
6259b56
Refactor buoyancy and add moist anelastic option.
AMLattanzi Oct 28, 2024
e39c0f6
fix typo
AMLattanzi Oct 28, 2024
d6dedda
fix omp capture and change name.
AMLattanzi Oct 28, 2024
999d116
fix dry default case.
AMLattanzi Oct 28, 2024
c513a1a
Merge branch 'development' into BuoyancyRefactor
AMLattanzi Oct 28, 2024
b04479a
Fix EOS call and make device functions not host.
AMLattanzi Oct 28, 2024
353d500
Forgot one device function.
AMLattanzi Oct 28, 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
253 changes: 253 additions & 0 deletions Source/SourceTerms/ERF_buoyancy_utils.H
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,253 @@
#ifndef ERF_BUOYANCY_UTILS_H_
#define ERF_BUOYANCY_UTILS_H_

#include <ERF_EOS.H>
#include <ERF_Constants.H>

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_dry_anelastic (int& i,
int& j,
int& k,
amrex::Real const& grav_gpu,
amrex::Real const& rd_over_cp,
const amrex::Array4<const amrex::Real>& p0_arr,
const amrex::Array4<const amrex::Real>& r0_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);
}


AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
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>& p0_arr,
const amrex::Array4<const amrex::Real>& r0_arr,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1);
amrex::Real qv_lo = cell_data(i,j,k-1,RhoQ1_comp)/r0_arr(i,j,k-1);
amrex::Real qc_lo = cell_data(i,j,k-1,RhoQ2_comp)/r0_arr(i,j,k-1);
amrex::Real qt_lo = qv_lo + qc_lo;
amrex::Real theta_v_lo = theta_d_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo);

amrex::Real theta_d_hi = cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k);
amrex::Real qv_hi = cell_data(i,j,k,RhoQ1_comp)/r0_arr(i,j,k);
amrex::Real qc_hi = cell_data(i,j,k,RhoQ2_comp)/r0_arr(i,j,k);
amrex::Real qt_hi = qv_hi + qc_hi;
amrex::Real theta_v_hi = theta_d_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi);

amrex::Real theta_v_wface = amrex::Real(0.5) * (theta_v_lo + theta_v_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_v_0_lo = theta_d_0_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo);
amrex::Real theta_v_0_hi = theta_d_0_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi);

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);
}

AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
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>& p0_arr,
const amrex::Array4<const amrex::Real>& r0_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 = 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 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 = 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_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)
{
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);
for (int q_offset(0); q_offset<n_qstate; ++q_offset) {
rhop_hi += cell_data(i,j,k ,RhoQ1_comp+q_offset);
rhop_lo += cell_data(i,j,k-1,RhoQ1_comp+q_offset);
}
return( grav_gpu * amrex::Real(0.5) * ( rhop_hi + rhop_lo ) );
}


AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_type2 (int& i,
int& j,
int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
amrex::Real* rho_d_ptr,
amrex::Real* theta_d_ptr,
amrex::Real* qv_d_ptr,
amrex::Real* qc_d_ptr,
amrex::Real* qp_d_ptr,
const amrex::Array4<const amrex::Real>& cell_prim,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ], qv_d_ptr[k ]);
amrex::Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]);

amrex::Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp),
cell_data(i,j,k ,RhoTheta_comp),
cell_data(i,j,k ,RhoQ1_comp)/cell_data(i,j,k ,Rho_comp));
amrex::Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp),
cell_data(i,j,k-1,RhoTheta_comp),
cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp));

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

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

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

amrex::Real qplus = 0.61 * ( qv_plus - qv_d_ptr[k] ) -
( qc_plus - qc_d_ptr[k] +
qp_plus - qp_d_ptr[k] )
+ (tempp3d-tempp1d)/tempp1d*(amrex::Real(1.0) + amrex::Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qp_d_ptr[k]);

amrex::Real 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*(amrex::Real(1.0) + amrex::Real(0.61)*qv_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]);

amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus);
amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]);

return (-qavg * r0avg * grav_gpu);
}


AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_type3 (int& i,
int& j,
int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
amrex::Real* rho_d_ptr,
amrex::Real* theta_d_ptr,
amrex::Real* qv_d_ptr,
const amrex::Array4<const amrex::Real>& cell_prim,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ], qv_d_ptr[k ]);
amrex::Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]);

amrex::Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp),
cell_data(i,j,k ,RhoTheta_comp),
cell_data(i,j,k ,RhoQ1_comp)/cell_data(i,j,k ,Rho_comp));
amrex::Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp),
cell_data(i,j,k-1,RhoTheta_comp),
cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp));

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

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

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

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

amrex::Real qminus = 0.61 * qv_minus - (qc_minus + qp_minus) + (tempm3d-tempm1d)/tempm1d;

amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus);
amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]);

return ( -qavg * r0avg * grav_gpu );
}


AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
buoyancy_type4 (int& i,
int& j,
int& k,
const int& n_qstate,
amrex::Real const& grav_gpu,
amrex::Real* rho_d_ptr,
amrex::Real* theta_d_ptr,
amrex::Real* qv_d_ptr,
amrex::Real* qc_d_ptr,
amrex::Real* qp_d_ptr,
const amrex::Array4<const amrex::Real>& cell_prim,
const amrex::Array4<const amrex::Real>& cell_data)
{
amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0;
amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0;

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

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

amrex::Real qplus = amrex::Real(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 ];

amrex::Real qminus = amrex::Real(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];

amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus);
amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]);

return (-qavg * r0avg * grav_gpu);
}
#endif
Loading
Loading