diff --git a/Source/SourceTerms/ERF_buoyancy_utils.H b/Source/SourceTerms/ERF_buoyancy_utils.H new file mode 100644 index 000000000..91f2a8ea0 --- /dev/null +++ b/Source/SourceTerms/ERF_buoyancy_utils.H @@ -0,0 +1,253 @@ +#ifndef ERF_BUOYANCY_UTILS_H_ +#define ERF_BUOYANCY_UTILS_H_ + +#include +#include + +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& p0_arr, + const amrex::Array4& r0_arr, + const amrex::Array4& 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& p0_arr, + const amrex::Array4& r0_arr, + const amrex::Array4& 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& p0_arr, + const amrex::Array4& r0_arr, + const amrex::Array4& 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& r0_arr, + const amrex::Array4& 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& cell_prim, + const amrex::Array4& 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& cell_prim, + const amrex::Array4& 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& cell_prim, + const amrex::Array4& 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 diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp index d9ec30a53..11a3743a9 100644 --- a/Source/SourceTerms/ERF_make_buoyancy.cpp +++ b/Source/SourceTerms/ERF_make_buoyancy.cpp @@ -7,7 +7,7 @@ #include #include #include -//#include +#include using namespace amrex; @@ -48,6 +48,7 @@ void make_buoyancy (Vector& S_data, const int khi = geom.Domain().bigEnd()[2] + 1; Real rd_over_cp = solverChoice.rdOcp; + Real rv_over_rd = R_v/R_d; if (anelastic == 1) { #ifdef _OPENMP @@ -68,21 +69,24 @@ void make_buoyancy (Vector& S_data, const Array4& r0_arr = r0->const_array(mfi); const Array4& p0_arr = p0->const_array(mfi); - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k)); - Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp); - Real t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp); - Real qplus = (t_hi-t0_hi)/t0_hi; - - Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)); - Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp); - 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); - Real qminus = (t_lo-t0_lo)/t0_lo; - - Real r0_q_avg = Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus); - buoyancy_fab(i, j, k) = -r0_q_avg * grav_gpu[2]; - }); + if (solverChoice.moisture_type == MoistureType::None) { + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k, + grav_gpu[2],rd_over_cp, + p0_arr,r0_arr,cell_data); + }); + } else { + // NOTE: For decomposition in the vertical direction, klo may not + // reside in the valid box and this call will yield an out + // of bounds error since it depends upon the surface theta_l + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k, + grav_gpu[2],rv_over_rd, + p0_arr,r0_arr,cell_data); + }); + } } // mfi } else @@ -92,75 +96,66 @@ void make_buoyancy (Vector& S_data, // ****************************************************************************************** if (solverChoice.moisture_type == MoistureType::None) { - if (solverChoice.buoyancy_type == 1) { + int n_q_dry = 0; + if (solverChoice.buoyancy_type == 1) { #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); + 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); + // 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 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); + const Array4 & cell_data = S_data[IntVars::cons].array(mfi); + const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - // Base state density - const Array4& r0_arr = r0->const_array(mfi); + // Base state density + const Array4& r0_arr = r0->const_array(mfi); - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = grav_gpu[2] * 0.5 * ( (cell_data(i,j,k ) - r0_arr(i,j,k )) - +(cell_data(i,j,k-1) - r0_arr(i,j,k-1)) ); - }); - } // mfi + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + buoyancy_fab(i, j, k) = buoyancy_type1(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data); + }); + } // mfi - } - else // (buoyancy_type != 1) - { - // We now use the base state rather than planar average because - // 1) we don't want to average over the limited region of the fine level if doing multilevel. - // 2) it's cheaper to use the base state than to compute the horizontal averages - // 3) when running in a smallish domain, the horizontal average may evolve over time, - // which is not necessarily the intended behavior - // + } + else // (buoyancy_type != 1) + { + // We now use the base state rather than planar average because + // 1) we don't want to average over the limited region of the fine level if doing multilevel. + // 2) it's cheaper to use the base state than to compute the horizontal averages + // 3) when running in a smallish domain, the horizontal average may evolve over time, + // which is not necessarily the intended behavior + // #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); + 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 boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); + // We don't compute a source term for z-momentum on the bottom or top boundary + if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); + if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - // Base state density and pressure - const Array4& r0_arr = r0->const_array(mfi); - const Array4& p0_arr = p0->const_array(mfi); + // Base state density and pressure + const Array4& r0_arr = r0->const_array(mfi); + const Array4& p0_arr = p0->const_array(mfi); - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); + const Array4 & cell_data = S_data[IntVars::cons].array(mfi); + const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k)); - Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp); - Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); - Real qplus = (t_hi-t0_hi)/t0_hi; - - Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)); - Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp); - Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - Real qminus = (t_lo-t0_lo)/t0_lo; - - Real r0_q_avg = Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus); - buoyancy_fab(i, j, k) = -r0_q_avg * grav_gpu[2]; - }); - } // mfi - } // buoyancy_type + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + buoyancy_fab(i, j, k) = buoyancy_dry_default(i,j,k, + grav_gpu[2],rd_over_cp, + p0_arr,r0_arr,cell_data); + }); + } // mfi + } // buoyancy_type } // moisture type else { @@ -196,23 +191,15 @@ void make_buoyancy (Vector& S_data, 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); - - for (int q_offset(2); q_offset& S_data, const Array4 & cell_prim = S_prim.array(mfi); // TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp) - ParallelFor(tbz, [=, buoyancy_type=solverChoice.buoyancy_type] AMREX_GPU_DEVICE (int i, int j, int k) - { - Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ], qv_d_ptr[k ]); - Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]); - - 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)); - 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)); - - Real qplus, qminus; - - Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; - Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; - - Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; - Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; - - Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; - Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; - - if (buoyancy_type == 2) { - 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]); - - 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]); - - } else { // (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_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); - Real r0avg = Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); - - buoyancy_fab(i, j, k) = -qavg * r0avg * grav_gpu[2]; - }); + if (solverChoice.buoyancy_type == 2) { + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + buoyancy_fab(i, j, k) = buoyancy_type2(i,j,k,n_qstate,grav_gpu[2], + rho_d_ptr,theta_d_ptr, + qv_d_ptr,qc_d_ptr,qp_d_ptr, + cell_prim,cell_data); + }); + } else { + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + buoyancy_fab(i, j, k) = buoyancy_type4(i,j,k,n_qstate,grav_gpu[2], + rho_d_ptr,theta_d_ptr, + qv_d_ptr,qc_d_ptr,qp_d_ptr, + cell_prim,cell_data); + }); + } } // mfi } else if (solverChoice.buoyancy_type == 3) { @@ -344,33 +297,9 @@ void make_buoyancy (Vector& S_data, 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 ], qv_d_ptr[k ]); - Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]); - - 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)); - 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)); - - Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; - Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; - - Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; - Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; - - Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; - Real qp_minus = (n_qstate >= 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_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]); - - buoyancy_fab(i, j, k) = -qavg * r0avg * grav_gpu[2]; + buoyancy_fab(i, j, k) = buoyancy_type3(i,j,k,n_qstate,grav_gpu[2], + rho_d_ptr,theta_d_ptr,qv_d_ptr, + cell_prim,cell_data); }); } // mfi } // buoyancy_type diff --git a/Source/SourceTerms/Make.package b/Source/SourceTerms/Make.package index 32cf2443a..2ed03c5a6 100644 --- a/Source/SourceTerms/Make.package +++ b/Source/SourceTerms/Make.package @@ -12,4 +12,4 @@ endif CEXE_headers += ERF_NumericalDiffusion.H CEXE_headers += ERF_Src_headers.H - +CEXE_headers += ERF_buoyancy_utils.H