From d35fda2c4ae0df657a1dc274d8c09ab3d98f5c56 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 13 Jun 2026 20:01:40 -0700 Subject: [PATCH 1/2] ERF_Constants.H --- Source/BoundaryConditions/ERF_SurfaceLayer.H | 17 ++++++++----- Source/EB/ERF_EBMOSTStress.H | 16 +++++++----- Source/ERF.H | 12 ++++----- Source/ERF.cpp | 26 ++++++++++---------- Source/IO/ERF_ReadFromWRFBdy.H | 8 +++--- Source/IO/ERF_ReadFromWRFBdy.cpp | 9 +++---- Source/IO/ERF_ReadFromWRFLow.cpp | 4 +-- Source/TimeIntegration/ERF_TimeStep.cpp | 2 +- 8 files changed, 51 insertions(+), 43 deletions(-) diff --git a/Source/BoundaryConditions/ERF_SurfaceLayer.H b/Source/BoundaryConditions/ERF_SurfaceLayer.H index b714cb8fd..db8a93e5b 100644 --- a/Source/BoundaryConditions/ERF_SurfaceLayer.H +++ b/Source/BoundaryConditions/ERF_SurfaceLayer.H @@ -439,23 +439,28 @@ public: // 2D MFs for U*, T*, T_surf //-------------------------------------------------------- +#ifdef AMREX_USE_FLOAT + amrex::Real bogus_val = amrex::Real(1.234e10); +#else + amrex::Real bogus_val = amrex::Real(1.234e30); +#endif u_star[lev] = std::make_unique(ba_flux, dm, ncomp, ng); - u_star[lev]->setVal(1.E34); + u_star[lev]->setVal(bogus_val); w_star[lev] = std::make_unique(ba_flux, dm, ncomp, ng); - w_star[lev]->setVal(1.E34); + w_star[lev]->setVal(bogus_val); t_star[lev] = std::make_unique(ba_flux, dm, ncomp, ng); - t_star[lev]->setVal(0.0); // default to neutral + t_star[lev]->setVal(zero); // default to neutral q_star[lev] = std::make_unique(ba_flux, dm, ncomp, ng); - q_star[lev]->setVal(0.0); // default to dry + q_star[lev]->setVal(zero); // default to dry olen[lev] = std::make_unique(ba_flux, dm, ncomp, ng); - olen[lev]->setVal(1.E34); + olen[lev]->setVal(bogus_val); pblh[lev] = std::make_unique(ba_flux, dm, ncomp, ng); - pblh[lev]->setVal(1.E34); + pblh[lev]->setVal(bogus_val); t_surf[lev] = std::make_unique(ba_flux, dm, ncomp, ng); t_surf[lev]->setVal(default_land_surf_temp); diff --git a/Source/EB/ERF_EBMOSTStress.H b/Source/EB/ERF_EBMOSTStress.H index d4fbab247..0c655b64a 100644 --- a/Source/EB/ERF_EBMOSTStress.H +++ b/Source/EB/ERF_EBMOSTStress.H @@ -337,7 +337,7 @@ struct moeng_flux_eb amrex::Real wsp_mean = umm_arr(i,j,0); wsp_mean = std::max(wsp_mean, WSMIN); - amrex::Real wsp = sqrt(velx*velx+vely*vely); + amrex::Real wsp = std::sqrt(velx*velx+vely*vely); amrex::Real num1 = wsp * (qv_mean-qv_surf); amrex::Real num2 = wsp_mean * (qv-qv_mean); @@ -407,7 +407,7 @@ struct moeng_flux_eb wsp_mean = std::max(wsp_mean, WSMIN); // Use tangential velocity magnitude instead of Cartesian - amrex::Real wsp = sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent); + amrex::Real wsp = std::sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent); amrex::Real num1 = wsp * (theta_mean-theta_surf); amrex::Real num2 = wsp_mean * (theta-theta_mean); @@ -613,7 +613,7 @@ struct moeng_flux_eb // multiplying the modeled shear stress (rho*ustar^2) with // a factor of umean/wsp_mean for directionality; this factor // modifies the denominator from what is in Moeng amrex::Real(1984.) - amrex::Real wsp = sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent); + amrex::Real wsp = std::sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent); amrex::Real num1 = wsp * umean; amrex::Real num2 = wsp_mean * (velx_tangent-umean); @@ -818,7 +818,7 @@ struct moeng_flux_eb // multiplying the modeled shear stress (rho*ustar^2) with // a factor of vmean/wsp_mean for directionality; this factor // modifies the denominator from what is in Moeng amrex::Real(1984.) - amrex::Real wsp = sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent); + amrex::Real wsp = std::sqrt(velx_tangent*velx_tangent+vely_tangent*vely_tangent); amrex::Real num1 = wsp * vmean; amrex::Real num2 = wsp_mean * (vely_tangent-vmean); @@ -829,7 +829,11 @@ struct moeng_flux_eb } private: - const amrex::Real eps = 1e-15; +#ifdef AMREX_USE_FLOAT + const amrex::Real eps = amrex::Real(1e-6); +#else + const amrex::Real eps = amrex::Real(1e-12); +#endif const amrex::Real WSMIN = amrex::Real(0.1); // minimum wind speed }; -#endif \ No newline at end of file +#endif diff --git a/Source/ERF.H b/Source/ERF.H index 076fdf8d7..497e65ec6 100644 --- a/Source/ERF.H +++ b/Source/ERF.H @@ -758,14 +758,14 @@ private: amrex::Vector z_vel_lateral, amrex::Vector T_lateral); - static amrex::Real start_bdy_time; - static amrex::Real final_bdy_time; + static double start_bdy_time; + static double final_bdy_time; - static amrex::Real start_low_time; - static amrex::Real final_low_time; + static double start_low_time; + static double final_low_time; - static amrex::Real bdy_time_interval; - static amrex::Real low_time_interval; + static double bdy_time_interval; + static double low_time_interval; // *** *** FArrayBox's for holding the SURFACE data // amrex::IArrayBox NC_IVGTYP_fab; // Vegetation type (IVGTYP); Discrete numbers; diff --git a/Source/ERF.cpp b/Source/ERF.cpp index 182b43ce0..0eba08d64 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -46,14 +46,14 @@ Real ERF::start_time = zero; Real ERF::stop_time = std::numeric_limits::max(); #ifdef ERF_USE_NETCDF -Real ERF::start_bdy_time = zero; -Real ERF::final_bdy_time = -one; +double ERF::start_bdy_time = zero; +double ERF::final_bdy_time = -one; -Real ERF::start_low_time = zero; -Real ERF::final_low_time = -one; +double ERF::start_low_time = zero; +double ERF::final_low_time = -one; -Real ERF::bdy_time_interval = std::numeric_limits::max(); -Real ERF::low_time_interval = std::numeric_limits::max(); +double ERF::bdy_time_interval = std::numeric_limits::max(); +double ERF::low_time_interval = std::numeric_limits::max(); #endif // Time step control @@ -871,7 +871,7 @@ ERF::post_timestep (int nstep, double time, Real dt_lev0) if ( rad_datalog_int > 0 && (((nstep+1) % rad_datalog_int == 0) || (nstep==0)) ) { if (rad[0]->hasDatalog()) { - rad[0]->WriteDataLog(time+start_time); + rad[0]->WriteDataLog(static_cast(time+start_time)); } } } @@ -1125,7 +1125,7 @@ ERF::InitData_post () bdy_data_xlo,bdy_data_xhi,bdy_data_ylo,bdy_data_yhi, start_bdy_time, final_bdy_time); - Real time_since_start_bdy = t_new[0] + start_time - start_bdy_time; + double time_since_start_bdy = t_new[0] + start_time - start_bdy_time; int n_time_old = static_cast(time_since_start_bdy / bdy_time_interval); MultiFab r_hse(base_state[0], make_alias, BaseState::r0_comp, 1); Array area_vec = {ax[0].get(), ay[0].get(), az[0].get()}; @@ -1164,7 +1164,7 @@ ERF::InitData_post () sst_lev[lev].resize(low_data_zlo.size()); tsk_lev[lev].resize(low_data_zlo.size()); - Real time_since_start_low = t_new[0] + start_time - start_low_time; + double time_since_start_low = t_new[0] + start_time - start_low_time; int n_time_old = static_cast(time_since_start_low / low_time_interval); int ntimes = std::min(n_time_old+3, static_cast(low_data_zlo.size())); @@ -2003,7 +2003,7 @@ ERF::Interp2DArrays (int lev, const BoxArray& my_ba2d, const DistributionMapping sst_lev[lev].resize(sst_lev[lev-1].size()); } #ifdef ERF_USE_NETCDF - Real time_since_start_low = t_new[0] + start_time - start_low_time; + double time_since_start_low = t_new[0] + start_time - start_low_time; int n_time_old = static_cast(time_since_start_low / low_time_interval); int ntimes_to_interp = std::min(n_time_old+3, static_cast(sst_lev[lev-1].size())); #else @@ -2030,7 +2030,7 @@ ERF::Interp2DArrays (int lev, const BoxArray& my_ba2d, const DistributionMapping tsk_lev[lev].resize(tsk_lev[lev-1].size()); } #ifdef ERF_USE_NETCDF - Real time_since_start_low = t_new[0] + start_time - start_low_time; + double time_since_start_low = t_new[0] + start_time - start_low_time; int n_time_old = static_cast(time_since_start_low / low_time_interval); int ntimes_to_interp = std::min(n_time_old+3, static_cast(tsk_lev[lev-1].size())); #else @@ -2107,7 +2107,7 @@ ERF::Interp2DArrays (int lev, const BoxArray& my_ba2d, const DistributionMapping if (sst_lev[lev][0]) { // Call FillPatchTwoLevels which ASSUMES that all ghost cells at lev-1 have already been filled #ifdef ERF_USE_NETCDF - Real time_since_start_low = t_new[0] + start_time - start_low_time; + double time_since_start_low = t_new[0] + start_time - start_low_time; int n_time_old = static_cast(time_since_start_low / low_time_interval); int ntimes_to_interp = std::min(n_time_old+3, static_cast(sst_lev[lev-1].size())); #else @@ -2131,7 +2131,7 @@ ERF::Interp2DArrays (int lev, const BoxArray& my_ba2d, const DistributionMapping if (tsk_lev[lev][0]) { // Call FillPatchTwoLevels which ASSUMES that all ghost cells at lev-1 have already been filled #ifdef ERF_USE_NETCDF - Real time_since_start_low = t_new[0] + start_time - start_low_time; + double time_since_start_low = t_new[0] + start_time - start_low_time; int n_time_old = static_cast(time_since_start_low / low_time_interval); int ntimes_to_interp = std::min(n_time_old+3, static_cast(tsk_lev[lev-1].size())); #else diff --git a/Source/IO/ERF_ReadFromWRFBdy.H b/Source/IO/ERF_ReadFromWRFBdy.H index 0b0239e7d..42d116404 100644 --- a/Source/IO/ERF_ReadFromWRFBdy.H +++ b/Source/IO/ERF_ReadFromWRFBdy.H @@ -12,13 +12,13 @@ #ifdef ERF_USE_NETCDF -amrex::Real +double read_times_from_wrfbdy (const std::string& nc_bdy_file, amrex::Vector>& bdy_data_xlo, amrex::Vector>& bdy_data_xhi, amrex::Vector>& bdy_data_ylo, amrex::Vector>& bdy_data_yhi, - amrex::Real& start_bdy_time, amrex::Real& final_bdy_time); + double& start_bdy_time, double& final_bdy_time); void read_and_convert_from_wrfbdy (const int itime, const std::string& nc_bdy_file, @@ -33,13 +33,13 @@ read_and_convert_from_wrfbdy (const int itime, const std::string& nc_bdy_file, amrex::Array& area_vec, const amrex::Geometry& geom, const bool& use_moist, const amrex::Vector& domain_bcs_type_h, - int real_width, amrex::Real bdy_time_interval, + int real_width, double bdy_time_interval, bool is_anelastic, bool do_conversion = true); amrex::Real read_times_from_wrflow (const std::string& nc_low_file, amrex::Vector>& low_data_zlo, - amrex::Real& start_low_time, amrex::Real& final_low_time); + double& start_low_time, double& final_low_time); void read_from_wrflow (const int itime, diff --git a/Source/IO/ERF_ReadFromWRFBdy.cpp b/Source/IO/ERF_ReadFromWRFBdy.cpp index 14840cf50..5a66a73da 100644 --- a/Source/IO/ERF_ReadFromWRFBdy.cpp +++ b/Source/IO/ERF_ReadFromWRFBdy.cpp @@ -26,14 +26,14 @@ namespace WRFBdyTypes { }; } -Real +double read_times_from_wrfbdy (const std::string& nc_bdy_file, Vector>& bdy_data_xlo, Vector>& bdy_data_xhi, Vector>& bdy_data_ylo, Vector>& bdy_data_yhi, - Real& start_bdy_time, - Real& final_bdy_time) + double& start_bdy_time, + double& final_bdy_time) { Print() << "Loading boundary data from NetCDF file " << std::endl; @@ -163,7 +163,6 @@ convert_wrfbdy_data (const int itime, int ihi = domain.bigEnd()[0]; int jlo = domain.smallEnd()[1]; int jhi = domain.bigEnd()[1]; - int klo = domain.smallEnd()[2]; int khi = domain.bigEnd()[2]; // PH bounds limiting @@ -468,7 +467,7 @@ read_and_convert_from_wrfbdy (const int itime, const std::string& nc_bdy_file, const Geometry& geom, const bool& use_moist, const Vector& domain_bcs_type_h, - int real_width, Real bdy_time_interval, + int real_width, double bdy_time_interval, bool is_anelastic, bool do_conversion) { int ioproc = ParallelDescriptor::IOProcessorNumber(); // I/O rank diff --git a/Source/IO/ERF_ReadFromWRFLow.cpp b/Source/IO/ERF_ReadFromWRFLow.cpp index c195f66c7..87dfdb25b 100644 --- a/Source/IO/ERF_ReadFromWRFLow.cpp +++ b/Source/IO/ERF_ReadFromWRFLow.cpp @@ -17,10 +17,10 @@ using namespace amrex; #ifdef ERF_USE_NETCDF -Real +double read_times_from_wrflow (const std::string& nc_low_file, Vector>& low_data_zlo, - Real& start_low_time, Real& final_low_time) + double& start_low_time, double& final_low_time) { Print() << "Loading lower boundary data from NetCDF file " << std::endl; diff --git a/Source/TimeIntegration/ERF_TimeStep.cpp b/Source/TimeIntegration/ERF_TimeStep.cpp index 9fa8727f4..a358e3ac6 100644 --- a/Source/TimeIntegration/ERF_TimeStep.cpp +++ b/Source/TimeIntegration/ERF_TimeStep.cpp @@ -37,7 +37,7 @@ ERF::timeStep (int lev, double time, int /*iteration*/) Array area_vec = {ax[lev].get(), ay[lev].get(), az[lev].get()}; int ntimes = bdy_data_xlo.size(); - Real time_since_start_bdy = time + start_time - start_bdy_time; + double time_since_start_bdy = time + start_time - start_bdy_time; int n_time_old = std::min(static_cast( (time_since_start_bdy ) / bdy_time_interval), ntimes-1); int n_time_new = std::min(static_cast( (time_since_start_bdy+dt[lev]) / bdy_time_interval), ntimes-1); From 2279f8342126e627a24f47b36e959d18a761b202 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Sat, 13 Jun 2026 20:02:13 -0700 Subject: [PATCH 2/2] Make more of the time variables double --- Exec/ERF_Prob.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Exec/ERF_Prob.cpp b/Exec/ERF_Prob.cpp index 4a6178667..67a46591d 100644 --- a/Exec/ERF_Prob.cpp +++ b/Exec/ERF_Prob.cpp @@ -263,7 +263,7 @@ Problem::update_rhotheta_sources (const Real& time, { if (src->empty()) return; - const int khi = geom.Domain().bigEnd()[2]; + const int khi = geom.Domain().bigEnd()[2]; // If the z coordinate varies in time and or space, then the the height // needs to be calculated at each time step. Here, we assume that only