From 72bf21232d655e574176043b5e608cea90adad24 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Tue, 5 Nov 2024 16:51:46 -0800 Subject: [PATCH] more cleanup; no change in functionality expected (#1930) --- CMakeLists.txt | 2 +- Exec/MoistRegTests/Bubble/ERF_prob.cpp | 40 ++----- Exec/MoistRegTests/SquallLine_2D/ERF_prob.cpp | 46 ++------ Exec/MoistRegTests/SuperCell_3D/ERF_prob.cpp | 42 ++----- Source/ERF_make_new_level.cpp | 6 - ...TI_fast_rhs_fun.H => ERF_TI_substep_fun.H} | 0 Source/TimeIntegration/ERF_advance_dycore.cpp | 2 +- Source/TimeIntegration/Make.package | 2 +- Source/Utils/ERF_EOS.H | 103 ++++++++++++------ 9 files changed, 104 insertions(+), 139 deletions(-) rename Source/TimeIntegration/{ERF_TI_fast_rhs_fun.H => ERF_TI_substep_fun.H} (100%) diff --git a/CMakeLists.txt b/CMakeLists.txt index e045d2eb6..3b8d6b340 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -247,7 +247,7 @@ target_link_libraries(erf_api PUBLIC erf_srclib) add_library(${PROJECT_NAME}::erf_api ALIAS erf_srclib) # Collect all headers and make them installable with the target -set(ERF_INCLUDES "Source/ERF.H;Source/ERF_Constants.H;Source/WindFarmParametrization/SimpleActuatorDisk/ERF_SimpleAD.H;Source/WindFarmParametrization/EWP/ERF_EWP.H;Source/WindFarmParametrization/Null/ERF_NullWindFarm.H;Source/WindFarmParametrization/ERF_WindFarm.H;Source/WindFarmParametrization/Fitch/ERF_Fitch.H;Source/BoundaryConditions/ERF_PhysBCFunct.H;Source/BoundaryConditions/ERF_MOSTAverage.H;Source/BoundaryConditions/ERF_MOSTRoughness.H;Source/BoundaryConditions/ERF_ABLMost.H;Source/BoundaryConditions/ERF_FillPatcher.H;Source/BoundaryConditions/ERF_MOSTStress.H;Source/BoundaryConditions/ERF_TimeInterpolatedData.H;Source/Utils/ERF_Interpolation.H;Source/Utils/ERF_TileNoZ.H;Source/Utils/ERF_PlaneAverage.H;Source/Utils/ERF_Interpolation_WENO.H;Source/Utils/ERF_DirectionSelector.H;Source/Utils/ERF_ParFunctions.H;Source/Utils/ERF_Wstar.H;Source/Utils/ERF_Microphysics_Utils.H;Source/Utils/ERF_Sat_methods.H;Source/Utils/ERF_Interpolation_1D.H;Source/Utils/ERF_Interpolation_UPW.H;Source/Utils/ERF_TerrainMetrics.H;Source/Utils/ERF_Interpolation_WENO_Z.H;Source/Utils/ERF_Thetav.H;Source/Utils/ERF_Water_vapor_saturation.H;Source/Utils/ERF_Utils.H;Source/Utils/ERF_Orbit.H;Source/Utils/ERF_EOS.H;Source/Utils/ERF_HSE_utils.H;Source/EB/ERF_TerrainIF.H;Source/EB/ERF_FlowerIF.H;Source/EB/ERF_eb_if.H;Source/Particles/ERFPC.H;Source/Particles/ERF_ParticleData.H;Source/Prob/ERF_init_density_hse_dry.H;Source/Prob/ERF_init_rayleigh_damping.H;Source/Prob/ERF_init_constant_density_hse.H;Source/ERF_prob_common.H;Source/ERF_Derive.H;Source/Radiation/ERF_Mam4_constituents.H;Source/Radiation/ERF_Mam4_aero.H;Source/Radiation/ERF_Optics.H;Source/Radiation/ERF_Modal_aero_wateruptake.H;Source/Radiation/ERF_Cloud_rad_props.H;Source/Radiation/ERF_Phys_prop.H;Source/Radiation/ERF_Radiation.H;Source/Radiation/ERF_Albedo.H;Source/Radiation/ERF_Parameterizations.H;Source/Radiation/ERF_Rad_constants.H;Source/Radiation/ERF_Aero_rad_props.H;Source/Radiation/ERF_m2005_effradius.H;Source/Radiation/ERF_Linear_interpolate.H;Source/Radiation/ERF_Slingo.H;Source/Radiation/ERF_Rrtmgp.H;Source/Radiation/ERF_Ebert_curry.H;Source/SourceTerms/ERF_NumericalDiffusion.H;Source/SourceTerms/ERF_Src_headers.H;Source/IO/ERF_NCInterface.H;Source/IO/ERF_NCWpsFile.H;Source/IO/ERF_NCPlotFile.H;Source/IO/ERF_ReadBndryPlanes.H;Source/IO/ERF_WriteBndryPlanes.H;Source/PBL/ERF_MYNNStruct.H;Source/PBL/ERF_PBLModels.H;Source/PBL/ERF_PBLHeight.H;Source/TimeIntegration/ERF_TI_fast_rhs_fun.H;Source/TimeIntegration/ERF_TI_slow_headers.H;Source/TimeIntegration/ERF_TI_slow_rhs_fun.H;Source/TimeIntegration/ERF_TI_fast_headers.H;Source/TimeIntegration/ERF_TI_utils.H;Source/TimeIntegration/ERF_MRI.H;Source/TimeIntegration/ERF_TI_no_substep_fun.H;Source/LandSurfaceModel/Null/ERF_NullSurf.H;Source/LandSurfaceModel/ERF_LandSurface.H;Source/LandSurfaceModel/MM5/ERF_MM5.H;Source/LandSurfaceModel/SLM/ERF_SLM.H;Source/ERF_IndexDefines.H;Source/Advection/ERF_AdvectionSrcForMom_N.H;Source/Advection/ERF_AdvectionSrcForScalars.H;Source/Advection/ERF_AdvectionSrcForMom_T.H;Source/Advection/ERF_Advection.H;Source/MultiBlock/ERF_MultiBlockContainer.H;Source/Initialization/ERF_Metgrid_utils.H;Source/Diffusion/ERF_EddyViscosity.H;Source/Diffusion/ERF_Diffusion.H;Source/Microphysics/Null/ERF_NullMoistLagrangian.H;Source/Microphysics/Null/ERF_NullMoist.H;Source/Microphysics/ERF_Microphysics.H;Source/Microphysics/ERF_LagrangianMicrophysics.H;Source/Microphysics/ERF_EulerianMicrophysics.H;Source/Microphysics/Kessler/ERF_Kessler.H;Source/Microphysics/SAM/ERF_SAM.H;Source/DataStructs/ERF_InputSpongeData.H;Source/DataStructs/ERF_TurbPertStruct.H;Source/DataStructs/ERF_SpongeStruct.H;Source/DataStructs/ERF_AdvStruct.H;Source/DataStructs/ERF_DataStruct.H;Source/DataStructs/ERF_InputSoundingData.H;Source/DataStructs/ERF_DiffStruct.H;Source/DataStructs/ERF_TurbStruct.H") +set(ERF_INCLUDES "Source/ERF.H;Source/ERF_Constants.H;Source/WindFarmParametrization/SimpleActuatorDisk/ERF_SimpleAD.H;Source/WindFarmParametrization/EWP/ERF_EWP.H;Source/WindFarmParametrization/Null/ERF_NullWindFarm.H;Source/WindFarmParametrization/ERF_WindFarm.H;Source/WindFarmParametrization/Fitch/ERF_Fitch.H;Source/BoundaryConditions/ERF_PhysBCFunct.H;Source/BoundaryConditions/ERF_MOSTAverage.H;Source/BoundaryConditions/ERF_MOSTRoughness.H;Source/BoundaryConditions/ERF_ABLMost.H;Source/BoundaryConditions/ERF_FillPatcher.H;Source/BoundaryConditions/ERF_MOSTStress.H;Source/BoundaryConditions/ERF_TimeInterpolatedData.H;Source/Utils/ERF_Interpolation.H;Source/Utils/ERF_TileNoZ.H;Source/Utils/ERF_PlaneAverage.H;Source/Utils/ERF_Interpolation_WENO.H;Source/Utils/ERF_DirectionSelector.H;Source/Utils/ERF_ParFunctions.H;Source/Utils/ERF_Wstar.H;Source/Utils/ERF_Microphysics_Utils.H;Source/Utils/ERF_Sat_methods.H;Source/Utils/ERF_Interpolation_1D.H;Source/Utils/ERF_Interpolation_UPW.H;Source/Utils/ERF_TerrainMetrics.H;Source/Utils/ERF_Interpolation_WENO_Z.H;Source/Utils/ERF_Thetav.H;Source/Utils/ERF_Water_vapor_saturation.H;Source/Utils/ERF_Utils.H;Source/Utils/ERF_Orbit.H;Source/Utils/ERF_EOS.H;Source/Utils/ERF_HSE_utils.H;Source/EB/ERF_TerrainIF.H;Source/EB/ERF_FlowerIF.H;Source/EB/ERF_eb_if.H;Source/Particles/ERFPC.H;Source/Particles/ERF_ParticleData.H;Source/Prob/ERF_init_density_hse_dry.H;Source/Prob/ERF_init_rayleigh_damping.H;Source/Prob/ERF_init_constant_density_hse.H;Source/ERF_prob_common.H;Source/ERF_Derive.H;Source/Radiation/ERF_Mam4_constituents.H;Source/Radiation/ERF_Mam4_aero.H;Source/Radiation/ERF_Optics.H;Source/Radiation/ERF_Modal_aero_wateruptake.H;Source/Radiation/ERF_Cloud_rad_props.H;Source/Radiation/ERF_Phys_prop.H;Source/Radiation/ERF_Radiation.H;Source/Radiation/ERF_Albedo.H;Source/Radiation/ERF_Parameterizations.H;Source/Radiation/ERF_Rad_constants.H;Source/Radiation/ERF_Aero_rad_props.H;Source/Radiation/ERF_m2005_effradius.H;Source/Radiation/ERF_Linear_interpolate.H;Source/Radiation/ERF_Slingo.H;Source/Radiation/ERF_Rrtmgp.H;Source/Radiation/ERF_Ebert_curry.H;Source/SourceTerms/ERF_NumericalDiffusion.H;Source/SourceTerms/ERF_Src_headers.H;Source/IO/ERF_NCInterface.H;Source/IO/ERF_NCWpsFile.H;Source/IO/ERF_NCPlotFile.H;Source/IO/ERF_ReadBndryPlanes.H;Source/IO/ERF_WriteBndryPlanes.H;Source/PBL/ERF_MYNNStruct.H;Source/PBL/ERF_PBLModels.H;Source/PBL/ERF_PBLHeight.H;Source/TimeIntegration/ERF_TI_substep_fun.H;Source/TimeIntegration/ERF_TI_slow_headers.H;Source/TimeIntegration/ERF_TI_slow_rhs_fun.H;Source/TimeIntegration/ERF_TI_fast_headers.H;Source/TimeIntegration/ERF_TI_utils.H;Source/TimeIntegration/ERF_MRI.H;Source/TimeIntegration/ERF_TI_no_substep_fun.H;Source/LandSurfaceModel/Null/ERF_NullSurf.H;Source/LandSurfaceModel/ERF_LandSurface.H;Source/LandSurfaceModel/MM5/ERF_MM5.H;Source/LandSurfaceModel/SLM/ERF_SLM.H;Source/ERF_IndexDefines.H;Source/Advection/ERF_AdvectionSrcForMom_N.H;Source/Advection/ERF_AdvectionSrcForScalars.H;Source/Advection/ERF_AdvectionSrcForMom_T.H;Source/Advection/ERF_Advection.H;Source/MultiBlock/ERF_MultiBlockContainer.H;Source/Initialization/ERF_Metgrid_utils.H;Source/Diffusion/ERF_EddyViscosity.H;Source/Diffusion/ERF_Diffusion.H;Source/Microphysics/Null/ERF_NullMoistLagrangian.H;Source/Microphysics/Null/ERF_NullMoist.H;Source/Microphysics/ERF_Microphysics.H;Source/Microphysics/ERF_LagrangianMicrophysics.H;Source/Microphysics/ERF_EulerianMicrophysics.H;Source/Microphysics/Kessler/ERF_Kessler.H;Source/Microphysics/SAM/ERF_SAM.H;Source/DataStructs/ERF_InputSpongeData.H;Source/DataStructs/ERF_TurbPertStruct.H;Source/DataStructs/ERF_SpongeStruct.H;Source/DataStructs/ERF_AdvStruct.H;Source/DataStructs/ERF_DataStruct.H;Source/DataStructs/ERF_InputSoundingData.H;Source/DataStructs/ERF_DiffStruct.H;Source/DataStructs/ERF_TurbStruct.H") set_target_properties( erf_srclib PROPERTIES PUBLIC_HEADER "${ERF_INCLUDES}") diff --git a/Exec/MoistRegTests/Bubble/ERF_prob.cpp b/Exec/MoistRegTests/Bubble/ERF_prob.cpp index 87f12a7c0..5edf0e671 100644 --- a/Exec/MoistRegTests/Bubble/ERF_prob.cpp +++ b/Exec/MoistRegTests/Bubble/ERF_prob.cpp @@ -1,6 +1,7 @@ #include "ERF_prob.H" #include "ERF_Microphysics_Utils.H" #include "ERF_Constants.H" +#include "ERF_EOS.H" using namespace amrex; @@ -52,13 +53,6 @@ Real compute_relative_humidity () return 1.0; } -AMREX_FORCE_INLINE -AMREX_GPU_HOST_DEVICE -Real compute_vapor_pressure (const Real p_s, const Real RH) -{ - return p_s*RH; -} - AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Real vapor_mixing_ratio (const Real p_b, const Real T_b, const Real RH) @@ -118,30 +112,20 @@ Real compute_dewpoint_temperature (const Real T_b, const Real RH) return T_dp; } -AMREX_FORCE_INLINE -AMREX_GPU_HOST_DEVICE -Real Problem::compute_theta(const Real T_b, const Real p_b) -{ - return T_b*std::pow(p_0/p_b, R_d/Cp_d); -} - -AMREX_FORCE_INLINE -AMREX_GPU_HOST_DEVICE -Real compute_temperature_from_theta(const Real theta, const Real p) -{ - return theta*std::pow(p/p_0, R_d/Cp_d); -} - AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE void Problem::compute_rho (const Real& pressure, Real& theta, Real& rho, Real& q_v, Real& T_dp, Real& T_b) { T_b = compute_temperature(pressure); - theta = compute_theta(T_b, pressure); + + theta = getThgivenPandT(T_b, pressure, (R_d/Cp_d)); + Real RH = compute_relative_humidity(); q_v = vapor_mixing_ratio(pressure, T_b, RH); - rho = pressure/(R_d*T_b*(1.0 + (R_v/R_d)*q_v)); + + rho = getRhogivenTandPress(T_b, pressure, q_v); rho = rho*(1.0 + parms.qt_init); // q_t = 0.02 a given constant for this test case + T_dp = compute_dewpoint_temperature(T_b, RH); } @@ -228,16 +212,10 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, Vector h_q_v(khi+2); Gpu::DeviceVector d_r(khi+2); - Gpu::DeviceVector d_p(khi+2); - Gpu::DeviceVector d_t(khi+2); - Gpu::DeviceVector d_q_v(khi+2); init_isentropic_hse_no_terrain(h_t.data(), h_r.data(),h_p.data(), h_q_v.data(), dz, khi); Gpu::copyAsync(Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin()); - Gpu::copyAsync(Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin()); - Gpu::copyAsync(Gpu::hostToDevice, h_t.begin(), h_t.end(), d_t.begin()); - Gpu::copyAsync(Gpu::hostToDevice, h_q_v.begin(), h_q_v.end(), d_q_v.begin()); Real* r = d_r.data(); @@ -392,13 +370,13 @@ Problem::init_custom_pert( } theta_total = theta_back[k]*(delta_theta/300.0 + 1); - Real T = compute_temperature_from_theta(theta_total, p_back[k]); + Real T = getTgivenPandTh(theta_total, p_back[k], (R_d/Cp_d)); rho = p_back[k]/(R_d*T*(1.0 + (R_v/R_d)*q_v_back[k])); RH = compute_relative_humidity(); Real q_v_hot = vapor_mixing_ratio(p_back[k], T, RH); // Compute background quantities - Real T_back = compute_temperature_from_theta(theta_back[k], p_back[k]); + Real T_back = getTgivenPandTh(theta_back[k], p_back[k], (R_d/Cp_d)); Real rho_back = p_back[k]/(R_d*T_back*(1.0 + (R_v/R_d)*q_v_back[k])); // This version perturbs rho but not p diff --git a/Exec/MoistRegTests/SquallLine_2D/ERF_prob.cpp b/Exec/MoistRegTests/SquallLine_2D/ERF_prob.cpp index a0d555625..01a96d0b2 100644 --- a/Exec/MoistRegTests/SquallLine_2D/ERF_prob.cpp +++ b/Exec/MoistRegTests/SquallLine_2D/ERF_prob.cpp @@ -69,13 +69,6 @@ Real compute_relative_humidity (const Real z, const Real height, const Real z_tr } } -AMREX_FORCE_INLINE -AMREX_GPU_HOST_DEVICE -Real compute_vapor_pressure (const Real p_s, const Real RH) -{ - return p_s*RH; -} - AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Real vapor_mixing_ratio (const Real z, const Real height, const Real p_b, const Real T_b, const Real RH) @@ -92,13 +85,6 @@ Real vapor_mixing_ratio (const Real z, const Real height, const Real p_b, const } } -AMREX_FORCE_INLINE -AMREX_GPU_HOST_DEVICE -Real compute_temperature (const Real p_b, const Real theta_b) -{ - return theta_b*std::pow(p_b/p_0,R_d/Cp_d); -} - AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Real compute_dewpoint_temperature (const Real T_b, const Real RH) @@ -119,10 +105,10 @@ void Problem::compute_rho (const Real& z, const Real& pressure, Real& theta, Rea { theta = compute_theta(z); - T_b = compute_temperature(pressure, theta); + T_b = getTgivenPandTh(theta, pressure, (R_d/Cp_d)); Real RH = compute_relative_humidity(z, parms.height, parms.z_tr, pressure, T_b); q_v = vapor_mixing_ratio(z, parms.height, pressure, T_b, RH); - rho = pressure/(R_d*T_b*(1.0 + (R_v/R_d)*q_v)); + rho = getRhogivenTandPress(T_b, pressure, q_v); rho = rho*(1.0 + q_v); T_dp = compute_dewpoint_temperature(T_b, RH); } @@ -168,26 +154,19 @@ Problem::init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v const Real& dz, const Real& prob_lo_z, const int& khi) { - - //FILE *file_IC; - //file_IC = fopen("input_sounding_probcpp.txt","w"); Real z, T_b, T_dp; // Compute the quantities at z = 0.5*dz (first cell center) z = prob_lo_z + 0.5*dz; p[0] = p_0; compute_p_k(p[0], p_0, theta[0], r[0], q_v[0], T_dp, T_b, dz, z, 0.0); - //fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[0], r[0], theta[0], q_v[0]); for (int k=1;k<=khi;k++){ z = prob_lo_z + (k+0.5)*dz; p[k] = p[k-1]; compute_p_k(p[k], p[k-1], theta[k], r[k], q_v[k], T_dp, T_b, dz, z, r[k-1]); - //fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[k], r[k], theta[k], q_v[k]); } - //fclose(file_IC); - r[khi+1] = r[khi]; } @@ -239,17 +218,11 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, Vector h_t(khi+2); Vector h_q_v(khi+2); - amrex::Gpu::DeviceVector d_r(khi+2); - amrex::Gpu::DeviceVector d_p(khi+2); - amrex::Gpu::DeviceVector d_t(khi+2); - amrex::Gpu::DeviceVector d_q_v(khi+2); - init_isentropic_hse_no_terrain(h_t.data(), h_r.data(),h_p.data(), h_q_v.data(), dz,prob_lo_z,khi); + amrex::Gpu::DeviceVector d_r(khi+2); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_t.begin(), h_t.end(), d_t.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_q_v.begin(), h_q_v.end(), d_q_v.begin()); Real* r = d_r.data(); @@ -352,15 +325,18 @@ Problem::init_custom_pert ( } theta_total = t[k] + delta_theta; - temperature = compute_temperature(p[k], theta_total); - Real T_b = compute_temperature(p[k], t[k]); + + temperature = getTgivenPandTh(theta_total, p[k], (R_d/Cp_d)); + Real T_b = getTgivenPandTh(t[k] , p[k], (R_d/Cp_d)); + RH = compute_relative_humidity(z, height, z_tr, p[k], T_b); Real q_v_hot = vapor_mixing_ratio(z, height, p[k], T_b, RH); rho = p[k]/(R_d*temperature*(1.0 + (R_v/R_d)*q_v_hot)); // Compute background quantities - Real temperature_back = compute_temperature(p[k], t[k]); - Real T_back = compute_temperature(p[k], t[k]); + Real temperature_back = getTgivenPandTh(t[k], p[k], (R_d/Cp_d)); + Real T_back = getTgivenPandTh(t[k], p[k], (R_d/Cp_d)); + Real RH_back = compute_relative_humidity(z, height, z_tr, p[k], T_back); Real q_v_back = vapor_mixing_ratio(z, height, p[k], T_back, RH_back); Real rho_back = p[k]/(R_d*temperature_back*(1.0 + (R_v/R_d)*q_v_back)); diff --git a/Exec/MoistRegTests/SuperCell_3D/ERF_prob.cpp b/Exec/MoistRegTests/SuperCell_3D/ERF_prob.cpp index 1b596b69a..2aaff6787 100644 --- a/Exec/MoistRegTests/SuperCell_3D/ERF_prob.cpp +++ b/Exec/MoistRegTests/SuperCell_3D/ERF_prob.cpp @@ -71,13 +71,6 @@ Real compute_relative_humidity (const Real z, const Real height, const Real z_tr } } -AMREX_FORCE_INLINE -AMREX_GPU_HOST_DEVICE -Real compute_vapor_pressure (const Real p_s, const Real RH) -{ - return p_s*RH; -} - AMREX_FORCE_INLINE AMREX_GPU_HOST_DEVICE Real vapor_mixing_ratio (const Real z, const Real height, const Real p_b, const Real T_b, const Real RH) @@ -121,11 +114,13 @@ void Problem::compute_rho (const Real& z, const Real& pressure, Real& theta, Rea { theta = compute_theta(z); - T_b = compute_temperature(pressure, theta); + T_b = getTgivenPandTh(theta, pressure, (R_d/Cp_d)); Real RH = compute_relative_humidity(z, parms.height, parms.z_tr, pressure, T_b); q_v = vapor_mixing_ratio(z, parms.height, pressure, T_b, RH); - rho = pressure/(R_d*T_b*(1.0 + (R_v/R_d)*q_v)); + + rho = getRhogivenTandPress(T_b, pressure, q_v); rho = rho*(1.0 + q_v); + T_dp = compute_dewpoint_temperature(T_b, RH); } @@ -170,26 +165,19 @@ Problem::init_isentropic_hse_no_terrain(Real *theta, Real* r, Real* p, Real *q_v const Real& dz, const Real& prob_lo_z, const int& khi) { - - //FILE *file_IC; - //file_IC = fopen("input_sounding_probcpp.txt","w"); Real z, T_b, T_dp; // Compute the quantities at z = 0.5*dz (first cell center) z = prob_lo_z + 0.5*dz; p[0] = p_0; compute_p_k(p[0], p_0, theta[0], r[0], q_v[0], T_dp, T_b, dz, z, 0.0); - //fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[0], r[0], theta[0], q_v[0]); for (int k=1;k<=khi;k++){ z = prob_lo_z + (k+0.5)*dz; p[k] = p[k-1]; compute_p_k(p[k], p[k-1], theta[k], r[k], q_v[k], T_dp, T_b, dz, z, r[k-1]); - //fprintf(file_IC, "%0.15g %0.15g %0.15g %0.15g %0.15g %0.15g %0.15g\n", z, T_b-273.15, T_dp, p[k], r[k], theta[k], q_v[k]); } - //fclose(file_IC); - r[khi+1] = r[khi]; } @@ -199,8 +187,6 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, std::unique_ptr& z_phys_nd, Geometry const& geom) { - - const Real prob_lo_z = geom.ProbLo()[2]; const Real dz = geom.CellSize()[2]; const int khi = geom.Domain().bigEnd()[2]; @@ -241,17 +227,11 @@ Problem::erf_init_dens_hse_moist (MultiFab& rho_hse, Vector h_t(khi+2); Vector h_q_v(khi+2); - amrex::Gpu::DeviceVector d_r(khi+2); - amrex::Gpu::DeviceVector d_p(khi+2); - amrex::Gpu::DeviceVector d_t(khi+2); - amrex::Gpu::DeviceVector d_q_v(khi+2); - init_isentropic_hse_no_terrain(h_t.data(), h_r.data(),h_p.data(), h_q_v.data(), dz,prob_lo_z,khi); + amrex::Gpu::DeviceVector d_r(khi+2); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_r.begin(), h_r.end(), d_r.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_p.begin(), h_p.end(), d_p.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_t.begin(), h_t.end(), d_t.begin()); - amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, h_q_v.begin(), h_q_v.end(), d_q_v.begin()); Real* r = d_r.data(); @@ -355,18 +335,18 @@ Problem::init_custom_pert ( } theta_total = t[k] + delta_theta; - temperature = compute_temperature(p[k], theta_total); - Real T_b = compute_temperature(p[k], t[k]); + temperature = getTgivenPandTh(theta_total, p[k], (R_d/Cp_d)); + Real T_b = getTgivenPandTh(t[k] , p[k], (R_d/Cp_d)); RH = compute_relative_humidity(z, height, z_tr, p[k], T_b); Real q_v_hot = vapor_mixing_ratio(z, height, p[k], T_b, RH); rho = p[k]/(R_d*temperature*(1.0 + (R_v/R_d)*q_v_hot)); // Compute background quantities - Real temperature_back = compute_temperature(p[k], t[k]); - Real T_back = compute_temperature(p[k], t[k]); + Real temperature_back = getTgivenPandTh(t[k], p[k], (R_d/Cp_d)); + Real T_back = getTgivenPandTh(t[k], p[k], (R_d/Cp_d)); Real RH_back = compute_relative_humidity(z, height, z_tr, p[k], T_back); Real q_v_back = vapor_mixing_ratio(z, height, p[k], T_back, RH_back); - Real rho_back = p[k]/(R_d*temperature_back*(1.0 + (R_v/R_d)*q_v_back)); + Real rho_back = getRhogivenTandPress(temperature_back, p[k], q_v_back); // This version perturbs rho but not p state_pert(i, j, k, RhoTheta_comp) = rho*theta_total - rho_back*t[k]*(1.0 + (R_v/R_d)*q_v_back);// rho*d_t[k]*(1.0 + R_v_by_R_d*q_v_hot); diff --git a/Source/ERF_make_new_level.cpp b/Source/ERF_make_new_level.cpp index f9d29ce2b..dbe285ec8 100644 --- a/Source/ERF_make_new_level.cpp +++ b/Source/ERF_make_new_level.cpp @@ -63,11 +63,6 @@ void ERF::MakeNewLevelFromScratch (int lev, Real time, const BoxArray& ba_in, m_factory[lev] = std::make_unique(); #endif - // The number of ghost cells for density must be 1 greater than that for velocity - // so that we can go back in forth between velocity and momentum on all faces - // int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; - // int ngrow_vels = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff); - auto& lev_new = vars_new[lev]; auto& lev_old = vars_old[lev]; @@ -304,7 +299,6 @@ ERF::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMapp int ncomp_cons = vars_new[lev][Vars::cons].nComp(); IntVect ngrow_state = vars_new[lev][Vars::cons].nGrowVect(); - // int ngrow_state = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff) + 1; int ngrow_vels = ComputeGhostCells(solverChoice.advChoice, solverChoice.use_NumDiff); Vector temp_lev_new(Vars::NumTypes); diff --git a/Source/TimeIntegration/ERF_TI_fast_rhs_fun.H b/Source/TimeIntegration/ERF_TI_substep_fun.H similarity index 100% rename from Source/TimeIntegration/ERF_TI_fast_rhs_fun.H rename to Source/TimeIntegration/ERF_TI_substep_fun.H diff --git a/Source/TimeIntegration/ERF_advance_dycore.cpp b/Source/TimeIntegration/ERF_advance_dycore.cpp index 64313e133..e211ff464 100644 --- a/Source/TimeIntegration/ERF_advance_dycore.cpp +++ b/Source/TimeIntegration/ERF_advance_dycore.cpp @@ -276,8 +276,8 @@ void ERF::advance_dycore(int level, cons_to_prim(state_old[IntVars::cons], state_old[IntVars::cons].nGrow()); #include "ERF_TI_no_substep_fun.H" +#include "ERF_TI_substep_fun.H" #include "ERF_TI_slow_rhs_fun.H" -#include "ERF_TI_fast_rhs_fun.H" // *************************************************************************************** // Setup the integrator and integrate for a single timestep diff --git a/Source/TimeIntegration/Make.package b/Source/TimeIntegration/Make.package index d117f5e87..7333d3676 100644 --- a/Source/TimeIntegration/Make.package +++ b/Source/TimeIntegration/Make.package @@ -12,9 +12,9 @@ CEXE_sources += ERF_fast_rhs_N.cpp CEXE_sources += ERF_fast_rhs_T.cpp CEXE_sources += ERF_fast_rhs_MT.cpp -CEXE_headers += ERF_TI_fast_rhs_fun.H CEXE_headers += ERF_TI_slow_rhs_fun.H CEXE_headers += ERF_TI_no_substep_fun.H +CEXE_headers += ERF_TI_substep_fun.H CEXE_headers += ERF_TI_fast_headers.H CEXE_headers += ERF_TI_slow_headers.H CEXE_headers += ERF_TI_utils.H diff --git a/Source/Utils/ERF_EOS.H b/Source/Utils/ERF_EOS.H index 3d47639b8..dcae76733 100644 --- a/Source/Utils/ERF_EOS.H +++ b/Source/Utils/ERF_EOS.H @@ -9,34 +9,38 @@ /** * Function to return potential temperature given pressure and temperature * - * @params[in] pressure - * @params[in] temperature - * @params[in] rdOcp ratio of R_d to c_p + * @params[in ] P pressure + * @params[in ] T temperature + * @params[in ] rdOcp ratio of R_d to c_p + * @params[ out] potential temperature */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenPandT(const amrex::Real T, const amrex::Real P, const amrex::Real rdOcp) { - return T*std::pow(p_0/P, rdOcp); + return T * std::pow(p_0/P, rdOcp); } /** * Function to return temperature given pressure and potential temperature * - * @params[in] pressure - * @params[in] potential temperature - * @params[in] rdOcp ratio of R_d to c_p + * @params[in ] P pressure + * @params[in ] th potential temperature + * @params[in ] rdOcp ratio of R_d to c_p + * @params[ out] temperature */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getTgivenPandTh(const amrex::Real P, const amrex::Real Th, const amrex::Real rdOcp) +amrex::Real getTgivenPandTh(const amrex::Real th, const amrex::Real P, const amrex::Real rdOcp) { - return Th / std::pow(p_0/P, rdOcp); + return th / std::pow(p_0/P, rdOcp); } /** * Function to return temperature given density and potential temperature * - * @params[in] rho density - * @params[in] rhotheta density times potential temperature theta + * @params[in ] rho density + * @params[in ] rhotheta (density times potential temperature) + * @params[in ] qv water vapor + * @params[ out] temperature */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, const amrex::Real qv=0.0) @@ -53,15 +57,18 @@ amrex::Real getTgivenRandRTh(const amrex::Real rho, const amrex::Real rhotheta, /** * Function to return potential temperature given density and temperature * - * @params[in] rho density - * @params[in] T temperature - * @params[in] rdOcp ratio of R_d to c_p + * @params[in ] rho density + * @params[in ] T temperature + * @params[in ] rdOcp ratio of R_d to c_p + * @params[in ] qv water vapor + * @params[ out] potential temperature */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const amrex::Real rdOcp, const amrex::Real qv=0.0) { // p = rho_d * R_d * T_moist amrex::Real p_loc = rho * R_d * T * (1.0 + R_v/R_d*qv); + // theta_d = T * (p0/p)^(R_d/C_p) return T * std::pow((p_0/p_loc),rdOcp); } @@ -69,49 +76,69 @@ amrex::Real getThgivenRandT(const amrex::Real rho, const amrex::Real T, const am /** * Function to return pressure given density times theta * - * @params[in] rhotheta density times potential temperature - * @params[in] qv water vapor + * @params[in ] rhotheta density times potential temperature + * @params[in ] qv water vapor + * @params[ out] pressure */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getPgivenRTh(const amrex::Real rhotheta, const amrex::Real qv = 0.) { - return p_0 * std::pow(R_d * rhotheta * (1.0+(R_v/R_d)*qv) * ip_0, Gamma); + return p_0 * std::pow(R_d * rhotheta * ( amrex::Real(1.0) + (R_v/R_d)*qv) * ip_0, Gamma); } /** * Function to return density given theta and pressure * - * @params[in] theta potential temperature - * @params[in] p pressure - * @params[in] rdOcp ratio of R_d to c_p + * @params[in ] theta potential temperature + * @params[in ] p pressure + * @params[in ] rdOcp ratio of R_d to c_p + * @params[in ] qv water vapor + * @params[ out] density */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -amrex::Real getRhogivenThetaPress (const amrex::Real theta, const amrex::Real p, const amrex::Real rdOcp, const amrex::Real qv=0.0) +amrex::Real getRhogivenThetaPress (const amrex::Real th, const amrex::Real p, const amrex::Real rdOcp, const amrex::Real qv=0.0) { // We should be using moist value of theta when using moisture // theta_m = theta * (1 + R_v/R_d*qv) - return std::pow(p_0, rdOcp) * std::pow(p, iGamma) / (R_d * theta * (1.0 + R_v/R_d*qv) ); + return std::pow(p_0, rdOcp) * std::pow(p, iGamma) / (R_d * th * (amrex::Real(1.0) + R_v/R_d*qv) ); +} + +/** + * Function to return density given temperature and pressure + * + * @params[in ] theta potential temperature + * @params[in ] p pressure + * @params[in ] rdOcp ratio of R_d to c_p + * @params[in ] qv water vapor + * @params[ out] density +*/ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real getRhogivenTandPress (const amrex::Real T, const amrex::Real p, const amrex::Real qv=0.0) +{ + return p / ( R_d * T * (amrex::Real(1.0) + (R_v/R_d)*qv) ); } /** * Function to return dP/drho at constant theta * - * @params[in] theta potential temperature - * @params[in] p pressure - * @params[in] rdOcp ratio of R_d to c_p + * @params[in ] rho density + * @params[in ] theta potential temperature + * @params[in ] qv water vapor + * @params[ out] pressure */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getdPdRgivenConstantTheta(const amrex::Real rho, const amrex::Real theta, const amrex::Real qv=0.0) { // We should be using moist value of theta when using moisture // theta_m = theta * (1 + R_v/R_d*qv) - return Gamma * p_0 * std::pow( (R_d * theta * (1.0 + R_v/R_d*qv) * ip_0), Gamma) * std::pow(rho, Gamma-1.0) ; + return Gamma * p_0 * std::pow( (R_d * theta * (amrex::Real(1.0) + R_v/R_d*qv) * ip_0), Gamma) * std::pow(rho, Gamma-1.0) ; } /** * Function to return the Exner function pi given pressure - * @params[in] p pressure - * @params[in] rdOcp ratio of R_d to c_p + * @params[in ] p pressure + * @params[in ] rdOcp ratio of R_d to c_p + * @params[ out] Exner function */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getExnergivenP(const amrex::Real P, const amrex::Real rdOcp) @@ -121,10 +148,12 @@ amrex::Real getExnergivenP(const amrex::Real P, const amrex::Real rdOcp) } /** - * Function to return the Exner function pi given densith times potential temperature + * Function to return the Exner function pi given density times potential temperature * - * @params[in] rhotheta density times potential temperature - * @params[in] rdOcp ratio of R_d to c_p + * @params[in ] rhotheta density times potential temperature + * @params[in ] rdOcp ratio of R_d to c_p + * @params[in ] qv water vapor + * @params[ out] Exner function */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp, const amrex::Real qv=0.0 ) @@ -136,9 +165,11 @@ amrex::Real getExnergivenRTh(const amrex::Real rhotheta, const amrex::Real rdOcp } /** - * Function to return the density given pressure + * Function to return (rho theta) given pressure * - * @params[in] p pressure + * @params[in ] p pressure + * @params[in ] qv water vapor + * @params[ out] density times potential temperature */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=0.0) @@ -149,5 +180,11 @@ amrex::Real getRhoThetagivenP(const amrex::Real p, const amrex::Real qv=0.0) return std::pow(p*std::pow(p_0, Gamma-1), iGamma) * iR_d / (1.0 + R_v/R_d*qv) ; } +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real compute_vapor_pressure (const amrex::Real p_s, const amrex::Real RH) +{ + return p_s*RH; +} + #endif