Skip to content

Commit

Permalink
Merge branch 'development' into SAM_Cloud
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi authored Sep 26, 2024
2 parents 522ed91 + ce5fd12 commit 2a8ba24
Show file tree
Hide file tree
Showing 6 changed files with 117 additions and 34 deletions.
41 changes: 26 additions & 15 deletions Source/DataStructs/ERF_InputSoundingData.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include <ERF_EOS.H>
#include <ERF_Constants.H>
#include <ERF_Interpolation_1D.H>
#include <ERF_HSE_utils.H>

/**
* Data structure storing input sounding data. Also
Expand Down Expand Up @@ -179,7 +180,7 @@ public:
* from the surface up through the air column to get the dry density
* and moist pressure.
*/
const int maxiter = 10;
const amrex::Real tol = 1.0e-12;
const int Ninp = size(itime);
pm_integ.resize(Ninp);
rhod_integ.resize(Ninp);
Expand All @@ -199,8 +200,6 @@ public:
// p_tot = rho_m R_d T_v
// = rho_d(1 + q_v) R_d T_v

// integrate from surface to domain top
amrex::Real dz;
#if 0 // Printing
// In this absence of moisture, this moist profile will match the
// following dry profile
Expand All @@ -214,21 +213,33 @@ public:
<< " " << V_inp_sound[itime][0]
<< std::endl;
#endif

// integrate from surface to domain top
amrex::Real dz, F, C;
amrex::Real rho_tot_hi, rho_tot_lo;
for (int k=1; k < size(itime); ++k)
{
// Vertical grid spacing
dz = z_inp_sound[itime][k] - z_inp_sound[itime][k-1];
rhod_integ[k] = rhod_integ[k-1]; // guess
for (int it=0; it < maxiter; ++it)
{
amrex::Real rho_tot_hi = rhod_integ[k ] * (1. + qv_inp_sound[itime][k ]);
amrex::Real rho_tot_lo = rhod_integ[k-1] * (1. + qv_inp_sound[itime][k-1]);
pm_integ[k] = pm_integ[k-1] - 0.5*dz*(rho_tot_hi + rho_tot_lo)*CONST_GRAV;
AMREX_ALWAYS_ASSERT(pm_integ[k] > 0);
rhod_integ[k] = getRhogivenThetaPress(theta_inp_sound[itime][k],
pm_integ[k],
R_d/Cp_d,
qv_inp_sound[itime][k]);
}

// Establish known constant
rho_tot_lo = rhod_integ[k-1] * (1. + qv_inp_sound[itime][k-1]);
C = -pm_integ[k-1] + 0.5*rho_tot_lo*CONST_GRAV*dz;

// Initial guess and residual
pm_integ[k] = pm_integ[k-1];
rhod_integ[k] = getRhogivenThetaPress(theta_inp_sound[itime][k],
pm_integ[k],
R_d/Cp_d,
qv_inp_sound[itime][k]);
rho_tot_hi = rhod_integ[k] * (1. + qv_inp_sound[itime][k]);
F = pm_integ[k] + 0.5*rho_tot_hi*CONST_GRAV*dz + C;

// Do iterations
if (std::abs(F)>tol) HSEutils::Newton_Raphson_hse(tol, R_d/Cp_d, dz,
CONST_GRAV, C, theta_inp_sound[itime][k],
qv_inp_sound[itime][k], qv_inp_sound[itime][k],
pm_integ[k], rhod_integ[k], F);
#if 0 // Printing
amrex::Print() << z_inp_sound[itime][k]
<< " " << pm_integ[k]
Expand Down
39 changes: 28 additions & 11 deletions Source/Initialization/ERF_Metgrid_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include <ERF_Constants.H>
#include <ERF_Utils.H>
#include <ERF_prob_common.H>
#include <ERF_HSE_utils.H>

void
read_from_metgrid (int lev,
Expand Down Expand Up @@ -116,6 +117,7 @@ calc_rho_p (const int& kmax,
amrex::Real* Pm_vec)
{
const int maxiter = 10;
const amrex::Real tol = 1.0e-12;

// Calculate or use moist pressure at the surface.
amrex::Real Psurf;
Expand All @@ -141,17 +143,32 @@ calc_rho_p (const int& kmax,

// Integrate from the first CC point to the top boundary.
for (int k=1; k<=kmax; k++) {
amrex::Real dz = z_vec[k]-z_vec[k-1];

// Assignment only for downward integration in following loop
amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[k];
Thetam_vec[k] = Thetad_vec[k]*qvf;
Rhom_vec[k] = Rhom_vec[k-1]; // an initial guess.
for (int it=0; it<maxiter; it++) {
amrex::Real Rho_tot_hi = Rhom_vec[k ] * (1.0+Q_vec[k ]);
amrex::Real Rho_tot_lo = Rhom_vec[k-1] * (1.0+Q_vec[k-1]);
Pm_vec[k] = Pm_vec[k-1]-0.5*dz*(Rho_tot_hi + Rho_tot_lo)*CONST_GRAV;
if (Pm_vec[k] < 0.0) Pm_vec[k] = 0.0;
Rhom_vec[k] = (p_0/(R_d*Thetam_vec[k]))*std::pow(Pm_vec[k]/p_0, iGamma);
} // it
Thetam_vec[k] = Thetad_vec[k]*qvf;

// Vertical grid spacing
amrex::Real dz = z_vec[k]-z_vec[k-1];

// Establish known constant
amrex::Real rho_tot_lo = Rhom_vec[k-1] * (1. + Q_vec[k-1]);
amrex::Real C = -Pm_vec[k-1] + 0.5*rho_tot_lo*CONST_GRAV*dz;

// Initial guess and residual
Pm_vec[k] = Pm_vec[k-1];
Rhom_vec[k] = getRhogivenThetaPress(Thetad_vec[k],
Pm_vec[k],
R_d/Cp_d,
Q_vec[k]);
amrex::Real rho_tot_hi = Rhom_vec[k] * (1. + Q_vec[k]);
amrex::Real F = Pm_vec[k] + 0.5*rho_tot_hi*CONST_GRAV*dz + C;

// Do iterations
if (std::abs(F)>tol) HSEutils::Newton_Raphson_hse(tol, R_d/Cp_d, dz,
CONST_GRAV, C, Thetad_vec[k],
Q_vec[k], Q_vec[k],
Pm_vec[k], Rhom_vec[k], F);
} // k

// integrate from the top back down to get dry pressure and density.
Expand All @@ -160,7 +177,7 @@ calc_rho_p (const int& kmax,
for (int k=kmax-1; k>=0; k--) {
amrex::Real dz = z_vec[k+1]-z_vec[k];
Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess.
for (int it=0; it < maxiter; it++) {
for (int it=0; it<maxiter; it++) {
Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV;
if (Pd_vec[k] < 0.0) Pd_vec[k] = 0.0;
Rhod_vec[k] = (p_0/(R_d*Thetam_vec[k]))*std::pow(Pd_vec[k]/p_0, iGamma);
Expand Down
55 changes: 55 additions & 0 deletions Source/Utils/ERF_HSE_utils.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,60 @@ namespace HSEutils
const int MAX_ITER = 10;
const amrex::Real TOL = 1.e-8;

/**
* Function to calculate the hydrostatic density and pressure
* from Newton-Raphson iterations
*
* @param[in] m_tol iteration tolerance
* @param[in] RdOCp Rd/Cp
* @param[out] dz change in vertical height
* @param[out] g magnitude of gravity
* @param[in] C sum of known terms in HSE balance
* @param[in] T theta at the current cell center
* @param[in] qt total moisture (non-precip and precip)
* @param[in] qv water vapor
* @param[in] P pressure at cell center
* @param[in] rd dry density at cell center
* @param[in] F starting residual of non-linear eq
*/
AMREX_GPU_HOST_DEVICE
AMREX_FORCE_INLINE
void
Newton_Raphson_hse (const Real& m_tol,
const Real& RdOCp,
const Real& dz,
const Real& g,
const Real& C,
const Real& T,
const Real& qt,
const Real& qv,
Real& P,
Real& rd,
Real& F)
{
int iter=0;
Real eps = 1.0e-6;
Real ieps = 1.0e6;
do {
// Compute change in pressure
Real hi = getRhogivenThetaPress(T, P + eps, RdOCp, qv);
Real lo = getRhogivenThetaPress(T, P - eps, RdOCp, qv);
Real dFdp = 1.0 + 0.25*ieps*(hi - lo)*g*dz;
P -= F/dFdp;
AMREX_ALWAYS_ASSERT(P > 0.0);

// Diagnose density and residual
rd = getRhogivenThetaPress(T, P, RdOCp, qv);
AMREX_ALWAYS_ASSERT(rd > 0.0);
Real r_tot = rd * (1. + qt);
F = P + 0.5*r_tot*g*dz + C;
++iter;
}
while (std::abs(F)>m_tol && iter<25);
if (iter>=25) printf("WARNING: HSE Newton iterations did not converge to tolerance!\n");
}


/**
* Function to calculate the hydrostatic density and pressure
*
Expand Down Expand Up @@ -123,6 +177,7 @@ namespace HSEutils
r[khi+1] = r[khi];
}


/**
* Function to calculate the hydrostatic density and pressure over terrain
*
Expand Down
Binary file modified Tests/ERFGoldFiles/ABL_MYNN_PBL/Level_0/Cell_D_00000
View file
Binary file not shown.
4 changes: 2 additions & 2 deletions Tests/ERFGoldFiles/ABL_MYNN_PBL/Level_0/Cell_H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,8 +9,8 @@
FabOnDisk: Cell_D_00000 0

1,9
1.26070047102876259e+00,9.99999999999999980e-13,4.46574933942103414e+00,-1.12232652875381444e-08,-7.69073596531594396e-07,2.64998234491973051e+02,9.57655146925804438e+04,9.99044297871173762e-12,1.51954705865492993e-12,
1.26070047102875904e+00,9.99999999999999980e-13,4.46574933952489417e+00,-1.12232648780406732e-08,-7.69073595650813326e-07,2.64998234491972653e+02,9.57655146925807348e+04,9.99044297877987339e-12,1.51954705866529331e-12,

1,9
1.32197155095610142e+00,1.22162865379599972e+00,8.00000067867113174e+00,2.08085771756629026e-02,1.28715032183049716e-05,2.67968750004691060e+02,1.00759448324675424e+05,4.06335629766707118e+00,3.80301353815072885e+00,
1.32197155095610364e+00,1.22162865373526963e+00,8.00000067867087772e+00,2.08085771745764106e-02,1.28715031999786871e-05,2.67968750004692424e+02,1.00759448324675468e+05,4.06335629715192681e+00,3.80301353767852346e+00,

12 changes: 6 additions & 6 deletions Tests/ERFGoldFiles/ABL_MYNN_PBL/job_info
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,19 +5,19 @@ inputs file: ABL_MYNN_PBL.i

number of MPI processes: 1

CPU time used since start of simulation (CPU-hours): 0.000230912
CPU time used since start of simulation (CPU-hours): 0.000225806

===============================================================================
Plotfile Information
===============================================================================
output data / time: Fri Sep 13 11:17:15 2024
output dir: /home/alattanz/git/ERF/Exec/ABL/HOLD
output data / time: Tue Sep 24 15:03:42 2024
output dir: /home/alattanz/git/ERF/Exec/ABL


===============================================================================
Build Information
===============================================================================
build date: 2024-09-13 11:05:52.807635
build date: 2024-09-24 14:52:58.564129
build machine: Linux LAPTOP-PO22H9R0 4.4.0-22621-Microsoft #3672-Microsoft Fri Jan 01 08:00:00 PST 2016 x86_64 x86_64 x86_64 GNU/Linux
build dir: /home/alattanz/git/ERF/Exec/ABL
AMReX dir: /home/alattanz/git/amrex
Expand All @@ -26,8 +26,8 @@ COMP: gnu
COMP version: 9.4.0


ERF git hash: bea13e5a
AMReX git hash: 22.05-998-g4460afbbc
ERF git hash: dd122aab-dirty
AMReX git hash: 22.05-1011-g1bfdb84c1


===============================================================================
Expand Down

0 comments on commit 2a8ba24

Please sign in to comment.