HSEutils Namespace Reference

Functions
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)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	init_isentropic_hse (const amrex::Real &r_sfc, const amrex::Real &theta, amrex::Real *r, amrex::Real *p, const amrex::Real &dz, const int klo, const int khi)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	init_isentropic_hse_terrain (int i, int j, const amrex::Real &r_sfc, const amrex::Real &theta, amrex::Real *r, amrex::Real *p, const amrex::Array4< amrex::Real const > z_cc, const int &klo, const int &khi)

Variables
const int	MAX_ITER = 10
const amrex::Real	TOL = 1.e-8

Detailed Description

Utility functions for calculating a hydrostatic equilibrium (HSE) base state

Function Documentation

init_isentropic_hse()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void HSEutils::init_isentropic_hse	(	const amrex::Real &	r_sfc,
		const amrex::Real &	theta,
		amrex::Real *	r,
		amrex::Real *	p,
		const amrex::Real &	dz,
		const int	klo,
		const int	khi
	)

Function to calculate the hydrostatic density and pressure

Parameters

[in]	r_sfc	surface density
[in]	theta	surface potential temperature
[out]	r	hydrostatically balanced density profile
[out]	p	hydrostatically balanced pressure profile
[in]	dz	vertical grid spacing (constant)
[in]	prob_lo_z	surface height
[in]	khi	z-index corresponding to the big end of the domain

    {
      int kstart;
 
      // r_sfc / p_0 are the density / pressure at klo
 
      // Initial guess
      Real half_dz = 0.5*dz;
      r[klo] = r_sfc;
      p[klo] = p_0 - half_dz * r[klo] * CONST_GRAV;
 
      if (klo == 0)
      {
          kstart = 1;
 
          // We do a Newton iteration to satisfy the EOS & HSE (with constant theta)
          bool converged_hse = false;
          Real p_hse;
          Real p_eos;
 
          for (int iter = 0; iter < MAX_ITER && !converged_hse; iter++)
          {
              p_hse = p_0 - half_dz * r[klo] * CONST_GRAV;
              p_eos = getPgivenRTh(r[klo]*theta);
 
              Real A = p_hse - p_eos;
 
              Real dpdr = getdPdRgivenConstantTheta(r[klo],theta);
 
              Real drho = A / (dpdr + half_dz * CONST_GRAV);
 
              r[klo] = r[klo] + drho;
              p[klo] = getPgivenRTh(r[klo]*theta);
 
              if (std::abs(drho) < TOL)
              {
                  converged_hse = true;
                  break;
              }
          }
 
          // if (!converged_hse) amrex::Print() << "DOING ITERATIONS AT K = " << klo << std::endl;
          // if (!converged_hse) amrex::Error("Didn't converge the iterations in init");
      } else {
          kstart = klo; // because we need to solve for r[klo] here; r[klo-1] is what was passed in r_sfc
          r[klo-1] = r_sfc; // r_sfc is really the interpolated density in the ghost cell in this case
          p[klo-1] = getPgivenRTh(r[klo-1]*theta);
      }
 
      // To get values at k > 0 we do a Newton iteration to satisfy the EOS (with constant theta) and
      for (int k = kstart; k <= khi; k++)
      {
          // To get values at k > 0 we do a Newton iteration to satisfy the EOS (with constant theta) and
          // to discretely satisfy HSE -- here we assume spatial_order = 2 -- we can generalize this later if needed
          bool converged_hse = false;
 
          r[k] = r[k-1];
 
          Real p_eos = getPgivenRTh(r[k]*theta);
          Real p_hse;
 
          for (int iter = 0; iter < MAX_ITER && !converged_hse; iter++)
          {
              Real r_avg = 0.5 * (r[k-1]+r[k]);
              p_hse = p[k-1] -  dz * r_avg * CONST_GRAV;
              p_eos = getPgivenRTh(r[k]*theta);
 
              Real A = p_hse - p_eos;
 
              Real dpdr = getdPdRgivenConstantTheta(r[k],theta);
              // Gamma * p_0 * std::pow( (R_d * theta / p_0), Gamma) * std::pow(r[k], Gamma-1.0) ;
 
              Real drho = A / (dpdr + dz * CONST_GRAV);
 
              r[k] = r[k] + drho;
              p[k] = getPgivenRTh(r[k]*theta);
 
              if (std::abs(drho) < TOL * r[k-1])
              {
                  converged_hse = true;
                  // amrex::Print() << " converged " << std::endl;
                  break;
              }
          }
 
          // if (!converged_hse) amrex::Print() << "DOING ITERATIONS AT K = " << k << std::endl;
          // if (!converged_hse) amrex::Error("Didn't converge the iterations in init");
      }
      r[khi+1] = r[khi];
    }

Referenced by erf_init_dens_hse().

Here is the call graph for this function:

Here is the caller graph for this function:

init_isentropic_hse_terrain()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void HSEutils::init_isentropic_hse_terrain	(	int	i,
		int	j,
		const amrex::Real &	r_sfc,
		const amrex::Real &	theta,
		amrex::Real *	r,
		amrex::Real *	p,
		const amrex::Array4< amrex::Real const >	z_cc,
		const int &	klo,
		const int &	khi
	)

Function to calculate the hydrostatic density and pressure over terrain

Parameters

[in]	i	x-index
[in]	j	y-index
[in]	r_sfc	surface density
[in]	theta	surface potential temperature
[out]	r	hydrostatically balanced density profile
[out]	p	hydrostatically balanced pressure profile
[in]	z_cc	cell-center heights
[in]	khi	z-index corresponding to the big end of the domain

{
    int kstart;
 
    if (klo == 0)  {
        //
        // r_sfc / p_0 are the density / pressure at the surface
        //
        // Initial guess
        int k0 = 0;
 
        // Where we start the lower iteration
        kstart = 1;
 
        Real half_dz = z_cc(i,j,k0);
        r[k0] = r_sfc;
        p[k0] = p_0 - half_dz * r[k0] * CONST_GRAV;
        {
            // We do a Newton iteration to satisfy the EOS & HSE (with constant theta)
            bool converged_hse = false;
            Real p_hse;
            Real p_eos;
 
            for (int iter = 0; iter < MAX_ITER && !converged_hse; iter++)
            {
                p_hse = p_0 - half_dz * r[k0] * CONST_GRAV;
                p_eos = getPgivenRTh(r[k0]*theta);
 
                Real A = p_hse - p_eos;
 
                Real dpdr = getdPdRgivenConstantTheta(r[k0],theta);
 
                Real drho = A / (dpdr + half_dz * CONST_GRAV);
 
                r[k0] = r[k0] + drho;
                p[k0] = getPgivenRTh(r[k0]*theta);
 
                if (std::abs(drho) < TOL)
                {
                    converged_hse = true;
                    break;
                }
            }
 
            //if (!converged_hse) amrex::Print() << "DOING ITERATIONS AT K = " << k0 << std::endl;
            //if (!converged_hse) amrex::Error("Didn't converge the iterations in init");
        }
    } else {
        kstart = klo; // because we need to solve for r[klo] here; r[klo-1] is what was passed in r_sfc
        r[klo-1] = r_sfc; // r_sfc is really the interpolated density in the ghost cell in this case
        p[klo-1] = getPgivenRTh(r[klo-1]*theta);
    }
 
     // To get values at k > 0 we do a Newton iteration to satisfy the EOS (with constant theta) and
     for (int k = kstart; k <= khi; k++)
     {
         // To get values at k > 0 we do a Newton iteration to satisfy the EOS (with constant theta) and
         // to discretely satisfy HSE -- here we assume spatial_order = 2 -- we can generalize this later if needed
         bool converged_hse = false;
 
         Real dz_loc = (z_cc(i,j,k) - z_cc(i,j,k-1));
 
         r[k] = r[k-1];
 
         Real p_eos = getPgivenRTh(r[k]*theta);
         Real p_hse;
 
         for (int iter = 0; iter < MAX_ITER && !converged_hse; iter++)
         {
             p_hse = p[k-1] -  dz_loc * 0.5 * (r[k-1]+r[k]) * CONST_GRAV;
             p_eos = getPgivenRTh(r[k]*theta);
 
             Real A = p_hse - p_eos;
 
             Real dpdr = getdPdRgivenConstantTheta(r[k],theta);
 
             Real drho = A / (dpdr + dz_loc * CONST_GRAV);
 
             r[k] = r[k] + drho;
             p[k] = getPgivenRTh(r[k]*theta);
 
             if (std::abs(drho) < TOL * r[k-1])
             {
                 converged_hse = true;
                 //amrex::Print() << " converged " << std::endl;
                 break;
             }
         }
 
         // if (!converged_hse) amrex::Print() << "DOING ITERATIONS AT K = " << k << std::endl;
         // if (!converged_hse) amrex::Error("Didn't converge the iterations in init");
    }
}

Referenced by erf_init_dens_hse().

Here is the call graph for this function:

Here is the caller graph for this function:

Newton_Raphson_hse()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void HSEutils::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
	)

Function to calculate the hydrostatic density and pressure from Newton-Raphson iterations

Parameters

[in]	m_tol	iteration tolerance
[in]	RdOCp	Rd/Cp
[out]	dz	change in vertical height
[out]	g	magnitude of gravity
[in]	C	sum of known terms in HSE balance
[in]	T	theta at the current cell center
[in]	qt	total moisture (non-precip and precip)
[in]	qv	water vapor
[in]	P	pressure at cell center
[in]	rd	dry density at cell center
[in]	F	starting residual of non-linear eq

    {
        int iter=0;
        int max_iter=30;
        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;
            if (P < 1.0e3) amrex::Warning("P < 1000 [Pa]; Domain height may be too large...");
            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<max_iter);
        if (iter>=max_iter) {
            printf("WARNING: HSE Newton iterations did not converge to tolerance!\n");
            printf("HSE Newton tol: %e %e\n",F,m_tol);
        }
    }

Referenced by InputSoundingData::calc_rho_p().

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

MAX_ITER

const int HSEutils::MAX_ITER = 10

Referenced by init_isentropic_hse(), and init_isentropic_hse_terrain().

TOL

const amrex::Real HSEutils::TOL = 1.e-8

Referenced by init_isentropic_hse(), and init_isentropic_hse_terrain().