Skip to content

Commit

Permalink
Fix bugs.
Browse files Browse the repository at this point in the history
  • Loading branch information
Aaron Lattanzi committed Sep 26, 2024
1 parent e2cf47a commit 9a83c9f
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 35 deletions.
73 changes: 38 additions & 35 deletions Source/Initialization/ERF_init_from_metgrid.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
/**
* \file ERF_init_from_metgrid.cpp
*/

#include <ERF_Constants.H>
#include <ERF_Metgrid_utils.H>

using namespace amrex;
Expand Down Expand Up @@ -847,6 +847,9 @@ init_base_state_from_metgrid (const bool use_moisture,
Arena_Used = The_Pinned_Arena();
#endif

// Copy to access on GPU
Real grav = CONST_GRAV;

{ // set pressure and density at initialization.
const Array4<Real>& r_hse_arr = r_hse_fab.array();
const Array4<Real>& p_hse_arr = p_hse_fab.array();
Expand All @@ -863,7 +866,7 @@ init_base_state_from_metgrid (const bool use_moisture,
auto const new_z = z_phys_cc_fab.const_array();

const int maxiter = 10;
const amrex::Real tol = 1.0e-12;
const Real tol = 1.0e-12;
ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
{
for (int k=0; k<=kmax; k++) {
Expand All @@ -874,22 +877,22 @@ init_base_state_from_metgrid (const bool use_moisture,
z_vec[kmax+1] = new_z(i,j,kmax+1);

// Get Psurf
amrex::Real Psurf;
Real Psurf;
if (flag_psfc_vec[0] == 1) {
Psurf = orig_psfc(i,j,0);
} else {
amrex::Real t_0 = 290.0; // WRF's model_config_rec%base_temp
amrex::Real a = 50.0; // WRF's model_config_rec%base_lapse
Psurf = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*CONST_GRAV*z_vec[0]/(a*R_d)), 0.5));
Real t_0 = 290.0; // WRF's model_config_rec%base_temp
Real a = 50.0; // WRF's model_config_rec%base_lapse
Psurf = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*grav*z_vec[0]/(a*R_d)), 0.5));
}

// Iterations for the first CC point that is 1/2 dz off the surface
{
amrex::Real half_dz = z_vec[0];
amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[0];
Real half_dz = z_vec[0];
Real qvf = 1.0+(R_v/R_d)*Q_vec[0];
Thetam_vec[0] = Thetad_vec[0]*qvf;
for (int iter=0; iter<maxiter; iter++) {
Pm_vec[0] = Psurf-half_dz*(Rhom_vec[0])*(1.0+Q_vec[0])*CONST_GRAV;
Pm_vec[0] = Psurf-half_dz*(Rhom_vec[0])*(1.0+Q_vec[0])*grav;
if (Pm_vec[0] < 0.0) Pm_vec[0] = 0.0;
Rhom_vec[0] = (p_0/(R_d*Thetam_vec[0]))*std::pow(Pm_vec[0]/p_0, iGamma);
} // it
Expand All @@ -899,28 +902,28 @@ init_base_state_from_metgrid (const bool use_moisture,
for (int k=1; k<=kmax; k++) {

// Assignment only for downward integration in following loop
amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[k];
Real qvf = 1.0+(R_v/R_d)*Q_vec[k];
Thetam_vec[k] = Thetad_vec[k]*qvf;

// Vertical grid spacing
amrex::Real dz = z_vec[k]-z_vec[k-1];
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;
Real rho_tot_lo = Rhom_vec[k-1] * (1. + Q_vec[k-1]);
Real C = -Pm_vec[k-1] + 0.5*rho_tot_lo*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;
Real rho_tot_hi = Rhom_vec[k] * (1. + Q_vec[k]);
Real F = Pm_vec[k] + 0.5*rho_tot_hi*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],
grav, C, Thetad_vec[k],
Q_vec[k], Q_vec[k],
Pm_vec[k], Rhom_vec[k], F);
} // k
Expand All @@ -929,10 +932,10 @@ init_base_state_from_metgrid (const bool use_moisture,
Pd_vec[kmax] = Pm_vec[kmax];
Rhod_vec[kmax] = (p_0/(R_d*Thetam_vec[kmax]))*std::pow(Pd_vec[kmax]/p_0, iGamma);
for (int k=kmax-1; k>=0; k--) {
amrex::Real dz = z_vec[k+1]-z_vec[k];
Real dz = z_vec[k+1]-z_vec[k];
Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess.
for (int iter=0; iter<maxiter; iter++) {
Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV;
Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*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);
} // it
Expand Down Expand Up @@ -1030,7 +1033,7 @@ init_base_state_from_metgrid (const bool use_moisture,
auto p_hse_arr = p_hse_bcs_fab.array();

const int maxiter = 10;
const amrex::Real tol = 1.0e-12;
const Real tol = 1.0e-12;
ParallelFor(valid_bx2d, [=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
{
for (int k=0; k<=kmax; k++) {
Expand All @@ -1041,22 +1044,22 @@ init_base_state_from_metgrid (const bool use_moisture,
z_vec[kmax+1] = new_z(i,j,kmax+1);

// Get Psurf
amrex::Real Psurf;
Real Psurf;
if (flag_psfc_vec[it] == 1) {
Psurf = orig_psfc(i,j,0);
} else {
amrex::Real t_0 = 290.0; // WRF's model_config_rec%base_temp
amrex::Real a = 50.0; // WRF's model_config_rec%base_lapse
Psurf = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*CONST_GRAV*z_vec[0]/(a*R_d)), 0.5));
Real t_0 = 290.0; // WRF's model_config_rec%base_temp
Real a = 50.0; // WRF's model_config_rec%base_lapse
Psurf = p_0*exp(-t_0/a+std::pow((std::pow(t_0/a, 2)-2.0*grav*z_vec[0]/(a*R_d)), 0.5));
}

// Iterations for the first CC point that is 1/2 dz off the surface
{
amrex::Real half_dz = z_vec[0];
amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[0];
Real half_dz = z_vec[0];
Real qvf = 1.0+(R_v/R_d)*Q_vec[0];
Thetam_vec[0] = Thetad_vec[0]*qvf;
for (int iter=0; iter<maxiter; iter++) {
Pm_vec[0] = Psurf-half_dz*(Rhom_vec[0])*(1.0+Q_vec[0])*CONST_GRAV;
Pm_vec[0] = Psurf-half_dz*(Rhom_vec[0])*(1.0+Q_vec[0])*grav;
if (Pm_vec[0] < 0.0) Pm_vec[0] = 0.0;
Rhom_vec[0] = (p_0/(R_d*Thetam_vec[0]))*std::pow(Pm_vec[0]/p_0, iGamma);
} // it
Expand All @@ -1066,28 +1069,28 @@ init_base_state_from_metgrid (const bool use_moisture,
for (int k=1; k<=kmax; k++) {

// Assignment only for downward integration in following loop
amrex::Real qvf = 1.0+(R_v/R_d)*Q_vec[k];
Real qvf = 1.0+(R_v/R_d)*Q_vec[k];
Thetam_vec[k] = Thetad_vec[k]*qvf;

// Vertical grid spacing
amrex::Real dz = z_vec[k]-z_vec[k-1];
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;
Real rho_tot_lo = Rhom_vec[k-1] * (1. + Q_vec[k-1]);
Real C = -Pm_vec[k-1] + 0.5*rho_tot_lo*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;
Real rho_tot_hi = Rhom_vec[k] * (1. + Q_vec[k]);
Real F = Pm_vec[k] + 0.5*rho_tot_hi*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],
grav, C, Thetad_vec[k],
Q_vec[k], Q_vec[k],
Pm_vec[k], Rhom_vec[k], F);
} // k
Expand All @@ -1096,10 +1099,10 @@ init_base_state_from_metgrid (const bool use_moisture,
Pd_vec[kmax] = Pm_vec[kmax];
Rhod_vec[kmax] = (p_0/(R_d*Thetam_vec[kmax]))*std::pow(Pd_vec[kmax]/p_0, iGamma);
for (int k=kmax-1; k>=0; k--) {
amrex::Real dz = z_vec[k+1]-z_vec[k];
Real dz = z_vec[k+1]-z_vec[k];
Rhod_vec[k] = Rhod_vec[k+1]; // an initial guess.
for (int iter=0; iter<maxiter; iter++) {
Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*CONST_GRAV;
Pd_vec[k] = Pd_vec[k+1]+0.5*dz*(Rhod_vec[k]+Rhod_vec[k+1])*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);
} // it
Expand Down
2 changes: 2 additions & 0 deletions Source/Utils/ERF_HSE_utils.H
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
#ifndef ERF_HSEUTIL_H_
#define ERF_HSEUTIL_H_

#include <ERF_Constants.H>

/**
* Utility functions for calculating a hydrostatic equilibrium (HSE) base state
*/
Expand Down

0 comments on commit 9a83c9f

Please sign in to comment.