From 9ad5c827de38f8b827112a48029dd4ec2b1ec712 Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Thu, 22 Jan 2026 10:47:16 +0000 Subject: [PATCH 1/5] Experimental set up * Adds profiling tags, standardise seed and set num threads to 1 in issl_rhi140 example Co-authored-by: Coco Yeung --- Code.v05-00/CMakeLists.txt | 7 ++++++- examples/issl_rhi140/input.yaml | 4 ++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/Code.v05-00/CMakeLists.txt b/Code.v05-00/CMakeLists.txt index 895dad2a..11d0ee0e 100644 --- a/Code.v05-00/CMakeLists.txt +++ b/Code.v05-00/CMakeLists.txt @@ -11,6 +11,7 @@ option(DEBUG "Enable debug mode" OFF) SET(RINGS 0) SET(OMP 1) option(BUILD_TEST ON) +option(ENABLE_PROFILING "Enable perf-friendly build (-g -O2)" OFF) if (CMAKE_BUILD_TYPE STREQUAL "Debug") include(CheckCXXCompilerFlag) @@ -35,6 +36,10 @@ elseif (NOT CMAKE_BUILD_TYPE OR CMAKE_BUILD_TYPE STREQUAL "") #-Michael #set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ffast-math -fno-trapping-math -funsafe-math-optimizations") endif() +if(ENABLE_PROFILING) + message(STATUS "Profiling enabled: adding debug symbols and frame pointers") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O2 -fno-omit-frame-pointer") +endif() set(CMAKE_EXPORT_COMPILE_COMMANDS ON) execute_process(COMMAND git rev-parse --short HEAD @@ -132,4 +137,4 @@ add_subdirectory(tests) if (DEBUG) message(STATUS "DEBUG mode is enabled") -endif() +endif() \ No newline at end of file diff --git a/examples/issl_rhi140/input.yaml b/examples/issl_rhi140/input.yaml index 686704eb..b874d3ed 100644 --- a/examples/issl_rhi140/input.yaml +++ b/examples/issl_rhi140/input.yaml @@ -2,7 +2,7 @@ SIMULATION MENU: # Only one of parameter sweep or MC simulation can be on. # Parameter sweep lets you specify an arbitrary number of custom values for each parameter; Monte Carlo simulation is self-explanatory. # At the moment, you cannot mix and match MC sim and param sweep on each individual parameter. - OpenMP Num Threads (positive int): 8 + OpenMP Num Threads (positive int): 1 PARAM SWEEP SUBMENU: Parameter sweep (T/F): T #-OR--------------- @@ -32,7 +32,7 @@ SIMULATION MENU: Run box model (T/F): F netCDF filename format (string): APCEMM_BOX_CASE_* RANDOM NUMBER GENERATION SUBMENU: - Force seed value (T/F): F + Force seed value (T/F): T Seed value (positive int): 0 EPM type (original/external/new): original From 8f5cd64302f6b066d4b7fb2007ed7e71abdfff1e Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Thu, 22 Jan 2026 12:43:57 +0000 Subject: [PATCH 2/5] Add input files for different weather conditions Add input files for different weather conditions --------- Co-authored-by: Coco Yeung --- examples/issl_rhi140/default.yaml | 216 ++++++++++++++++++++++++++ examples/issl_rhi140/rhi.yaml | 216 ++++++++++++++++++++++++++ examples/issl_rhi140/shear.yaml | 216 ++++++++++++++++++++++++++ examples/issl_rhi140/temperature.yaml | 216 ++++++++++++++++++++++++++ examples/issl_rhi140/w.yaml | 216 ++++++++++++++++++++++++++ 5 files changed, 1080 insertions(+) create mode 100644 examples/issl_rhi140/default.yaml create mode 100644 examples/issl_rhi140/rhi.yaml create mode 100644 examples/issl_rhi140/shear.yaml create mode 100644 examples/issl_rhi140/temperature.yaml create mode 100644 examples/issl_rhi140/w.yaml diff --git a/examples/issl_rhi140/default.yaml b/examples/issl_rhi140/default.yaml new file mode 100644 index 00000000..8eafca27 --- /dev/null +++ b/examples/issl_rhi140/default.yaml @@ -0,0 +1,216 @@ +SIMULATION MENU: + # Only one of parameter sweep or MC simulation can be on. + # Parameter sweep lets you specify an arbitrary number of custom values for each parameter; Monte Carlo simulation is self-explanatory. + # At the moment, you cannot mix and match MC sim and param sweep on each individual parameter. + OpenMP Num Threads (positive int): 1 + PARAM SWEEP SUBMENU: + Parameter sweep (T/F): T + #-OR--------------- + Run Monte Carlo (T/F): F + Num Monte Carlo runs (int): 2 + # Where APCEMM output for this set of runs will go + OUTPUT SUBMENU: + Output folder (string): APCEMM_out/ + Overwrite if folder exists (T/F): T + # FFT options (for spectral solver) + Use threaded FFT (T/F): F + FFTW WISDOM SUBMENU: + Use FFTW WISDOM (T/F): T + Dir w/ write permission (string): ./ + # This mostly contains information on background aerosol concentration; Not too relevant for contrail behavior + Input background condition (string): ../../input_data/init.txt + # All parameters here are overwritten in EMISSION INDICES SUBMENU + Input engine emissions (string): ../../input_data/ENG_EI.txt + # Ignore/Don't change these, these are deprecated features. + SAVE FORWARD RESULTS SUBMENU: + Save forward results (T/F): F + netCDF filename format (string): APCEMM_Case_* + ADJOINT OPTIMIZATION SUBMENU: + Turn on adjoint optim. (T/F): F + netCDF filename format (string): APCEMM_ADJ_Case_* + BOX MODEL SUBMENU: + Run box model (T/F): F + netCDF filename format (string): APCEMM_BOX_CASE_* + RANDOM NUMBER GENERATION SUBMENU: + Force seed value (T/F): T + Seed value (positive int): 0 + EPM type (original/external/new): original + +# Format of parameter items: +# Param name [unit] (Variable type) +PARAMETER MENU: + # Parameter sweep format : Format is either: x1 x2 x3 or start:increment:end + # : Example: 200 220 240 and 200:20:240 are identical + + # Monte Carlo simulation : min:max + # : Example: 200:240 will generate values for the parameter in between 200 and 240 + + # Maximum simulation time if contrail isn't gone by then: + Plume Process [hr] (double): 12 + + # Temperature, RH, and wind shear can be overwritten if using meteorological input files + METEOROLOGICAL PARAMETERS SUBMENU: + # Pressure altitude at which the contrail is initialized + Pressure [hPa] (double): 265 + Horiz. diff. coeff. [m^2/s] (double): 15.0 + # Can be overwritten if met file is passed + Verti. diff. [m^2/s] (double): 0.15 + Brunt-Vaisala Frequency [s^-1] (double): 0.013 + LOCATION AND TIME SUBMENU: + # Affects for solar zenith angle and photolysis which affect EPM/diurnal variations of temperature + LON [deg] (double): -15 + LAT [deg] (double): 60 + Emission day [1-365] (int): 81 + Emission time [hr] (double) : 8 + BACKGROUND MIXING RATIOS SUBMENU: + # Affects the EPM only + NOx [ppt] (double): 5100 + HNO3 [ppt] (double): 81.5 + O3 [ppb] (double): 100 + CO [ppb] (double): 40 + CH4 [ppm] (double): 1.76 + SO2 [ppt] (double): 7.25 + EMISSION INDICES SUBMENU: + # Engine specific parameters + # Default parameters are estimates for a 737-800 at 35k ft + # Affects EPM mostly: + NOx [g(NO2)/kg_fuel] (double): 10 # EDB for CFM56-5B3 + CO [g/kg_fuel] (double): 1 + UHC [g/kg_fuel] (double): 0.6 + # Affects contrail diffusion model: + SO2 [g/kg_fuel] (double): 1.2 # Assuming 600ppm of Sulfur + SO2 to SO4 conv [%] (double): 2 # AEDT paper, Barrett et al. (2010) + Soot [g/kg_fuel] (double): 0.008 # Assuming EI(nvPM) = 10^14 + Soot Radius [m] (double): 20.0E-09 + # Total fuel flow is then divided per engine + Total fuel flow [kg/s] (double) : 0.7 + Aircraft mass [kg] (double): 1.00e+05 # At beginning of cruise (35k ft), NPSS + Flight speed [m/s] (double): 265.42 # NPSS + Num. of engines [2/4] (int): 2 + Wingspan [m] (double): 34.32 # https://www.skybrary.aero/aircraft/b738 + # Values depend on where we assume condensation/freezing to occur and engine type + Core exit temp. [K] (double): 553.65 # NPSS + Exit bypass area [m^2] (double): 0.9772 + +TRANSPORT MENU: + # Keep on + Turn on Transport (T/F): T + # Outdated, not used (was used by spectral solver) + Fill Negative Values (T/F): T + Transport Timestep [min] (double): 1 + # Keep off: not sure of the effect yet + met updraft is included (if met file input) + PLUME UPDRAFT SUBMENU: + Turn on plume updraft (T/F): F + Updraft timescale [s] (double): 3600 + Updraft veloc. [cm/s] (double): 5 + +# Chemistry component of APCEMM hasn't been touched in a long time; leave off if only interested in contrail simulation +CHEMISTRY MENU: + Turn on Chemistry (T/F): F + Perform hetero. chem. (T/F): F + Chemistry Timestep [min] (double): 10 + Photolysis rates folder (string): /path/to/input/ + +AEROSOL MENU: + # Keep on + Turn on grav. settling (T/F): T + # Keep on + Turn on solid coagulation (T/F): T + # Keep off + Turn on liquid coagulation (T/F): F + Coag. timestep [min] (double): 60 + # Keep on + Turn on ice growth (T/F): T + Ice growth timestep [min] (double): 1 + +# At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true +# Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate + +# If using met. input: +# Exactly one of "Init temp. from met." and "Impose lapse rate" must be true +# Exactly one of "Init RH from met." and "Impose moist layer depth" must be true +METEOROLOGY MENU: + # --- MET INPUT OPTIONS ---+ + METEOROLOGICAL INPUT SUBMENU: + Use met. input (T/F): T + Met input file path (string): example_met_file.nc + # Frequency of met data availability + Time series data timestep [hr] (double): 1.0 + # If off, uses the parameter specified in METEOROLOGICAL PARAMETERS SUBMENU + Init temp. from met. (T/F): T + Temp. time series input (T/F): T + # Always interpolates in time, this controls spatial interpolation + Interpolate temp. met. data (T/F): T + # Same options as the ones for Temperature: + Init RH from met. (T/F): T + RH time series input (T/F): T + Interpolate RH met. data (T/F): T + Init wind shear from met. (T/F): T + Wind shear time series input (T/F): T + Interpolate shear met. data (T/F): F + Init vert. veloc. from met. data (T/F): T + Vert. veloc. time series input (T/F): F + Interpolate vert. veloc. met. data (T/F): F + # Option to modify NWP RH data + + # Perturbations to the temperature field can be added to the contrail simulation to account for some effects of + # atmospheric turbulence and gravity waves. + # Every N minutes, a temperature perturbation of e1*e2*T_amp where e1 and e2 are random variables + # uniformly distributed between -1 and 1 generated individually for each grid cell. + # The importance of turbulence and grav. waves can be increased by increasing T_amp + # The relative importance of grav. waves vs. turb. is increased by increasing the timescale + # See Lewellen, Persistent Contrails and Contrail Cirrus (2014) for full details. + TEMPERATURE PERTURBATION SUBMENU: + Enable Temp. Pert. (T/F): F + Temp. Perturb. Amplitude (double): 1.0 + Temp. Perturb. Timescale (min): 10 + +# The only thing that should be changed here is the save frequency +DIAGNOSTIC MENU: + netCDF filename format (string): trac_avg.apcemm.hhmm + # Leave "save species timeseries" off. It will do nothing without also turning chemistry on. + SPECIES TIMESERIES SUBMENU: + Save species timeseries (T/F): F + Inst timeseries file (string): ts_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Species indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Index 1 (ice) is the only relevant thing to save + AEROSOL TIMESERIES SUBMENU: + Save aerosol timeseries (T/F): T + Inst timeseries file (string): ts_aerosol_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Aerosol indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Keep off if chemistry is also off + PRODUCTION & LOSS SUBMENU: + Turn on P/L diag (T/F): F + Save O3 P/L (T/F): F + +# Sometimes you have to change YLIM_DOWN here if the supersaturated layer is very thick +# because YLIM_DOWN must be larger than the layer thickness. +ADVANCED OPTIONS MENU: + # Domain is defined as: X [-XLIM_LEFT, XLIM_RIGHT], Y [-YLIM_DOWN, YLIM_UP] + GRID SUBMENU: + NX (positive int): 200 + NY (positive int): 180 + XLIM_RIGHT (positive double): 1.0e+3 + XLIM_LEFT (positive double): 1.0e+3 + YLIM_UP (positive double): 300 + YLIM_DOWN (positive double): 1.5e+3 + INITIAL CONTRAIL SIZE SUBMENU: + #Depth = BaseDepth + DepthScalingFactor * Default_Depth + #Same formula for width + Base Contrail Depth [m] (double): 0.0 + Contrail Depth Scaling Factor [-] (double): 1.0 + Base Contrail Width [m] (double): 0.0 + Contrail Width Scaling Factor [-] (double): 1.0 + Ambient Lapse Rate [K/km] (double): -3.0 + Tropopause Pressure [Pa] (double): 2.0e+4 + EARLY PLUME SUBMENU: + Reference ice crystal count [#/m] (double): 3.38e12 + Reference wingspan [m] (double): 60.3 + # The values below can be overriden if required. + Override post-jet ice crystal count (T/F): F + Post-jet ice crystal count [#/m] (double): 1e15 + diff --git a/examples/issl_rhi140/rhi.yaml b/examples/issl_rhi140/rhi.yaml new file mode 100644 index 00000000..98e816b6 --- /dev/null +++ b/examples/issl_rhi140/rhi.yaml @@ -0,0 +1,216 @@ +SIMULATION MENU: + # Only one of parameter sweep or MC simulation can be on. + # Parameter sweep lets you specify an arbitrary number of custom values for each parameter; Monte Carlo simulation is self-explanatory. + # At the moment, you cannot mix and match MC sim and param sweep on each individual parameter. + OpenMP Num Threads (positive int): 1 + PARAM SWEEP SUBMENU: + Parameter sweep (T/F): T + #-OR--------------- + Run Monte Carlo (T/F): F + Num Monte Carlo runs (int): 2 + # Where APCEMM output for this set of runs will go + OUTPUT SUBMENU: + Output folder (string): APCEMM_out/ + Overwrite if folder exists (T/F): T + # FFT options (for spectral solver) + Use threaded FFT (T/F): F + FFTW WISDOM SUBMENU: + Use FFTW WISDOM (T/F): T + Dir w/ write permission (string): ./ + # This mostly contains information on background aerosol concentration; Not too relevant for contrail behavior + Input background condition (string): ../../input_data/init.txt + # All parameters here are overwritten in EMISSION INDICES SUBMENU + Input engine emissions (string): ../../input_data/ENG_EI.txt + # Ignore/Don't change these, these are deprecated features. + SAVE FORWARD RESULTS SUBMENU: + Save forward results (T/F): F + netCDF filename format (string): APCEMM_Case_* + ADJOINT OPTIMIZATION SUBMENU: + Turn on adjoint optim. (T/F): F + netCDF filename format (string): APCEMM_ADJ_Case_* + BOX MODEL SUBMENU: + Run box model (T/F): F + netCDF filename format (string): APCEMM_BOX_CASE_* + RANDOM NUMBER GENERATION SUBMENU: + Force seed value (T/F): T + Seed value (positive int): 0 + EPM type (original/external/new): original + +# Format of parameter items: +# Param name [unit] (Variable type) +PARAMETER MENU: + # Parameter sweep format : Format is either: x1 x2 x3 or start:increment:end + # : Example: 200 220 240 and 200:20:240 are identical + + # Monte Carlo simulation : min:max + # : Example: 200:240 will generate values for the parameter in between 200 and 240 + + # Maximum simulation time if contrail isn't gone by then: + Plume Process [hr] (double): 12 + + # Temperature, RH, and wind shear can be overwritten if using meteorological input files + METEOROLOGICAL PARAMETERS SUBMENU: + # Pressure altitude at which the contrail is initialized + Pressure [hPa] (double): 265 + Horiz. diff. coeff. [m^2/s] (double): 15.0 + # Can be overwritten if met file is passed + Verti. diff. [m^2/s] (double): 0.15 + Brunt-Vaisala Frequency [s^-1] (double): 0.013 + LOCATION AND TIME SUBMENU: + # Affects for solar zenith angle and photolysis which affect EPM/diurnal variations of temperature + LON [deg] (double): -15 + LAT [deg] (double): 60 + Emission day [1-365] (int): 81 + Emission time [hr] (double) : 8 + BACKGROUND MIXING RATIOS SUBMENU: + # Affects the EPM only + NOx [ppt] (double): 5100 + HNO3 [ppt] (double): 81.5 + O3 [ppb] (double): 100 + CO [ppb] (double): 40 + CH4 [ppm] (double): 1.76 + SO2 [ppt] (double): 7.25 + EMISSION INDICES SUBMENU: + # Engine specific parameters + # Default parameters are estimates for a 737-800 at 35k ft + # Affects EPM mostly: + NOx [g(NO2)/kg_fuel] (double): 10 # EDB for CFM56-5B3 + CO [g/kg_fuel] (double): 1 + UHC [g/kg_fuel] (double): 0.6 + # Affects contrail diffusion model: + SO2 [g/kg_fuel] (double): 1.2 # Assuming 600ppm of Sulfur + SO2 to SO4 conv [%] (double): 2 # AEDT paper, Barrett et al. (2010) + Soot [g/kg_fuel] (double): 0.008 # Assuming EI(nvPM) = 10^14 + Soot Radius [m] (double): 20.0E-09 + # Total fuel flow is then divided per engine + Total fuel flow [kg/s] (double) : 0.7 + Aircraft mass [kg] (double): 1.00e+05 # At beginning of cruise (35k ft), NPSS + Flight speed [m/s] (double): 265.42 # NPSS + Num. of engines [2/4] (int): 2 + Wingspan [m] (double): 34.32 # https://www.skybrary.aero/aircraft/b738 + # Values depend on where we assume condensation/freezing to occur and engine type + Core exit temp. [K] (double): 553.65 # NPSS + Exit bypass area [m^2] (double): 0.9772 + +TRANSPORT MENU: + # Keep on + Turn on Transport (T/F): T + # Outdated, not used (was used by spectral solver) + Fill Negative Values (T/F): T + Transport Timestep [min] (double): 1 + # Keep off: not sure of the effect yet + met updraft is included (if met file input) + PLUME UPDRAFT SUBMENU: + Turn on plume updraft (T/F): F + Updraft timescale [s] (double): 3600 + Updraft veloc. [cm/s] (double): 5 + +# Chemistry component of APCEMM hasn't been touched in a long time; leave off if only interested in contrail simulation +CHEMISTRY MENU: + Turn on Chemistry (T/F): F + Perform hetero. chem. (T/F): F + Chemistry Timestep [min] (double): 10 + Photolysis rates folder (string): /path/to/input/ + +AEROSOL MENU: + # Keep on + Turn on grav. settling (T/F): T + # Keep on + Turn on solid coagulation (T/F): T + # Keep off + Turn on liquid coagulation (T/F): F + Coag. timestep [min] (double): 60 + # Keep on + Turn on ice growth (T/F): T + Ice growth timestep [min] (double): 1 + +# At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true +# Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate + +# If using met. input: +# Exactly one of "Init temp. from met." and "Impose lapse rate" must be true +# Exactly one of "Init RH from met." and "Impose moist layer depth" must be true +METEOROLOGY MENU: + # --- MET INPUT OPTIONS ---+ + METEOROLOGICAL INPUT SUBMENU: + Use met. input (T/F): T + Met input file path (string): example_met_file_rhi.nc + # Frequency of met data availability + Time series data timestep [hr] (double): 1.0 + # If off, uses the parameter specified in METEOROLOGICAL PARAMETERS SUBMENU + Init temp. from met. (T/F): T + Temp. time series input (T/F): T + # Always interpolates in time, this controls spatial interpolation + Interpolate temp. met. data (T/F): T + # Same options as the ones for Temperature: + Init RH from met. (T/F): T + RH time series input (T/F): T + Interpolate RH met. data (T/F): T + Init wind shear from met. (T/F): T + Wind shear time series input (T/F): T + Interpolate shear met. data (T/F): F + Init vert. veloc. from met. data (T/F): T + Vert. veloc. time series input (T/F): F + Interpolate vert. veloc. met. data (T/F): F + # Option to modify NWP RH data + + # Perturbations to the temperature field can be added to the contrail simulation to account for some effects of + # atmospheric turbulence and gravity waves. + # Every N minutes, a temperature perturbation of e1*e2*T_amp where e1 and e2 are random variables + # uniformly distributed between -1 and 1 generated individually for each grid cell. + # The importance of turbulence and grav. waves can be increased by increasing T_amp + # The relative importance of grav. waves vs. turb. is increased by increasing the timescale + # See Lewellen, Persistent Contrails and Contrail Cirrus (2014) for full details. + TEMPERATURE PERTURBATION SUBMENU: + Enable Temp. Pert. (T/F): F + Temp. Perturb. Amplitude (double): 1.0 + Temp. Perturb. Timescale (min): 10 + +# The only thing that should be changed here is the save frequency +DIAGNOSTIC MENU: + netCDF filename format (string): trac_avg.apcemm.hhmm + # Leave "save species timeseries" off. It will do nothing without also turning chemistry on. + SPECIES TIMESERIES SUBMENU: + Save species timeseries (T/F): F + Inst timeseries file (string): ts_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Species indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Index 1 (ice) is the only relevant thing to save + AEROSOL TIMESERIES SUBMENU: + Save aerosol timeseries (T/F): T + Inst timeseries file (string): ts_aerosol_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Aerosol indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Keep off if chemistry is also off + PRODUCTION & LOSS SUBMENU: + Turn on P/L diag (T/F): F + Save O3 P/L (T/F): F + +# Sometimes you have to change YLIM_DOWN here if the supersaturated layer is very thick +# because YLIM_DOWN must be larger than the layer thickness. +ADVANCED OPTIONS MENU: + # Domain is defined as: X [-XLIM_LEFT, XLIM_RIGHT], Y [-YLIM_DOWN, YLIM_UP] + GRID SUBMENU: + NX (positive int): 200 + NY (positive int): 180 + XLIM_RIGHT (positive double): 1.0e+3 + XLIM_LEFT (positive double): 1.0e+3 + YLIM_UP (positive double): 300 + YLIM_DOWN (positive double): 1.5e+3 + INITIAL CONTRAIL SIZE SUBMENU: + #Depth = BaseDepth + DepthScalingFactor * Default_Depth + #Same formula for width + Base Contrail Depth [m] (double): 0.0 + Contrail Depth Scaling Factor [-] (double): 1.0 + Base Contrail Width [m] (double): 0.0 + Contrail Width Scaling Factor [-] (double): 1.0 + Ambient Lapse Rate [K/km] (double): -3.0 + Tropopause Pressure [Pa] (double): 2.0e+4 + EARLY PLUME SUBMENU: + Reference ice crystal count [#/m] (double): 3.38e12 + Reference wingspan [m] (double): 60.3 + # The values below can be overriden if required. + Override post-jet ice crystal count (T/F): F + Post-jet ice crystal count [#/m] (double): 1e15 + diff --git a/examples/issl_rhi140/shear.yaml b/examples/issl_rhi140/shear.yaml new file mode 100644 index 00000000..4ffdb0e3 --- /dev/null +++ b/examples/issl_rhi140/shear.yaml @@ -0,0 +1,216 @@ +SIMULATION MENU: + # Only one of parameter sweep or MC simulation can be on. + # Parameter sweep lets you specify an arbitrary number of custom values for each parameter; Monte Carlo simulation is self-explanatory. + # At the moment, you cannot mix and match MC sim and param sweep on each individual parameter. + OpenMP Num Threads (positive int): 1 + PARAM SWEEP SUBMENU: + Parameter sweep (T/F): T + #-OR--------------- + Run Monte Carlo (T/F): F + Num Monte Carlo runs (int): 2 + # Where APCEMM output for this set of runs will go + OUTPUT SUBMENU: + Output folder (string): APCEMM_out/ + Overwrite if folder exists (T/F): T + # FFT options (for spectral solver) + Use threaded FFT (T/F): F + FFTW WISDOM SUBMENU: + Use FFTW WISDOM (T/F): T + Dir w/ write permission (string): ./ + # This mostly contains information on background aerosol concentration; Not too relevant for contrail behavior + Input background condition (string): ../../input_data/init.txt + # All parameters here are overwritten in EMISSION INDICES SUBMENU + Input engine emissions (string): ../../input_data/ENG_EI.txt + # Ignore/Don't change these, these are deprecated features. + SAVE FORWARD RESULTS SUBMENU: + Save forward results (T/F): F + netCDF filename format (string): APCEMM_Case_* + ADJOINT OPTIMIZATION SUBMENU: + Turn on adjoint optim. (T/F): F + netCDF filename format (string): APCEMM_ADJ_Case_* + BOX MODEL SUBMENU: + Run box model (T/F): F + netCDF filename format (string): APCEMM_BOX_CASE_* + RANDOM NUMBER GENERATION SUBMENU: + Force seed value (T/F): T + Seed value (positive int): 0 + EPM type (original/external/new): original + +# Format of parameter items: +# Param name [unit] (Variable type) +PARAMETER MENU: + # Parameter sweep format : Format is either: x1 x2 x3 or start:increment:end + # : Example: 200 220 240 and 200:20:240 are identical + + # Monte Carlo simulation : min:max + # : Example: 200:240 will generate values for the parameter in between 200 and 240 + + # Maximum simulation time if contrail isn't gone by then: + Plume Process [hr] (double): 12 + + # Temperature, RH, and wind shear can be overwritten if using meteorological input files + METEOROLOGICAL PARAMETERS SUBMENU: + # Pressure altitude at which the contrail is initialized + Pressure [hPa] (double): 265 + Horiz. diff. coeff. [m^2/s] (double): 15.0 + # Can be overwritten if met file is passed + Verti. diff. [m^2/s] (double): 0.15 + Brunt-Vaisala Frequency [s^-1] (double): 0.013 + LOCATION AND TIME SUBMENU: + # Affects for solar zenith angle and photolysis which affect EPM/diurnal variations of temperature + LON [deg] (double): -15 + LAT [deg] (double): 60 + Emission day [1-365] (int): 81 + Emission time [hr] (double) : 8 + BACKGROUND MIXING RATIOS SUBMENU: + # Affects the EPM only + NOx [ppt] (double): 5100 + HNO3 [ppt] (double): 81.5 + O3 [ppb] (double): 100 + CO [ppb] (double): 40 + CH4 [ppm] (double): 1.76 + SO2 [ppt] (double): 7.25 + EMISSION INDICES SUBMENU: + # Engine specific parameters + # Default parameters are estimates for a 737-800 at 35k ft + # Affects EPM mostly: + NOx [g(NO2)/kg_fuel] (double): 10 # EDB for CFM56-5B3 + CO [g/kg_fuel] (double): 1 + UHC [g/kg_fuel] (double): 0.6 + # Affects contrail diffusion model: + SO2 [g/kg_fuel] (double): 1.2 # Assuming 600ppm of Sulfur + SO2 to SO4 conv [%] (double): 2 # AEDT paper, Barrett et al. (2010) + Soot [g/kg_fuel] (double): 0.008 # Assuming EI(nvPM) = 10^14 + Soot Radius [m] (double): 20.0E-09 + # Total fuel flow is then divided per engine + Total fuel flow [kg/s] (double) : 0.7 + Aircraft mass [kg] (double): 1.00e+05 # At beginning of cruise (35k ft), NPSS + Flight speed [m/s] (double): 265.42 # NPSS + Num. of engines [2/4] (int): 2 + Wingspan [m] (double): 34.32 # https://www.skybrary.aero/aircraft/b738 + # Values depend on where we assume condensation/freezing to occur and engine type + Core exit temp. [K] (double): 553.65 # NPSS + Exit bypass area [m^2] (double): 0.9772 + +TRANSPORT MENU: + # Keep on + Turn on Transport (T/F): T + # Outdated, not used (was used by spectral solver) + Fill Negative Values (T/F): T + Transport Timestep [min] (double): 1 + # Keep off: not sure of the effect yet + met updraft is included (if met file input) + PLUME UPDRAFT SUBMENU: + Turn on plume updraft (T/F): F + Updraft timescale [s] (double): 3600 + Updraft veloc. [cm/s] (double): 5 + +# Chemistry component of APCEMM hasn't been touched in a long time; leave off if only interested in contrail simulation +CHEMISTRY MENU: + Turn on Chemistry (T/F): F + Perform hetero. chem. (T/F): F + Chemistry Timestep [min] (double): 10 + Photolysis rates folder (string): /path/to/input/ + +AEROSOL MENU: + # Keep on + Turn on grav. settling (T/F): T + # Keep on + Turn on solid coagulation (T/F): T + # Keep off + Turn on liquid coagulation (T/F): F + Coag. timestep [min] (double): 60 + # Keep on + Turn on ice growth (T/F): T + Ice growth timestep [min] (double): 1 + +# At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true +# Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate + +# If using met. input: +# Exactly one of "Init temp. from met." and "Impose lapse rate" must be true +# Exactly one of "Init RH from met." and "Impose moist layer depth" must be true +METEOROLOGY MENU: + # --- MET INPUT OPTIONS ---+ + METEOROLOGICAL INPUT SUBMENU: + Use met. input (T/F): T + Met input file path (string): example_met_file_shear.nc + # Frequency of met data availability + Time series data timestep [hr] (double): 1.0 + # If off, uses the parameter specified in METEOROLOGICAL PARAMETERS SUBMENU + Init temp. from met. (T/F): T + Temp. time series input (T/F): T + # Always interpolates in time, this controls spatial interpolation + Interpolate temp. met. data (T/F): T + # Same options as the ones for Temperature: + Init RH from met. (T/F): T + RH time series input (T/F): T + Interpolate RH met. data (T/F): T + Init wind shear from met. (T/F): T + Wind shear time series input (T/F): T + Interpolate shear met. data (T/F): F + Init vert. veloc. from met. data (T/F): T + Vert. veloc. time series input (T/F): F + Interpolate vert. veloc. met. data (T/F): F + # Option to modify NWP RH data + + # Perturbations to the temperature field can be added to the contrail simulation to account for some effects of + # atmospheric turbulence and gravity waves. + # Every N minutes, a temperature perturbation of e1*e2*T_amp where e1 and e2 are random variables + # uniformly distributed between -1 and 1 generated individually for each grid cell. + # The importance of turbulence and grav. waves can be increased by increasing T_amp + # The relative importance of grav. waves vs. turb. is increased by increasing the timescale + # See Lewellen, Persistent Contrails and Contrail Cirrus (2014) for full details. + TEMPERATURE PERTURBATION SUBMENU: + Enable Temp. Pert. (T/F): F + Temp. Perturb. Amplitude (double): 1.0 + Temp. Perturb. Timescale (min): 10 + +# The only thing that should be changed here is the save frequency +DIAGNOSTIC MENU: + netCDF filename format (string): trac_avg.apcemm.hhmm + # Leave "save species timeseries" off. It will do nothing without also turning chemistry on. + SPECIES TIMESERIES SUBMENU: + Save species timeseries (T/F): F + Inst timeseries file (string): ts_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Species indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Index 1 (ice) is the only relevant thing to save + AEROSOL TIMESERIES SUBMENU: + Save aerosol timeseries (T/F): T + Inst timeseries file (string): ts_aerosol_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Aerosol indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Keep off if chemistry is also off + PRODUCTION & LOSS SUBMENU: + Turn on P/L diag (T/F): F + Save O3 P/L (T/F): F + +# Sometimes you have to change YLIM_DOWN here if the supersaturated layer is very thick +# because YLIM_DOWN must be larger than the layer thickness. +ADVANCED OPTIONS MENU: + # Domain is defined as: X [-XLIM_LEFT, XLIM_RIGHT], Y [-YLIM_DOWN, YLIM_UP] + GRID SUBMENU: + NX (positive int): 200 + NY (positive int): 180 + XLIM_RIGHT (positive double): 1.0e+3 + XLIM_LEFT (positive double): 1.0e+3 + YLIM_UP (positive double): 300 + YLIM_DOWN (positive double): 1.5e+3 + INITIAL CONTRAIL SIZE SUBMENU: + #Depth = BaseDepth + DepthScalingFactor * Default_Depth + #Same formula for width + Base Contrail Depth [m] (double): 0.0 + Contrail Depth Scaling Factor [-] (double): 1.0 + Base Contrail Width [m] (double): 0.0 + Contrail Width Scaling Factor [-] (double): 1.0 + Ambient Lapse Rate [K/km] (double): -3.0 + Tropopause Pressure [Pa] (double): 2.0e+4 + EARLY PLUME SUBMENU: + Reference ice crystal count [#/m] (double): 3.38e12 + Reference wingspan [m] (double): 60.3 + # The values below can be overriden if required. + Override post-jet ice crystal count (T/F): F + Post-jet ice crystal count [#/m] (double): 1e15 + diff --git a/examples/issl_rhi140/temperature.yaml b/examples/issl_rhi140/temperature.yaml new file mode 100644 index 00000000..52e5b3ad --- /dev/null +++ b/examples/issl_rhi140/temperature.yaml @@ -0,0 +1,216 @@ +SIMULATION MENU: + # Only one of parameter sweep or MC simulation can be on. + # Parameter sweep lets you specify an arbitrary number of custom values for each parameter; Monte Carlo simulation is self-explanatory. + # At the moment, you cannot mix and match MC sim and param sweep on each individual parameter. + OpenMP Num Threads (positive int): 1 + PARAM SWEEP SUBMENU: + Parameter sweep (T/F): T + #-OR--------------- + Run Monte Carlo (T/F): F + Num Monte Carlo runs (int): 2 + # Where APCEMM output for this set of runs will go + OUTPUT SUBMENU: + Output folder (string): APCEMM_out/ + Overwrite if folder exists (T/F): T + # FFT options (for spectral solver) + Use threaded FFT (T/F): F + FFTW WISDOM SUBMENU: + Use FFTW WISDOM (T/F): T + Dir w/ write permission (string): ./ + # This mostly contains information on background aerosol concentration; Not too relevant for contrail behavior + Input background condition (string): ../../input_data/init.txt + # All parameters here are overwritten in EMISSION INDICES SUBMENU + Input engine emissions (string): ../../input_data/ENG_EI.txt + # Ignore/Don't change these, these are deprecated features. + SAVE FORWARD RESULTS SUBMENU: + Save forward results (T/F): F + netCDF filename format (string): APCEMM_Case_* + ADJOINT OPTIMIZATION SUBMENU: + Turn on adjoint optim. (T/F): F + netCDF filename format (string): APCEMM_ADJ_Case_* + BOX MODEL SUBMENU: + Run box model (T/F): F + netCDF filename format (string): APCEMM_BOX_CASE_* + RANDOM NUMBER GENERATION SUBMENU: + Force seed value (T/F): T + Seed value (positive int): 0 + EPM type (original/external/new): original + +# Format of parameter items: +# Param name [unit] (Variable type) +PARAMETER MENU: + # Parameter sweep format : Format is either: x1 x2 x3 or start:increment:end + # : Example: 200 220 240 and 200:20:240 are identical + + # Monte Carlo simulation : min:max + # : Example: 200:240 will generate values for the parameter in between 200 and 240 + + # Maximum simulation time if contrail isn't gone by then: + Plume Process [hr] (double): 12 + + # Temperature, RH, and wind shear can be overwritten if using meteorological input files + METEOROLOGICAL PARAMETERS SUBMENU: + # Pressure altitude at which the contrail is initialized + Pressure [hPa] (double): 265 + Horiz. diff. coeff. [m^2/s] (double): 15.0 + # Can be overwritten if met file is passed + Verti. diff. [m^2/s] (double): 0.15 + Brunt-Vaisala Frequency [s^-1] (double): 0.013 + LOCATION AND TIME SUBMENU: + # Affects for solar zenith angle and photolysis which affect EPM/diurnal variations of temperature + LON [deg] (double): -15 + LAT [deg] (double): 60 + Emission day [1-365] (int): 81 + Emission time [hr] (double) : 8 + BACKGROUND MIXING RATIOS SUBMENU: + # Affects the EPM only + NOx [ppt] (double): 5100 + HNO3 [ppt] (double): 81.5 + O3 [ppb] (double): 100 + CO [ppb] (double): 40 + CH4 [ppm] (double): 1.76 + SO2 [ppt] (double): 7.25 + EMISSION INDICES SUBMENU: + # Engine specific parameters + # Default parameters are estimates for a 737-800 at 35k ft + # Affects EPM mostly: + NOx [g(NO2)/kg_fuel] (double): 10 # EDB for CFM56-5B3 + CO [g/kg_fuel] (double): 1 + UHC [g/kg_fuel] (double): 0.6 + # Affects contrail diffusion model: + SO2 [g/kg_fuel] (double): 1.2 # Assuming 600ppm of Sulfur + SO2 to SO4 conv [%] (double): 2 # AEDT paper, Barrett et al. (2010) + Soot [g/kg_fuel] (double): 0.008 # Assuming EI(nvPM) = 10^14 + Soot Radius [m] (double): 20.0E-09 + # Total fuel flow is then divided per engine + Total fuel flow [kg/s] (double) : 0.7 + Aircraft mass [kg] (double): 1.00e+05 # At beginning of cruise (35k ft), NPSS + Flight speed [m/s] (double): 265.42 # NPSS + Num. of engines [2/4] (int): 2 + Wingspan [m] (double): 34.32 # https://www.skybrary.aero/aircraft/b738 + # Values depend on where we assume condensation/freezing to occur and engine type + Core exit temp. [K] (double): 553.65 # NPSS + Exit bypass area [m^2] (double): 0.9772 + +TRANSPORT MENU: + # Keep on + Turn on Transport (T/F): T + # Outdated, not used (was used by spectral solver) + Fill Negative Values (T/F): T + Transport Timestep [min] (double): 1 + # Keep off: not sure of the effect yet + met updraft is included (if met file input) + PLUME UPDRAFT SUBMENU: + Turn on plume updraft (T/F): F + Updraft timescale [s] (double): 3600 + Updraft veloc. [cm/s] (double): 5 + +# Chemistry component of APCEMM hasn't been touched in a long time; leave off if only interested in contrail simulation +CHEMISTRY MENU: + Turn on Chemistry (T/F): F + Perform hetero. chem. (T/F): F + Chemistry Timestep [min] (double): 10 + Photolysis rates folder (string): /path/to/input/ + +AEROSOL MENU: + # Keep on + Turn on grav. settling (T/F): T + # Keep on + Turn on solid coagulation (T/F): T + # Keep off + Turn on liquid coagulation (T/F): F + Coag. timestep [min] (double): 60 + # Keep on + Turn on ice growth (T/F): T + Ice growth timestep [min] (double): 1 + +# At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true +# Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate + +# If using met. input: +# Exactly one of "Init temp. from met." and "Impose lapse rate" must be true +# Exactly one of "Init RH from met." and "Impose moist layer depth" must be true +METEOROLOGY MENU: + # --- MET INPUT OPTIONS ---+ + METEOROLOGICAL INPUT SUBMENU: + Use met. input (T/F): T + Met input file path (string): example_met_file_temperature.nc + # Frequency of met data availability + Time series data timestep [hr] (double): 1.0 + # If off, uses the parameter specified in METEOROLOGICAL PARAMETERS SUBMENU + Init temp. from met. (T/F): T + Temp. time series input (T/F): T + # Always interpolates in time, this controls spatial interpolation + Interpolate temp. met. data (T/F): T + # Same options as the ones for Temperature: + Init RH from met. (T/F): T + RH time series input (T/F): T + Interpolate RH met. data (T/F): T + Init wind shear from met. (T/F): T + Wind shear time series input (T/F): T + Interpolate shear met. data (T/F): F + Init vert. veloc. from met. data (T/F): T + Vert. veloc. time series input (T/F): F + Interpolate vert. veloc. met. data (T/F): F + # Option to modify NWP RH data + + # Perturbations to the temperature field can be added to the contrail simulation to account for some effects of + # atmospheric turbulence and gravity waves. + # Every N minutes, a temperature perturbation of e1*e2*T_amp where e1 and e2 are random variables + # uniformly distributed between -1 and 1 generated individually for each grid cell. + # The importance of turbulence and grav. waves can be increased by increasing T_amp + # The relative importance of grav. waves vs. turb. is increased by increasing the timescale + # See Lewellen, Persistent Contrails and Contrail Cirrus (2014) for full details. + TEMPERATURE PERTURBATION SUBMENU: + Enable Temp. Pert. (T/F): F + Temp. Perturb. Amplitude (double): 1.0 + Temp. Perturb. Timescale (min): 10 + +# The only thing that should be changed here is the save frequency +DIAGNOSTIC MENU: + netCDF filename format (string): trac_avg.apcemm.hhmm + # Leave "save species timeseries" off. It will do nothing without also turning chemistry on. + SPECIES TIMESERIES SUBMENU: + Save species timeseries (T/F): F + Inst timeseries file (string): ts_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Species indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Index 1 (ice) is the only relevant thing to save + AEROSOL TIMESERIES SUBMENU: + Save aerosol timeseries (T/F): T + Inst timeseries file (string): ts_aerosol_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Aerosol indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Keep off if chemistry is also off + PRODUCTION & LOSS SUBMENU: + Turn on P/L diag (T/F): F + Save O3 P/L (T/F): F + +# Sometimes you have to change YLIM_DOWN here if the supersaturated layer is very thick +# because YLIM_DOWN must be larger than the layer thickness. +ADVANCED OPTIONS MENU: + # Domain is defined as: X [-XLIM_LEFT, XLIM_RIGHT], Y [-YLIM_DOWN, YLIM_UP] + GRID SUBMENU: + NX (positive int): 200 + NY (positive int): 180 + XLIM_RIGHT (positive double): 1.0e+3 + XLIM_LEFT (positive double): 1.0e+3 + YLIM_UP (positive double): 300 + YLIM_DOWN (positive double): 1.5e+3 + INITIAL CONTRAIL SIZE SUBMENU: + #Depth = BaseDepth + DepthScalingFactor * Default_Depth + #Same formula for width + Base Contrail Depth [m] (double): 0.0 + Contrail Depth Scaling Factor [-] (double): 1.0 + Base Contrail Width [m] (double): 0.0 + Contrail Width Scaling Factor [-] (double): 1.0 + Ambient Lapse Rate [K/km] (double): -3.0 + Tropopause Pressure [Pa] (double): 2.0e+4 + EARLY PLUME SUBMENU: + Reference ice crystal count [#/m] (double): 3.38e12 + Reference wingspan [m] (double): 60.3 + # The values below can be overriden if required. + Override post-jet ice crystal count (T/F): F + Post-jet ice crystal count [#/m] (double): 1e15 + diff --git a/examples/issl_rhi140/w.yaml b/examples/issl_rhi140/w.yaml new file mode 100644 index 00000000..8aa05b06 --- /dev/null +++ b/examples/issl_rhi140/w.yaml @@ -0,0 +1,216 @@ +SIMULATION MENU: + # Only one of parameter sweep or MC simulation can be on. + # Parameter sweep lets you specify an arbitrary number of custom values for each parameter; Monte Carlo simulation is self-explanatory. + # At the moment, you cannot mix and match MC sim and param sweep on each individual parameter. + OpenMP Num Threads (positive int): 1 + PARAM SWEEP SUBMENU: + Parameter sweep (T/F): T + #-OR--------------- + Run Monte Carlo (T/F): F + Num Monte Carlo runs (int): 2 + # Where APCEMM output for this set of runs will go + OUTPUT SUBMENU: + Output folder (string): APCEMM_out/ + Overwrite if folder exists (T/F): T + # FFT options (for spectral solver) + Use threaded FFT (T/F): F + FFTW WISDOM SUBMENU: + Use FFTW WISDOM (T/F): T + Dir w/ write permission (string): ./ + # This mostly contains information on background aerosol concentration; Not too relevant for contrail behavior + Input background condition (string): ../../input_data/init.txt + # All parameters here are overwritten in EMISSION INDICES SUBMENU + Input engine emissions (string): ../../input_data/ENG_EI.txt + # Ignore/Don't change these, these are deprecated features. + SAVE FORWARD RESULTS SUBMENU: + Save forward results (T/F): F + netCDF filename format (string): APCEMM_Case_* + ADJOINT OPTIMIZATION SUBMENU: + Turn on adjoint optim. (T/F): F + netCDF filename format (string): APCEMM_ADJ_Case_* + BOX MODEL SUBMENU: + Run box model (T/F): F + netCDF filename format (string): APCEMM_BOX_CASE_* + RANDOM NUMBER GENERATION SUBMENU: + Force seed value (T/F): T + Seed value (positive int): 0 + EPM type (original/external/new): original + +# Format of parameter items: +# Param name [unit] (Variable type) +PARAMETER MENU: + # Parameter sweep format : Format is either: x1 x2 x3 or start:increment:end + # : Example: 200 220 240 and 200:20:240 are identical + + # Monte Carlo simulation : min:max + # : Example: 200:240 will generate values for the parameter in between 200 and 240 + + # Maximum simulation time if contrail isn't gone by then: + Plume Process [hr] (double): 12 + + # Temperature, RH, and wind shear can be overwritten if using meteorological input files + METEOROLOGICAL PARAMETERS SUBMENU: + # Pressure altitude at which the contrail is initialized + Pressure [hPa] (double): 265 + Horiz. diff. coeff. [m^2/s] (double): 15.0 + # Can be overwritten if met file is passed + Verti. diff. [m^2/s] (double): 0.15 + Brunt-Vaisala Frequency [s^-1] (double): 0.013 + LOCATION AND TIME SUBMENU: + # Affects for solar zenith angle and photolysis which affect EPM/diurnal variations of temperature + LON [deg] (double): -15 + LAT [deg] (double): 60 + Emission day [1-365] (int): 81 + Emission time [hr] (double) : 8 + BACKGROUND MIXING RATIOS SUBMENU: + # Affects the EPM only + NOx [ppt] (double): 5100 + HNO3 [ppt] (double): 81.5 + O3 [ppb] (double): 100 + CO [ppb] (double): 40 + CH4 [ppm] (double): 1.76 + SO2 [ppt] (double): 7.25 + EMISSION INDICES SUBMENU: + # Engine specific parameters + # Default parameters are estimates for a 737-800 at 35k ft + # Affects EPM mostly: + NOx [g(NO2)/kg_fuel] (double): 10 # EDB for CFM56-5B3 + CO [g/kg_fuel] (double): 1 + UHC [g/kg_fuel] (double): 0.6 + # Affects contrail diffusion model: + SO2 [g/kg_fuel] (double): 1.2 # Assuming 600ppm of Sulfur + SO2 to SO4 conv [%] (double): 2 # AEDT paper, Barrett et al. (2010) + Soot [g/kg_fuel] (double): 0.008 # Assuming EI(nvPM) = 10^14 + Soot Radius [m] (double): 20.0E-09 + # Total fuel flow is then divided per engine + Total fuel flow [kg/s] (double) : 0.7 + Aircraft mass [kg] (double): 1.00e+05 # At beginning of cruise (35k ft), NPSS + Flight speed [m/s] (double): 265.42 # NPSS + Num. of engines [2/4] (int): 2 + Wingspan [m] (double): 34.32 # https://www.skybrary.aero/aircraft/b738 + # Values depend on where we assume condensation/freezing to occur and engine type + Core exit temp. [K] (double): 553.65 # NPSS + Exit bypass area [m^2] (double): 0.9772 + +TRANSPORT MENU: + # Keep on + Turn on Transport (T/F): T + # Outdated, not used (was used by spectral solver) + Fill Negative Values (T/F): T + Transport Timestep [min] (double): 1 + # Keep off: not sure of the effect yet + met updraft is included (if met file input) + PLUME UPDRAFT SUBMENU: + Turn on plume updraft (T/F): F + Updraft timescale [s] (double): 3600 + Updraft veloc. [cm/s] (double): 5 + +# Chemistry component of APCEMM hasn't been touched in a long time; leave off if only interested in contrail simulation +CHEMISTRY MENU: + Turn on Chemistry (T/F): F + Perform hetero. chem. (T/F): F + Chemistry Timestep [min] (double): 10 + Photolysis rates folder (string): /path/to/input/ + +AEROSOL MENU: + # Keep on + Turn on grav. settling (T/F): T + # Keep on + Turn on solid coagulation (T/F): T + # Keep off + Turn on liquid coagulation (T/F): F + Coag. timestep [min] (double): 60 + # Keep on + Turn on ice growth (T/F): T + Ice growth timestep [min] (double): 1 + +# At least one of "Use met. input", "Impose moist layer depth", or "Impose lapse rate" must be true +# Imposing moist layer depth will automatically calculate the lapse rate and override the imposed lapse rate + +# If using met. input: +# Exactly one of "Init temp. from met." and "Impose lapse rate" must be true +# Exactly one of "Init RH from met." and "Impose moist layer depth" must be true +METEOROLOGY MENU: + # --- MET INPUT OPTIONS ---+ + METEOROLOGICAL INPUT SUBMENU: + Use met. input (T/F): T + Met input file path (string): example_met_file_w_01.nc + # Frequency of met data availability + Time series data timestep [hr] (double): 1.0 + # If off, uses the parameter specified in METEOROLOGICAL PARAMETERS SUBMENU + Init temp. from met. (T/F): T + Temp. time series input (T/F): T + # Always interpolates in time, this controls spatial interpolation + Interpolate temp. met. data (T/F): T + # Same options as the ones for Temperature: + Init RH from met. (T/F): T + RH time series input (T/F): T + Interpolate RH met. data (T/F): T + Init wind shear from met. (T/F): T + Wind shear time series input (T/F): T + Interpolate shear met. data (T/F): F + Init vert. veloc. from met. data (T/F): T + Vert. veloc. time series input (T/F): F + Interpolate vert. veloc. met. data (T/F): F + # Option to modify NWP RH data + + # Perturbations to the temperature field can be added to the contrail simulation to account for some effects of + # atmospheric turbulence and gravity waves. + # Every N minutes, a temperature perturbation of e1*e2*T_amp where e1 and e2 are random variables + # uniformly distributed between -1 and 1 generated individually for each grid cell. + # The importance of turbulence and grav. waves can be increased by increasing T_amp + # The relative importance of grav. waves vs. turb. is increased by increasing the timescale + # See Lewellen, Persistent Contrails and Contrail Cirrus (2014) for full details. + TEMPERATURE PERTURBATION SUBMENU: + Enable Temp. Pert. (T/F): F + Temp. Perturb. Amplitude (double): 1.0 + Temp. Perturb. Timescale (min): 10 + +# The only thing that should be changed here is the save frequency +DIAGNOSTIC MENU: + netCDF filename format (string): trac_avg.apcemm.hhmm + # Leave "save species timeseries" off. It will do nothing without also turning chemistry on. + SPECIES TIMESERIES SUBMENU: + Save species timeseries (T/F): F + Inst timeseries file (string): ts_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Species indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Index 1 (ice) is the only relevant thing to save + AEROSOL TIMESERIES SUBMENU: + Save aerosol timeseries (T/F): T + Inst timeseries file (string): ts_aerosol_hhmm.nc + #list input: separate by spaces. e.g. 1 2 3 4 5 + Aerosol indices to include (list of ints): 1 + Save frequency [min] (double): 10 + # Keep off if chemistry is also off + PRODUCTION & LOSS SUBMENU: + Turn on P/L diag (T/F): F + Save O3 P/L (T/F): F + +# Sometimes you have to change YLIM_DOWN here if the supersaturated layer is very thick +# because YLIM_DOWN must be larger than the layer thickness. +ADVANCED OPTIONS MENU: + # Domain is defined as: X [-XLIM_LEFT, XLIM_RIGHT], Y [-YLIM_DOWN, YLIM_UP] + GRID SUBMENU: + NX (positive int): 200 + NY (positive int): 180 + XLIM_RIGHT (positive double): 1.0e+3 + XLIM_LEFT (positive double): 1.0e+3 + YLIM_UP (positive double): 300 + YLIM_DOWN (positive double): 1.5e+3 + INITIAL CONTRAIL SIZE SUBMENU: + #Depth = BaseDepth + DepthScalingFactor * Default_Depth + #Same formula for width + Base Contrail Depth [m] (double): 0.0 + Contrail Depth Scaling Factor [-] (double): 1.0 + Base Contrail Width [m] (double): 0.0 + Contrail Width Scaling Factor [-] (double): 1.0 + Ambient Lapse Rate [K/km] (double): -3.0 + Tropopause Pressure [Pa] (double): 2.0e+4 + EARLY PLUME SUBMENU: + Reference ice crystal count [#/m] (double): 3.38e12 + Reference wingspan [m] (double): 60.3 + # The values below can be overriden if required. + Override post-jet ice crystal count (T/F): F + Post-jet ice crystal count [#/m] (double): 1e15 + From 1d358b5582bd3d179995719aa65b6c72aadc229c Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Mon, 16 Feb 2026 14:32:42 +0000 Subject: [PATCH 3/5] Reduced usage of pow() (and log()) * Add profiling flags * Replaced pow() * Aerosol * Added more optimisations * Code clean-up --------- Co-authored-by: Coco Yeung --- Code.v05-00/include/AIM/Aerosol.hpp | 6 +- Code.v05-00/include/Util/MetFunction.hpp | 10 +- Code.v05-00/src/AIM/Aerosol.cpp | 134 +++++++++++++++--- Code.v05-00/src/AIM/Settling.cpp | 5 +- Code.v05-00/src/AIM/buildKernel.cpp | 27 ++-- Code.v05-00/src/Core/Engine.cpp | 7 +- Code.v05-00/src/Core/LAGRIDPlumeModel.cpp | 5 +- Code.v05-00/src/Core/LiquidAer.cpp | 24 ++-- Code.v05-00/src/Core/Vortex.cpp | 6 +- .../src/EPM/Models/Original/Integrate.cpp | 11 +- Code.v05-00/src/EPM/Solution.cpp | 13 +- Code.v05-00/src/KPP/KPP_HetRates.cpp | 28 ++-- Code.v05-00/src/KPP/KPP_Integrator_ADJ.cpp | 9 +- Code.v05-00/src/KPP/KPP_LinearAlgebra.cpp | 5 +- Code.v05-00/src/KPP/KPP_Main_ADJ.cpp | 24 ++-- Code.v05-00/src/KPP/KPP_Rates.cpp | 10 +- Code.v05-00/src/LAGRID/RemappingFunctions.cpp | 4 +- Code.v05-00/src/Util/PhysFunction.cpp | 10 +- 18 files changed, 234 insertions(+), 104 deletions(-) diff --git a/Code.v05-00/include/AIM/Aerosol.hpp b/Code.v05-00/include/AIM/Aerosol.hpp index 77b5c56f..a1bc0299 100644 --- a/Code.v05-00/include/AIM/Aerosol.hpp +++ b/Code.v05-00/include/AIM/Aerosol.hpp @@ -62,7 +62,11 @@ class AIM::Aerosol /* Moments */ //NOTE: This gives the moment in [- / cm3]. You need to multiply by factors to get it in [ / m] or something. inline double binMoment(int iBin, int n = 0) const { - return (log(bin_Edges[iBin + 1]) - log(bin_Edges[iBin])) * pow(bin_Centers[iBin], n) * pdf[iBin]; + if (n == 0){ + return pdf[iBin]; // if n = 0, it can skip the pow() calculation + } + + return pow(bin_Centers[iBin], n) * pdf[iBin]; } //NOTE: This gives the moment in [- / cm3]. You need to multiply by factors to get it in [ / m] or something. double Moment( UInt n = 0 ) const; diff --git a/Code.v05-00/include/Util/MetFunction.hpp b/Code.v05-00/include/Util/MetFunction.hpp index 0c6cc4e6..cd20dc4e 100644 --- a/Code.v05-00/include/Util/MetFunction.hpp +++ b/Code.v05-00/include/Util/MetFunction.hpp @@ -65,12 +65,14 @@ namespace met Atmospheric Chemistry and Physics https://doi.org/10.5194/acp-22-10919-2022 */ + + double rhi_abs = rhi * 0.01; - double rhi_abs = rhi / 100.0; + double rhi_abs_a = rhi_abs / a; - double rhi_corr = (rhi_abs / a) <= 1 - ? rhi_abs / a - : std::min(std::pow(rhi_abs / a, b), 1.65); + double rhi_corr = (rhi_abs_a) <= 1 + ? rhi_abs_a + : std::min(std::pow(rhi_abs_a, b), 1.65); return rhi_corr * 100.0; } struct newXCoordsPair { diff --git a/Code.v05-00/src/AIM/Aerosol.cpp b/Code.v05-00/src/AIM/Aerosol.cpp index 79523ab3..d6aa6706 100644 --- a/Code.v05-00/src/AIM/Aerosol.cpp +++ b/Code.v05-00/src/AIM/Aerosol.cpp @@ -48,7 +48,7 @@ namespace AIM /* Compute size of each bin */ for (UInt iBin = 0; iBin < nBin; iBin++) { bin_Sizes[iBin] = bin_Edges[iBin + 1] - bin_Edges[iBin]; - bin_VCenters[iBin] = 4.0 / 3.0 * PI * (pow(bin_Edges_[iBin], 3) + pow(bin_Edges_[iBin + 1], 3)) * 0.5; + bin_VCenters[iBin] = 4.0 / 3.0 * PI * (bin_Edges_[iBin] * bin_Edges_[iBin] * bin_Edges_[iBin] + bin_Edges_[iBin + 1] * bin_Edges_[iBin + 1] * bin_Edges_[iBin + 1]) * 0.5; } /* Check mean and standard deviation */ @@ -84,8 +84,9 @@ namespace AIM if (alpha <= 0.0) { std::cout << "\nIn Aerosol::Aerosol: power law requires that alpha > 0 ( alpha = " << alpha << " )\n"; } + double inv_bin_Centers_0 = 1/bin_Centers[0]; for (UInt iBin = 0; iBin < bin_Centers.size(); iBin++) { - pdf[iBin] = nPart_ * alpha * pow(bin_Centers[iBin] / bin_Centers[0], -alpha); + pdf[iBin] = nPart_ * alpha * pow(bin_Centers[iBin] * inv_bin_Centers_0, -alpha); } } else if ((strcmp(type, "gam") == 0) || (strcmp(type, "gamma") == 0) || (strcmp(type, "generalized gamma") == 0)) { /* Gamma distribution: @@ -96,9 +97,10 @@ namespace AIM if ((alpha <= 0) || (std::fmod(alpha, 1.0) != 0)) { std::cout << "\nIn Aerosol::Aerosol: (generalized) gamma distribution requires that alpha is a positive integer ( alpha = " << alpha << " )\n"; } if (gamma <= 0) { std::cout << "\nIn Aerosol::Aerosol: (generalized) gamma distribution requires that gamma is positive ( gamma = " << gamma << " )\n"; } if (b <= 0) { std::cout << "\nIn Aerosol::Aerosol: (generalized) gamma distribution requires that b is positive ( b = " << b << " )\n"; } - + + double pow_value = pow(b, (alpha + 1) / gamma); for (UInt iBin = 0; iBin < bin_Centers.size(); iBin++) { - pdf[iBin] = nPart_ * gamma * pow(b, (alpha + 1) / gamma) / boost::math::tgamma((alpha + 1) / gamma) * pow(bin_Centers[iBin], alpha + 1) * exp(-b * pow(bin_Centers[iBin], gamma)); + pdf[iBin] = nPart_ * gamma * pow_value / boost::math::tgamma((alpha + 1) / gamma) * pow(bin_Centers[iBin], alpha + 1) * exp(-b * pow(bin_Centers[iBin], gamma)); } } else { std::cout << "\nIn Aerosol::Aerosol: distribution type must be either lognormal, normal, power or (generalized) gamma\n"; @@ -117,7 +119,7 @@ namespace AIM pdf(pdf) { for (UInt iBin = 0; iBin < nBin; iBin++) { bin_Sizes[iBin] = bin_Edges[iBin + 1] - bin_Edges[iBin]; - bin_VCenters[iBin] = 4.0 / 3.0 * PI * (pow(bin_Edges[iBin], 3) + pow(bin_Edges[iBin + 1], 3)) * 0.5; + bin_VCenters[iBin] = 4.0 / 3.0 * PI * (bin_Edges[iBin] * bin_Edges[iBin] * bin_Edges[iBin] + bin_Edges[iBin + 1] * bin_Edges[iBin + 1] * bin_Edges[iBin + 1]) * 0.5; } } @@ -304,10 +306,32 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { { double moment = 0; + double pow_value; + double log_current = log(bin_Edges[0]); + double log_next; for (UInt iBin = 0; iBin < nBin; iBin++) { - moment += (log(bin_Edges[iBin + 1]) - log(bin_Edges[iBin])) * pow(bin_Centers[iBin], n) * pdf[iBin]; + log_next = log(bin_Edges[iBin + 1]); + switch(n) { + case 0: + pow_value = 1; + break; + case 1: + pow_value = bin_Centers[iBin]; + break; + case 2: + pow_value = bin_Centers[iBin] * bin_Centers[iBin]; + break; + case 3: + pow_value = bin_Centers[iBin] * bin_Centers[iBin] * bin_Centers[iBin]; + break; + default: + pow_value = pow(bin_Centers[iBin], n); // fallback for safety + break; + } + moment += (log_next - log_current) * pow_value * pdf[iBin]; + log_current = log_next; // reuse for next iteration } return moment; @@ -355,7 +379,8 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { if (N > 0) { - return sqrt(Moment(2) / N - pow(Moment(1) / N, 2.0)); + double Moment_1_N = Moment(1) / N; + return sqrt(Moment(2) / N - Moment_1_N * Moment_1_N); } else { @@ -417,10 +442,10 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { /* Compute size of each bin */ for (UInt iBin = 0; iBin < nBin; iBin++) { - bin_VEdges[iBin] = 4.0 / 3.0 * PI * pow(bin_Edges[iBin], 3); + bin_VEdges[iBin] = 4.0 / 3.0 * PI * bin_Edges[iBin] * bin_Edges[iBin] * bin_Edges[iBin]; bin_Sizes[iBin] = bin_Edges[iBin + 1] - bin_Edges[iBin]; } - bin_VEdges[nBin] = 4.0 / 3.0 * PI * pow(bin_Edges[nBin], 3); + bin_VEdges[nBin] = 4.0 / 3.0 * PI * bin_Edges[nBin] * bin_Edges[nBin] * bin_Edges[nBin]; bin_VCenters.resize(nBin, Vector_2D(Ny, Vector_1D(Nx, 0.0E+00))); @@ -485,10 +510,11 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { if (alpha <= 0.0) { std::cout << "\nIn Grid_Aerosol::Grid_Aerosol: power law requires that alpha > 0 ( alpha = " << alpha << " )\n"; } - + + double inv_bin_Centers_0 = 1/ bin_Centers[0]; for (UInt iBin = 0; iBin < bin_Centers.size(); iBin++) { - pdf[iBin][0][0] = nPart_ * alpha * pow(bin_Centers[iBin] / bin_Centers[0], -alpha); + pdf[iBin][0][0] = nPart_ * alpha * pow(bin_Centers[iBin] * inv_bin_Centers_0, -alpha); for (UInt jNy = 0; jNy < Ny; jNy++) { for (UInt iNx = 0; iNx < Nx; iNx++) @@ -509,9 +535,10 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { if (gamma <= 0) { std::cout << "\nIn Grid_Aerosol::Grid_Aerosol: (generalized) gamma distribution requires that gamma is positive ( gamma = " << gamma << " )\n"; } if (b <= 0) { std::cout << "\nIn Grid_Aerosol::Grid_Aerosol: (generalized) gamma distribution requires that b is positive ( b = " << b << " )\n"; } + double pow_value = pow(b, (alpha + 1) / gamma); for (UInt iBin = 0; iBin < bin_Centers.size(); iBin++) { - pdf[iBin][0][0] = nPart_ * gamma * pow(b, (alpha + 1) / gamma) / boost::math::tgamma((alpha + 1) / gamma) * pow(bin_Centers[iBin], alpha + 1) * exp(-b * pow(bin_Centers[iBin], gamma)); + pdf[iBin][0][0] = nPart_ * gamma * pow_value / boost::math::tgamma((alpha + 1) / gamma) * pow(bin_Centers[iBin], alpha + 1) * exp(-b * pow(bin_Centers[iBin], gamma)); for (UInt jNy = 0; jNy < Ny; jNy++) { for (UInt iNx = 0; iNx < Nx; iNx++) @@ -1128,17 +1155,40 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { Vector_2D moment(Ny, Vector_1D(Nx, 0.0E+00)); const double FACTOR = 3.0 / double(4.0 * PI); + double log_current = log(bin_Edges[0]); + double log_next; + double pow_value; + #pragma omp parallel for default(shared) private(iNx, jNy, iBin) \ schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++) { + log_next = log(bin_Edges[iBin + 1]); for (jNy = 0; jNy < Ny; jNy++) { for (iNx = 0; iNx < Nx; iNx++) { - moment[jNy][iNx] += (log(bin_Edges[iBin + 1] / bin_Edges[iBin])) * pow(FACTOR * bin_VCenters[iBin][jNy][iNx], n / 3.0) * pdf[iBin][jNy][iNx]; + switch(n) { + case 0: + pow_value = 1; + break; + case 1: + pow_value = cbrt(FACTOR * bin_VCenters[iBin][jNy][iNx]); + break; + case 2: + pow_value = cbrt(FACTOR * FACTOR * bin_VCenters[iBin][jNy][iNx] * bin_VCenters[iBin][jNy][iNx]); + break; + case 3: + pow_value = FACTOR * bin_VCenters[iBin][jNy][iNx]; + break; + default: + pow_value = pow(FACTOR * bin_VCenters[iBin][jNy][iNx], n / 3.0); // fallback for safety + break; + } + moment[jNy][iNx] += (log_next - log_current) * pow_value * pdf[iBin][jNy][iNx]; } } + log_current = log_next; // reuse for next iteration } return moment; @@ -1619,13 +1669,38 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { double moment = 0.0E+00; + double log_current = log(bin_Edges[0]); + double log_next; + double pow_value; + #pragma omp parallel for default(shared) private(iBin) \ reduction(+ \ : moment) \ schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++) { - moment += (log(bin_Edges[iBin + 1] / bin_Edges[iBin])) * pow(bin_Centers[iBin], n) * PDF[iBin]; + log_next = log(bin_Edges[iBin+1]); + + switch(n) { + case 0: + pow_value = 1; + break; + case 1: + pow_value = bin_Centers[iBin]; + break; + case 2: + pow_value = bin_Centers[iBin] * bin_Centers[iBin]; + break; + case 3: + pow_value = bin_Centers[iBin] * bin_Centers[iBin] * bin_Centers[iBin]; + break; + default: + pow_value = pow(bin_Centers[iBin], n); // fallback for safety + break; + } + + moment += (log_next - log_current) * pow_value * PDF[iBin]; + log_current = log_next; } return moment; @@ -1640,12 +1715,36 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { double moment = 0.0E+00; const double FACTOR = 3.0 / double(4.0 * PI); + double pow_value; + double log_current = log(bin_Edges[0]); + double log_next; + #pragma omp parallel for default(shared) private(iBin) \ reduction(+ \ : moment) \ schedule(dynamic, 1) if (!PARALLEL_CASES) - for (iBin = 0; iBin < nBin; iBin++) - moment += (log(bin_Edges[iBin + 1] / bin_Edges[iBin])) * pow(FACTOR * bin_VCenters[iBin][jNy][iNx], n / double(3.0)) * pdf[iBin][jNy][iNx]; + for (iBin = 0; iBin < nBin; iBin++){ + log_next = log(bin_Edges[iBin+1]); + switch(n) { + case 0: + pow_value = 1; + break; + case 1: + pow_value = cbrt(FACTOR * bin_VCenters[iBin][jNy][iNx]); + break; + case 2: + pow_value = cbrt(FACTOR * FACTOR * bin_VCenters[iBin][jNy][iNx] * bin_VCenters[iBin][jNy][iNx]); + break; + case 3: + pow_value = FACTOR * bin_VCenters[iBin][jNy][iNx]; + break; + default: + pow_value = pow(FACTOR * bin_VCenters[iBin][jNy][iNx], n / 3.0); // fallback for safety + break; + } + moment += (log_next - log_current) * pow_value * pdf[iBin][jNy][iNx]; + log_current = log_next; + } return moment; @@ -1692,7 +1791,8 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { if (N > 0) { - return sqrt(Moment(2, jNy, iNx) / N - pow(Moment(1, jNy, iNx) / N, 2.0)); + double Moment_1_N = Moment(1, jNy, iNx) / N; + return sqrt(Moment(2, jNy, iNx) / N - Moment_1_N * Moment_1_N); } else { diff --git a/Code.v05-00/src/AIM/Settling.cpp b/Code.v05-00/src/AIM/Settling.cpp index e4b206f5..2fd3ccd6 100644 --- a/Code.v05-00/src/AIM/Settling.cpp +++ b/Code.v05-00/src/AIM/Settling.cpp @@ -97,8 +97,9 @@ namespace AIM } X = 8.0E+00 * physConst::g * rhoA / eta2 * binCenters[iBin] * binCenters[iBin] * mi_Ai; - - Re = C2 * pow(sqrt(1.0E+00 + C1 * sqrt(X)) - 1.0E+00, 2.0) - a0 * pow( X, b0 ); + + double nested_sqrt = sqrt(1.0E+00 + C1 * sqrt(X)) - 1.0E+00; + Re = C2 * nested_sqrt * nested_sqrt - a0 * pow( X, b0 ); vFall[iBin] = Re * eta / ( rhoA * 2.0E+00 * binCenters[iBin] ); diff --git a/Code.v05-00/src/AIM/buildKernel.cpp b/Code.v05-00/src/AIM/buildKernel.cpp index 8dee9e3b..ffb6cdcf 100644 --- a/Code.v05-00/src/AIM/buildKernel.cpp +++ b/Code.v05-00/src/AIM/buildKernel.cpp @@ -141,15 +141,15 @@ namespace AIM for ( unsigned int iBin = 0; iBin < bin_Centers.size(); iBin++ ) { if ( bin_Centers[iBin] <= bin_R ) { if ( physFunc::Reynolds_p( bin_R, rho_2, temperature_K, pressure_Pa ) <= 1 ) { - K_DE[iBin] = K_Brow[iBin] * 0.45 * pow( physFunc::Reynolds_p( bin_R, rho_2, temperature_K, pressure_Pa ), 1.0 / double(3.0) ) * pow( physFunc::Schmidt_p( bin_Centers[iBin], temperature_K, pressure_Pa ), 1.0 / double(3.0) ); + K_DE[iBin] = K_Brow[iBin] * 0.45 * cbrt( physFunc::Reynolds_p( bin_R, rho_2, temperature_K, pressure_Pa )) * cbrt( physFunc::Schmidt_p( bin_Centers[iBin], temperature_K, pressure_Pa )); } else { - K_DE[iBin] = K_Brow[iBin] * 0.45 * pow( physFunc::Reynolds_p( bin_R, rho_2, temperature_K, pressure_Pa ), 0.5 ) * pow( physFunc::Schmidt_p( bin_Centers[iBin], temperature_K, pressure_Pa ), 1.0 / double(3.0) ); + K_DE[iBin] = K_Brow[iBin] * 0.45 * sqrt( physFunc::Reynolds_p( bin_R, rho_2, temperature_K, pressure_Pa )) * cbrt( physFunc::Schmidt_p( bin_Centers[iBin], temperature_K, pressure_Pa )); } } else { if ( physFunc::Reynolds_p( bin_Centers[iBin], rho_1, temperature_K, pressure_Pa ) <= 1 ) { - K_DE[iBin] = K_Brow[iBin] * 0.45 * pow( physFunc::Reynolds_p( bin_Centers[iBin], rho_1, temperature_K, pressure_Pa ), 1.0 / double(3.0) ) * pow( physFunc::Schmidt_p( bin_R, temperature_K, pressure_Pa ), 1.0 / double(3.0) ); + K_DE[iBin] = K_Brow[iBin] * 0.45 * cbrt( physFunc::Reynolds_p( bin_Centers[iBin], rho_1, temperature_K, pressure_Pa )) * cbrt( physFunc::Schmidt_p( bin_R, temperature_K, pressure_Pa )); } else { - K_DE[iBin] = K_Brow[iBin] * 0.45 * pow( physFunc::Reynolds_p( bin_Centers[iBin], rho_1, temperature_K, pressure_Pa ), 0.5 ) * pow( physFunc::Schmidt_p( bin_R, temperature_K, pressure_Pa ), 1.0 / double(3.0) ); + K_DE[iBin] = K_Brow[iBin] * 0.45 * sqrt( physFunc::Reynolds_p( bin_Centers[iBin], rho_1, temperature_K, pressure_Pa )) * cbrt( physFunc::Schmidt_p( bin_R, temperature_K, pressure_Pa )); } } } @@ -174,17 +174,20 @@ namespace AIM for ( unsigned int iBin_1 = 0; iBin_1 < bin_Centers_1.size(); iBin_1++ ) { K_DE.push_back( Vector_1D( bin_Centers_2.size() ) ); for ( unsigned int iBin_2 = 0; iBin_2 < bin_Centers_2.size(); iBin_2++ ) { + double reynolds_p_1 = physFunc::Reynolds_p( bin_Centers_1[iBin_1], rho_1, temperature_K, pressure_Pa); + double reynolds_p_2 = physFunc::Reynolds_p( bin_Centers_2[iBin_2], rho_2, temperature_K, pressure_Pa); + if ( bin_Centers_1[iBin_1] <= bin_Centers_2[iBin_2]) { - if ( physFunc::Reynolds_p( bin_Centers_2[iBin_2], rho_2, temperature_K, pressure_Pa ) <= 1 ) { - K_DE[iBin_1][iBin_2] = K_Brow[iBin_1][iBin_2] * 0.45 * pow( physFunc::Reynolds_p( bin_Centers_2[iBin_2], rho_2, temperature_K, pressure_Pa ), 1.0 / double(3.0) ) * pow( physFunc::Schmidt_p( bin_Centers_1[iBin_1], temperature_K, pressure_Pa ), 1.0 / double(3.0) ); + if ( reynolds_p_2 <= 1 ) { + K_DE[iBin_1][iBin_2] = K_Brow[iBin_1][iBin_2] * 0.45 * cbrt( reynolds_p_2 ) * cbrt( physFunc::Schmidt_p( bin_Centers_1[iBin_1], temperature_K, pressure_Pa )); } else { - K_DE[iBin_1][iBin_2] = K_Brow[iBin_1][iBin_2] * 0.45 * pow( physFunc::Reynolds_p( bin_Centers_2[iBin_2], rho_2, temperature_K, pressure_Pa ), 0.5 ) * pow( physFunc::Schmidt_p( bin_Centers_1[iBin_1], temperature_K, pressure_Pa ), 1.0 / double(3.0) ); + K_DE[iBin_1][iBin_2] = K_Brow[iBin_1][iBin_2] * 0.45 * sqrt( reynolds_p_2 ) * cbrt( physFunc::Schmidt_p( bin_Centers_1[iBin_1], temperature_K, pressure_Pa )); } } else { - if ( physFunc::Reynolds_p( bin_Centers_1[iBin_1], rho_1, temperature_K, pressure_Pa ) <= 1 ) { - K_DE[iBin_1][iBin_2] = K_Brow[iBin_1][iBin_2] * 0.45 * pow( physFunc::Reynolds_p( bin_Centers_1[iBin_1], rho_1, temperature_K, pressure_Pa ), 1.0 / double(3.0) ) * pow( physFunc::Schmidt_p( bin_Centers_2[iBin_2], temperature_K, pressure_Pa ), 1.0 / double(3.0) ); + if ( reynolds_p_1 <= 1 ) { + K_DE[iBin_1][iBin_2] = K_Brow[iBin_1][iBin_2] * 0.45 * cbrt( reynolds_p_1 ) * cbrt( physFunc::Schmidt_p( bin_Centers_2[iBin_2], temperature_K, pressure_Pa )); } else { - K_DE[iBin_1][iBin_2] = K_Brow[iBin_1][iBin_2] * 0.45 * pow( physFunc::Reynolds_p( bin_Centers_1[iBin_1], rho_1, temperature_K, pressure_Pa ), 0.5 ) * pow( physFunc::Schmidt_p( bin_Centers_2[iBin_2], temperature_K, pressure_Pa ), 1.0 / double(3.0) ); + K_DE[iBin_1][iBin_2] = K_Brow[iBin_1][iBin_2] * 0.45 * sqrt( reynolds_p_1 ) * cbrt( physFunc::Schmidt_p( bin_Centers_2[iBin_2], temperature_K, pressure_Pa )); } } @@ -251,7 +254,7 @@ namespace AIM /* Initialize TI kernel */ for ( unsigned int iBin = 0; iBin < bin_Centers.size(); iBin++ ) { - K_TI[iBin] = physConst::PI * pow ( physConst::EPSILON , 0.75 ) / ( physConst::g * pow( physFunc::kinVisc( temperature_K, pressure_Pa ), 0.25 ) ) * ( bin_Centers[iBin] + bin_R ) * ( bin_Centers[iBin] + bin_R ) * std::abs( physFunc::vFall( bin_Centers[iBin], rho_1, temperature_K, pressure_Pa ) - physFunc::vFall( bin_R, rho_2, temperature_K, pressure_Pa ) ); + K_TI[iBin] = physConst::PI * pow ( physConst::EPSILON , 0.75 ) / ( physConst::g * sqrt(sqrt( physFunc::kinVisc( temperature_K, pressure_Pa ))) ) * ( bin_Centers[iBin] + bin_R ) * ( bin_Centers[iBin] + bin_R ) * std::abs( physFunc::vFall( bin_Centers[iBin], rho_1, temperature_K, pressure_Pa ) - physFunc::vFall( bin_R, rho_2, temperature_K, pressure_Pa ) ); } return K_TI; @@ -272,7 +275,7 @@ namespace AIM for ( unsigned int iBin_1 = 0; iBin_1 < bin_Centers_1.size(); iBin_1++ ) { K_TI.push_back( Vector_1D( bin_Centers_2.size() ) ); for ( unsigned int iBin_2 = 0; iBin_2 < bin_Centers_2.size(); iBin_2++ ) { - K_TI[iBin_1][iBin_2] = physConst::PI * pow ( physConst::EPSILON , 0.75 ) / ( physConst::g * pow( physFunc::kinVisc( temperature_K, pressure_Pa ), 0.25 ) ) * ( bin_Centers_1[iBin_1] + bin_Centers_2[iBin_2] ) * ( bin_Centers_1[iBin_1] + bin_Centers_2[iBin_2] ) * std::abs( physFunc::vFall( bin_Centers_1[iBin_1], rho_1, temperature_K, pressure_Pa ) - physFunc::vFall( bin_Centers_2[iBin_2], rho_2, temperature_K, pressure_Pa ) ); + K_TI[iBin_1][iBin_2] = physConst::PI * pow ( physConst::EPSILON , 0.75 ) / ( physConst::g * sqrt(sqrt( physFunc::kinVisc( temperature_K, pressure_Pa ))) ) * ( bin_Centers_1[iBin_1] + bin_Centers_2[iBin_2] ) * ( bin_Centers_1[iBin_1] + bin_Centers_2[iBin_2] ) * std::abs( physFunc::vFall( bin_Centers_1[iBin_1], rho_1, temperature_K, pressure_Pa ) - physFunc::vFall( bin_Centers_2[iBin_2], rho_2, temperature_K, pressure_Pa ) ); } } diff --git a/Code.v05-00/src/Core/Engine.cpp b/Code.v05-00/src/Core/Engine.cpp index 242f1cce..b7244441 100644 --- a/Code.v05-00/src/Core/Engine.cpp +++ b/Code.v05-00/src/Core/Engine.cpp @@ -241,7 +241,8 @@ Engine::Engine(std::string engineName, std::string engineFileName, H = -19.0 * (0.37318 * relHum_w / 100.0 * Pv) / (14.696 * delta - relHum_w / 100.0 * Pv); - EI_NOx *= exp(H) * pow(pow(delta, 1.02) / pow(theta, 3.3), 0.5); + double cruise_correction = pow(theta, 3.3) * pow(delta, -1.02); + EI_NOx *= exp(H) * sqrt(1 / cruise_correction); /* EI_NO in g(NO)/kg: * NOxtoNO is in mole of NO per mole of N @@ -299,7 +300,7 @@ Engine::Engine(std::string engineName, std::string engineFileName, EI_CO = pow(10.0, horzline); /* Cruise correction for CO */ - EI_CO *= pow(theta, 3.3) / pow(delta, 1.02); + EI_CO *= cruise_correction; /** Computing engine CO emission index **/ line1 = (log10(LTO_HC[1]) - log10(LTO_HC[0])) / @@ -340,7 +341,7 @@ Engine::Engine(std::string engineName, std::string engineFileName, EI_HC = pow(10.0, horzline); /* Cruise correction for HC */ - EI_HC *= pow(theta, 3.3) / pow(delta, 1.02); + EI_HC *= cruise_correction; /** Computing engine Soot emission index **/ EI_Soot = 0.02; /* [g/kg fuel] */ diff --git a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp index c0439aca..25b4bca9 100644 --- a/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp +++ b/Code.v05-00/src/Core/LAGRIDPlumeModel.cpp @@ -280,7 +280,8 @@ SimStatus LAGRIDPlumeModel::runFullModel() { std::variant LAGRIDPlumeModel::runEPM() { // Calculate the number of emitted soot particles - const double volParticle = 4.0 / 3.0 * physConst::PI * pow( EI_.getSootRad(), 3.0 ); //EI_SootRad in m -> volume in m3 + double EI_getSootRad = EI_.getSootRad(); + const double volParticle = 4.0 / 3.0 * physConst::PI * EI_getSootRad * EI_getSootRad * EI_getSootRad; //EI_SootRad in m -> volume in m3 const double massParticle = volParticle * physConst::RHO_SOOT * 1.0E+03; //Gives mass of a particle in grams const double EI_icenum = EI_.getSoot() / massParticle; /* [#/kg_fuel] */ const double N0 = EI_icenum * aircraft_.fuel_per_dist(); @@ -366,8 +367,8 @@ void LAGRIDPlumeModel::initializeGrid(const EPM::Output &epmOut) { const double ycenter = -aircraft_.vortex().z_center(); for (UInt n = 0; n < iceAerosol_.getNBin(); n++) { - double EPM_nPart_bin = epmIceAer.binMoment(n) * epmOut.area; double logBinRatio = log(iceAerosol_.getBinEdges()[n+1] / iceAerosol_.getBinEdges()[n]); + double EPM_nPart_bin = epmIceAer.binMoment(n) * epmOut.area * logBinRatio; pdf_init.push_back( LAGRID::initVarToGridRectangular(EPM_nPart_bin, xEdges_, yEdges_, xcenter, ycenter, initWidth, initDepth, logBinRatio) ); } diff --git a/Code.v05-00/src/Core/LiquidAer.cpp b/Code.v05-00/src/Core/LiquidAer.cpp index 85252672..f55770e0 100644 --- a/Code.v05-00/src/Core/LiquidAer.cpp +++ b/Code.v05-00/src/Core/LiquidAer.cpp @@ -155,13 +155,13 @@ unsigned int STRAT_AER( const double temperature_K , const double pressure_P /* Calculate conversion factors for SLA */ /* Factor to convert volume (m^3/m^3 air) to * surface area density (cm^2/cm^3 air) */ - const double SLA_VA = 8.406E-08 * pow( 10.0, 12.0 * 0.751E+00 ); + const double SLA_VA = 8.406E-08 * 1028016298.1264735114; /* Factor to convert effective radius to liquid radius (unitless) */ - const double SLA_RR = exp( -0.173E+00 ); + const double SLA_RR = 0.84113761484462320121169; /* Factor to convert volume (m^3/m^3) to effective radius (m) */ - const double SLA_VR = 0.357E-06 * pow( 10.0, 12.0 * 0.249E+00 ); + const double SLA_VR = 0.357E-06 * 972.747223776965; /* Reaction prefactors */ double KHET_COMMON; @@ -382,7 +382,7 @@ unsigned int STRAT_AER( const double temperature_K , const double pressure_P VOL_TOT = VOL_NAT + VOL_ICE; KG_AER_BOX = KG_NAT + KG_ICE; RAD_AER_BOX = MIN_RAD; - NDENS_AER_BOX = pow( ( 3.0E+00 * VOL_TOT / ( 4.0E+00 * physConst::PI * MAX_NDENS )), 1.0E+00 / 3.0E+00 ); + NDENS_AER_BOX = cbrt( ( 3.0E+00 * VOL_TOT / ( 4.0E+00 * physConst::PI * MAX_NDENS ))); } /* Prevent div zero */ @@ -1182,8 +1182,8 @@ void TERNARY( const double TIN_K , const double PIN_Pa , const double H2OSUM * Mole fraction of H2SO4 in binary solution */ X_H2SO4_BIN = 1.00E+00 / ( 2.00 * (KS[2]+KS[3]/T_K) ) * \ - (-KS[0]-KS[1]/T_K - pow( (KS[0]+KS[1]/T_K) * (KS[0]+KS[1]/T_K) \ - - 4.00E+00 * (KS[2]+KS[3]/T_K) * (KS[4]+KS[5]/T_K + KS[6]*log(T_K)-log(PATMH2O)), 0.50 ) ); + (-KS[0]-KS[1]/T_K - sqrt( (KS[0]+KS[1]/T_K) * (KS[0]+KS[1]/T_K) \ + - 4.00E+00 * (KS[2]+KS[3]/T_K) * (KS[4]+KS[5]/T_K + KS[6]*log(T_K)-log(PATMH2O))) ); /* Molality (mol H2SO4/kg H2O) in binary solution */ M_H2SO4_BIN = 5.551E+01 * X_H2SO4_BIN / ( 1.0E+00 - X_H2SO4_BIN ); @@ -1196,8 +1196,8 @@ void TERNARY( const double TIN_K , const double PIN_Pa , const double H2OSUM + ( QS[6] + QS[7]*TR + QS[8]*TR*TR )*PR*PR \ + ( QS[9]*TR )*PR*PR*PR ); X_HNO3_BIN = 1.00E+00 / ( 2.00 * (KN[2]+KN[3]/T_K) ) * \ - (-KN[0]-KN[1]/T_K - pow( (KN[0]+KN[1]/T_K) * (KN[0]+KN[1]/T_K) \ - -4.0E+00 * (KN[2]+KN[3]/T_K) * (KN[4]+KN[5]/T_K + KN[6]*log(T_K) - log(PATMH2O)), 0.50 ) ); + (-KN[0]-KN[1]/T_K - sqrt( (KN[0]+KN[1]/T_K) * (KN[0]+KN[1]/T_K) \ + -4.0E+00 * (KN[2]+KN[3]/T_K) * (KN[4]+KN[5]/T_K + KN[6]*log(T_K) - log(PATMH2O))) ); /* Molality */ M_HNO3_BIN = 5.551E+01 * X_HNO3_BIN / ( 1.0E+00 - X_HNO3_BIN ); H_HNO3_BIN = exp( QN[0] + QN[1]*TR*TR \ @@ -1371,12 +1371,12 @@ double CARSLAW_DENSITY( const double CS, const double CN, const double T_K ) /* CARSLAW_DENSITY begins here! */ - DENSS = 1.000E+03 + 1.2364E+02*CS - 5.600E-04*CS*T_K*T_K - 2.954E+01*pow( CS, 1.5 ) \ - + 1.814E-04*pow( CS, 1.5 )*T_K*T_K + 2.343E+00*CS*CS - 1.487E-03*CS*CS*T_K \ + DENSS = 1.000E+03 + 1.2364E+02*CS - 5.600E-04*CS*T_K*T_K - 2.954E+01*CS*sqrt(CS) \ + + 1.814E-04*CS*sqrt(CS)*T_K*T_K + 2.343E+00*CS*CS - 1.487E-03*CS*CS*T_K \ - 1.324E-05*CS*CS*T_K*T_K; - DENSN = 1.000E+03 + 8.5107E+01*CN - 5.043E-04*CN*T_K*T_K - 1.896E+01*pow( CN, 1.5 ) \ - + 1.427E-04*pow( CN, 1.5 )*T_K*T_K + 1.458E+00*CN*CN - 1.198E-03*CN*CN*T_K \ + DENSN = 1.000E+03 + 8.5107E+01*CN - 5.043E-04*CN*T_K*T_K - 1.896E+01*CN*sqrt(CN) \ + + 1.427E-04*CN*sqrt(CN)*T_K*T_K + 1.458E+00*CN*CN - 1.198E-03*CN*CN*T_K \ - 9.703E-06*CN*CN*T_K*T_K; return 1.00E+00/( (1.00E+00/DENSS*CS/(CS+CN) + 1.00E+00/DENSN*CN/(CS+CN)) ); diff --git a/Code.v05-00/src/Core/Vortex.cpp b/Code.v05-00/src/Core/Vortex.cpp index d50000b7..9fb29b77 100644 --- a/Code.v05-00/src/Core/Vortex.cpp +++ b/Code.v05-00/src/Core/Vortex.cpp @@ -74,7 +74,7 @@ Vortex::Vortex( double RHi_PC, double temperature_K, double pressure_Pa, \ N_BV_ = N_BV; /* Maximum downwash displacement, Eq. 5 in Lottermosser and Unterstrasser (2025) */ - z_desc_ = pow( 8.0 * gamma_ / ( physConst::PI * N_BV_ ), 0.5 ); + z_desc_ = sqrt( 8.0 * gamma_ / ( physConst::PI * N_BV_ ) ); /* Height an air parcel has to descend until it is no longer supersaturated */ s_ = RHi_PC / 100 - 1; // RHi in % -> excess supersaturation ratio, See S2 in U2016 @@ -85,8 +85,8 @@ Vortex::Vortex( double RHi_PC, double temperature_K, double pressure_Pa, \ r_p_ref_ = 1.5 + 0.314 * wingspan_ref; /* [m], from Eq. A6 in U2016 */ /* Plume area before vortex breakup*/ - plume_area_0_ = 2 * physConst::PI * pow(r_p_, 2); /* [m2], see Appendix 2 in LU2025 */ - plume_area_0_ref_ = 2 * physConst::PI * pow(r_p_ref_, 2); /* [m2], see Appendix 2 in LU2025 */ + plume_area_0_ = 2 * physConst::PI * r_p_ * r_p_; /* [m2], see Appendix 2 in LU2025 */ + plume_area_0_ref_ = 2 * physConst::PI * r_p_ref_ * r_p_ref_; /* [m2], see Appendix 2 in LU2025 */ /* Temperature - 205 K*/ T_205_ = temperature_K - 205.0; /* [K], from Eq. A3 in LU2025*/ diff --git a/Code.v05-00/src/EPM/Models/Original/Integrate.cpp b/Code.v05-00/src/EPM/Models/Original/Integrate.cpp index b95ec4b1..fc01be8e 100644 --- a/Code.v05-00/src/EPM/Models/Original/Integrate.cpp +++ b/Code.v05-00/src/EPM/Models/Original/Integrate.cpp @@ -96,7 +96,7 @@ namespace EPM::Models * radii * and volume ratio between two adjacent bins */ const UInt SO4_NBIN = static_cast( - std::floor(1 + log(pow(LA_R_HIG / LA_R_LOW, 3.0)) / log(LA_VRAT))); + std::floor(1 + 3.0*log(LA_R_HIG / LA_R_LOW) / log(LA_VRAT))); if (LA_VRAT <= 1.0) { std::cout << "\nVolume ratio of consecutive bins for SO4 has to be " @@ -110,7 +110,7 @@ namespace EPM::Models Vector_1D SO4_vJ( SO4_NBIN , 0.0 ); // bin centers, volume /* Adjacent bin radius ratio */ - const double LA_RRAT = pow( LA_VRAT, 1.0 / double(3.0) ); + const double LA_RRAT = cbrt(LA_VRAT); /* Initialize bin center and edge radii, as well as volume */ SO4_rE[0] = LA_R_LOW; @@ -125,10 +125,10 @@ namespace EPM::Models /* Number of ice size distribution bins based on specified min and max radii * * and volume ratio between two adjacent bins */ const UInt Ice_NBIN = static_cast( - std::floor(1 + log(pow(PA_R_HIG / PA_R_LOW, 3.0)) / log(PA_VRAT))); + std::floor(1 + 3.0*log(PA_R_HIG / PA_R_LOW) / log(PA_VRAT))); /* Adjacent bin radius ratio */ - const double PA_RRAT = pow( PA_VRAT, 1.0 / double(3.0) ); + const double PA_RRAT = cbrt( PA_VRAT ); /* Ice bin center and edge radii */ Vector_1D Ice_rJ( Ice_NBIN , 0.0 ); @@ -208,8 +208,9 @@ namespace EPM::Models double log10_timeInitial = log10(timeInitial); double timeFinal = 2.0E+03; double log10_timeFinal = log10(timeFinal); + double denominator = 1/double( nTime - 1.0 ); for ( iTime = 0; iTime < nTime; iTime++ ) { - timeArray[iTime] = pow( 10.0, log10_timeInitial + iTime * ( log10_timeFinal - log10_timeInitial ) / double( nTime - 1.0 ) ); + timeArray[iTime] = pow( 10.0, log10_timeInitial + iTime * ( log10_timeFinal - log10_timeInitial ) * denominator ); } UInt iTime_3mins; diff --git a/Code.v05-00/src/EPM/Solution.cpp b/Code.v05-00/src/EPM/Solution.cpp index 2304c749..b0cde6ec 100644 --- a/Code.v05-00/src/EPM/Solution.cpp +++ b/Code.v05-00/src/EPM/Solution.cpp @@ -100,14 +100,15 @@ void Solution::Initialize(std::string fileName, VectorUtils::set_value(sootDens, aer_Value[0][0]); VectorUtils::set_value(sootRadi, aer_Value[0][1]); VectorUtils::set_value(sootArea, 4.0 / double(3.0) * PI * aer_Value[0][0] * aer_Value[0][1] * aer_Value[0][1] * aer_Value[0][1]); - - nBin_LA = std::floor(1 + log(pow((LA_R_HIG/LA_R_LOW), 3.0)) / log(LA_VRAT)); + + double la_r_hig_low = LA_R_HIG/LA_R_LOW; + nBin_LA = std::floor(1 + log(la_r_hig_low * la_r_hig_low * la_r_hig_low) / log(LA_VRAT)); Vector_1D LA_rE( nBin_LA + 1, 0.0 ); /* Bin edges in m */ Vector_1D LA_rJ( nBin_LA , 0.0 ); /* Bin center radius in m */ Vector_1D LA_vJ( nBin_LA , 0.0 ); /* Bin volume centers in m^3 */ - const double LA_RRAT = pow( LA_VRAT, 1.0 / double(3.0) ); + const double LA_RRAT = cbrt( LA_VRAT ); LA_rE[0] = LA_R_LOW; for ( UInt iBin_LA = 1; iBin_LA < nBin_LA + 1; iBin_LA++ ) LA_rE[iBin_LA] = LA_rE[iBin_LA-1] * LA_RRAT; /* [m] */ @@ -142,14 +143,14 @@ void Solution::Initialize(std::string fileName, nBin_LA = 2; //dumb hardcoded Grid_Aerosol default constructor } - - nBin_PA = std::floor( 1 + log( pow( (PA_R_HIG/PA_R_LOW), 3.0 ) ) / log( PA_VRAT ) ); + double pa_r_hig_low = PA_R_HIG/PA_R_LOW; + nBin_PA = std::floor( 1 + log( pa_r_hig_low * pa_r_hig_low* pa_r_hig_low ) / log( PA_VRAT ) ); Vector_1D PA_rE( nBin_PA + 1, 0.0 ); /* Bin edges in m */ Vector_1D PA_rJ( nBin_PA , 0.0 ); /* Bin center radius in m */ Vector_1D PA_vJ( nBin_PA , 0.0 ); /* Bin volume centers in m^3 */ - const double PA_RRAT = pow( PA_VRAT, 1.0 / double(3.0) ); + const double PA_RRAT = cbrt( PA_VRAT ); PA_rE[0] = PA_R_LOW; for ( UInt iBin_PA = 1; iBin_PA < nBin_PA + 1; iBin_PA++ ) PA_rE[iBin_PA] = PA_rE[iBin_PA-1] * PA_RRAT; /* [m] */ diff --git a/Code.v05-00/src/KPP/KPP_HetRates.cpp b/Code.v05-00/src/KPP/KPP_HetRates.cpp index b4ab6285..641d567c 100644 --- a/Code.v05-00/src/KPP/KPP_HetRates.cpp +++ b/Code.v05-00/src/KPP/KPP_HetRates.cpp @@ -582,7 +582,7 @@ double HETNO3( const double A, const double B, const double AREA[NAERO], const d ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -693,7 +693,7 @@ double HETHO2( const double A, const double B, const double AREA[NAERO], const d ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -754,7 +754,7 @@ double HETHBr( const double A, const double B, const double KHETI_SLA[11], const ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -831,7 +831,7 @@ double HETN2O5( const double A, const double B, const double KHETI_SLA[11], cons if ( N == 1 ) { ADJUSTEDRATE = AREA[N] * XSTKCF; } else { - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -899,7 +899,7 @@ double HETBrNO3( const double A, const double B, const double KHETI_SLA[11], con ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -966,7 +966,7 @@ double HETHOBr( const double A, const double B, const double KHETI_SLA[11], cons ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -1465,7 +1465,7 @@ double HETN2O5_PSC( const double A, const double B, const double KHETI_SLA[11], ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -1538,7 +1538,7 @@ double HETClNO3_PSC1( const double A, const double B, const double KHETI_SLA[11] ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -1611,7 +1611,7 @@ double HETClNO3_PSC2( const double A, const double B, const double KHETI_SLA[11] ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -1684,7 +1684,7 @@ double HETClNO3_PSC3( const double A, const double B, const double KHETI_SLA[11] ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -1757,7 +1757,7 @@ double HETBrNO3_PSC( const double A, const double B, const double KHETI_SLA[11], ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -1830,7 +1830,7 @@ double HETHOCl_PSC1( const double A, const double B, const double KHETI_SLA[11], ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -1903,7 +1903,7 @@ double HETHOCl_PSC2( const double A, const double B, const double KHETI_SLA[11], ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { @@ -1976,7 +1976,7 @@ double HETHOBr_PSC( const double A, const double B, const double KHETI_SLA[11], ADJUSTEDRATE = AREA[N] * XSTKCF; } else { /* Reaction rate for surface of aerosol */ - ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, pow( TEMP, 0.5 ), pow( A, 0.5 ) ); + ADJUSTEDRATE = ARSL1K( AREA[N], RADI[N], AIRDENS, XSTKCF, sqrt( TEMP ), sqrt( A ) ); } if ( ( DO_EDUCT ) && ( N < 2 ) ) { diff --git a/Code.v05-00/src/KPP/KPP_Integrator_ADJ.cpp b/Code.v05-00/src/KPP/KPP_Integrator_ADJ.cpp index e5bd4ecd..5855b332 100644 --- a/Code.v05-00/src/KPP/KPP_Integrator_ADJ.cpp +++ b/Code.v05-00/src/KPP/KPP_Integrator_ADJ.cpp @@ -1837,7 +1837,7 @@ double ros_ErrorNorm ( double Y[], double Ynew[], double Yerr[], else Scale = AbsTol[0]+RelTol[0]*Ymax; - Err = Err+pow((Yerr[i]/Scale),2); + Err = Err+((Yerr[i]/Scale) * (Yerr[i]/Scale)); } Err = SQRT(Err/NVAR); @@ -2057,9 +2057,10 @@ void ros_Hermite3( double a, double b, double T, double Ya[], Tau = T - a; WCOPY(NVAR,&C[3][0],1,Y,1); - WSCAL(NVAR,pow(Tau,3),Y,1); - for(j=2; j>=0; j--) - WAXPY(NVAR,pow(Tau,(j-1)),&C[j][0],1,Y,1); + WSCAL(NVAR,Tau*Tau*Tau,Y,1); + WAXPY(NVAR,Tau,&C[2][0],1,Y,1); + WAXPY(NVAR,1,&C[1][0],1,Y,1); + WAXPY(NVAR,1/Tau,&C[0][0],1,Y,1); } /* End of ros_Hermite3 */ diff --git a/Code.v05-00/src/KPP/KPP_LinearAlgebra.cpp b/Code.v05-00/src/KPP/KPP_LinearAlgebra.cpp index 263aef1e..c5f5fdcd 100644 --- a/Code.v05-00/src/KPP/KPP_LinearAlgebra.cpp +++ b/Code.v05-00/src/KPP/KPP_LinearAlgebra.cpp @@ -1844,7 +1844,10 @@ double WLAMCH( char C ) if (First) { First = 0; - Eps = pow(HALF,16); + double HALF2 = HALF * HALF; + double HALF4 = HALF2 * HALF2; + double HALF8 = HALF4 * HALF4; + Eps = HALF8 * HALF8; for ( i = 17; i <= 80; i++ ) { Eps = Eps*HALF; Suma = WLAMCH_ADD(ONE,Eps); diff --git a/Code.v05-00/src/KPP/KPP_Main_ADJ.cpp b/Code.v05-00/src/KPP/KPP_Main_ADJ.cpp index 653ae2ec..45948b2c 100644 --- a/Code.v05-00/src/KPP/KPP_Main_ADJ.cpp +++ b/Code.v05-00/src/KPP/KPP_Main_ADJ.cpp @@ -539,9 +539,11 @@ int KPP_Main_ADJ( const double finalPlume[], const double initBackg[], \ METRIC = 0.0; for ( i = 0; i < NOPT - 1; i++ ) { if ( i < 1 ) - METRIC += pow( ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ), 2); + METRIC += ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ) + * ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ); else - METRIC += pow( ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ), 2); + METRIC += ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ) + * ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ); } if ( VERBOSE ) @@ -666,9 +668,11 @@ int KPP_Main_ADJ( const double finalPlume[], const double initBackg[], \ METRIC = 0.0; for ( i = 0; i < NOPT - 1; i++ ) { if ( i < 1 ) //NOx - METRIC += pow( ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ), 2); + METRIC += ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ) + * ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ); else - METRIC += pow( ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ), 2); + METRIC += ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ) + * ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ); } if ( VERBOSE ) { @@ -859,9 +863,11 @@ int KPP_Main_ADJ( const double finalPlume[], const double initBackg[], \ METRIC_ABS_MIN = 0.0; for ( i = 0; i < NOPT - 1; i++ ) { if ( i < 1 ) - METRIC_ABS_MIN += pow( ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ), 2); + METRIC_ABS_MIN += ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ) + * ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ); else - METRIC_ABS_MIN += pow( ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ), 2); + METRIC_ABS_MIN += ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ) + * ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ); } if ( VERBOSE ) @@ -921,9 +927,11 @@ int KPP_Main_ADJ( const double finalPlume[], const double initBackg[], \ METRIC = 0.0; for ( i = 0; i < NOPT - 1; i++ ) { if ( i < 1 ) - METRIC += pow( ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ), 2); + METRIC += ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ) + * ( DIAG[0] * ( VAR_RUN[ind_OPT[0]] + VAR_RUN[ind_OPT[1]] - finalPlume[ind_OPT[0]] - finalPlume[ind_OPT[1]] ) ); else - METRIC += pow( ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ), 2); + METRIC += ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ) + * ( DIAG[i+1] * ( VAR_RUN[ind_OPT[i+1]] - finalPlume[ind_OPT[i+1]] ) ); } if ( VERBOSE ) diff --git a/Code.v05-00/src/KPP/KPP_Rates.cpp b/Code.v05-00/src/KPP/KPP_Rates.cpp index afde4474..ffc360d6 100644 --- a/Code.v05-00/src/KPP/KPP_Rates.cpp +++ b/Code.v05-00/src/KPP/KPP_Rates.cpp @@ -315,14 +315,16 @@ double GC_OHCO( float A0, float B0, float C0, double PRESS, double AIRDENS, doub K0 = GCARR( (double)A0, (double)B0, (double)C0, TEMP ); K0 = K0 * (1.0E+00 + 0.6E+00*9.871E+07*PRESS); - KLO1 = 5.9E-33*pow( (300.0/TEMP), (1.4E+00) ); - KHI1 = 1.1E-12*pow( (300.0/TEMP), (-1.3E+00) ); + double temp_300 = 300.0/TEMP; + + KLO1 = 5.9E-33*pow( (temp_300), (1.4E+00) ); + KHI1 = 1.1E-12*pow( (temp_300), (-1.3E+00) ); XYRAT1 = KLO1*AIRDENS/KHI1; BLOG1 = log10(XYRAT1); FEXP1 = 1.0E+00/(1.0E+00 + BLOG1 * BLOG1); KCO1 = KLO1 * AIRDENS * pow( 0.6, (FEXP1) ) / (1.0E+00 + XYRAT1); - KLO2 = 1.5E-13*pow( (300.0/TEMP), (-0.6E+00) ); - KHI2 = 2.1E+09*pow( (300.0/TEMP), (-6.1E+00) ); + KLO2 = 1.5E-13*pow( (temp_300), (-0.6E+00) ); + KHI2 = 2.1E+09*pow( (temp_300), (-6.1E+00) ); XYRAT2 = KLO2*AIRDENS / KHI2; BLOG2 = log10(XYRAT2); FEXP2 = 1.0E+00 / (1.0E+00 + BLOG2*BLOG2); diff --git a/Code.v05-00/src/LAGRID/RemappingFunctions.cpp b/Code.v05-00/src/LAGRID/RemappingFunctions.cpp index 67a59575..f1b3564a 100644 --- a/Code.v05-00/src/LAGRID/RemappingFunctions.cpp +++ b/Code.v05-00/src/LAGRID/RemappingFunctions.cpp @@ -249,7 +249,7 @@ namespace LAGRID { double sigmaX, double sigmaY, double logBinRatio ) { auto gaussianFunc = [x0, y0, sigmaX, sigmaY] (double x, double y) -> double { - return exp(- (pow(x - x0, 2.0) / (2.0 * sigmaX * sigmaX) + pow(y - y0, 2.0) / (2.0 * sigmaY * sigmaY))); + return exp(- ( ( (x - x0)*(x - x0) ) / (2.0 * sigmaX * sigmaX) + ( (y - y0)*(y - y0) ) / (2.0 * sigmaY * sigmaY) ) ); }; return initVarToGrid(mass, xEdges, yEdges, gaussianFunc, logBinRatio); } @@ -260,7 +260,7 @@ namespace LAGRID { double sigmaX = width / 8; double omega = 2 * (physConst::PI / depth); auto func = [omega, x0, y0, sigmaX, depth] (double x, double y) -> double { - return exp(- (pow(x - x0, 2.0) / (2.0 * sigmaX * sigmaX) )) * ( std::abs(y - y0) < depth/2 ) * std::abs( std::sin(omega * (y - y0)) ); + return exp(- ( ( (x - x0)*(x - x0) ) / (2.0 * sigmaX * sigmaX) )) * ( std::abs(y - y0) < depth/2 ) * std::abs( std::sin(omega * (y - y0)) ); }; return initVarToGrid(mass, xEdges, yEdges, func, logBinRatio); } diff --git a/Code.v05-00/src/Util/PhysFunction.cpp b/Code.v05-00/src/Util/PhysFunction.cpp index 458932e9..cb306e29 100644 --- a/Code.v05-00/src/Util/PhysFunction.cpp +++ b/Code.v05-00/src/Util/PhysFunction.cpp @@ -149,7 +149,7 @@ namespace physFunc * OUTPUT PARAMETERS: * - double :: Dynamic viscosity in kg/m/s */ - return 1.8325E-05 * ( 416.16 / ( T + 120.0 ) ) * pow( T / 296.16, 1.5 ); + return 1.8325E-05 * ( 416.16 / ( T + 120.0 ) ) * ( T / 296.16 ) * sqrt ( T / 296.16); } /* End of dynVisc */ @@ -344,7 +344,7 @@ namespace physFunc double l = lambda_p( r, m, T, P ); - return ( pow( 2.0 * r + l , 3.0 ) - pow( 4.0 * r * r + l * l , 1.5 ) ) / ( 6.0 * r * l ) - 2.0 * r; + return ( (2.0 * r + l)*(2.0 * r + l)*(2.0 * r + l) - (4.0 * r * r + l * l) * sqrt(4.0 * r * r + l * l) ) / ( 6.0 * r * l ) - 2.0 * r; } /* End of delta_p */ @@ -468,7 +468,8 @@ namespace physFunc /* r_M = r_1 and r_m = r_2 */ s = Stokes_p( r_2, rho_2, r_1, rho_1, T, P ); if ( s > 1.214 ) { - E_V = pow(1 + 0.75*log(2*s)/(s - 1.214), -2); + double pow_var = 1 + 0.75*log(2*s)/(s - 1.214); + E_V = 1 / (pow_var * pow_var); } else { E_V = 0.0; } @@ -478,7 +479,8 @@ namespace physFunc /* r_M = r_2 and r_m = r_1 */ s = Stokes_p( r_1, rho_1, r_2, rho_2, T, P ); if ( s > 1.214 ) { - E_V = pow(1 + 0.75*log(2*s)/(s - 1.214), -2); + double pow_var = 1 + 0.75*log(2*s)/(s - 1.214); + E_V = 1 / (pow_var * pow_var); } else { E_V = 0.0; } From 310a0496394245fe01c9cfcfe5f6c82a79ed8a2e Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Mon, 16 Mar 2026 17:14:41 +0000 Subject: [PATCH 4/5] Fix parallelisation bug (#9) Co-authored-by: Coco Yeung --- Code.v05-00/include/AIM/Aerosol.hpp | 12 ++ Code.v05-00/src/AIM/Aerosol.cpp | 225 +++++++++++++++------------- 2 files changed, 129 insertions(+), 108 deletions(-) diff --git a/Code.v05-00/include/AIM/Aerosol.hpp b/Code.v05-00/include/AIM/Aerosol.hpp index a1bc0299..b03ef4fa 100644 --- a/Code.v05-00/include/AIM/Aerosol.hpp +++ b/Code.v05-00/include/AIM/Aerosol.hpp @@ -69,6 +69,8 @@ class AIM::Aerosol return pow(bin_Centers[iBin], n) * pdf[iBin]; } //NOTE: This gives the moment in [- / cm3]. You need to multiply by factors to get it in [ / m] or something. + template + double Moment() const; double Moment( UInt n = 0 ) const; double Radius( ) const; double EffRadius( ) const; @@ -83,6 +85,7 @@ class AIM::Aerosol inline const Vector_1D& getBinVCenters() const { return bin_VCenters; }; inline const Vector_1D& getBinEdges() const { return bin_Edges; }; inline const Vector_1D& getBinSizes() const { return bin_Sizes; }; + inline const Vector_1D& getLogBinEdges() const { return log_Bin_Edges; }; inline UInt getNBin() const { return nBin; }; inline const Vector_1D& getPDF() const { return pdf; }; @@ -92,6 +95,7 @@ class AIM::Aerosol Vector_1D bin_VCenters; Vector_1D bin_Edges; Vector_1D bin_Sizes; + Vector_1D log_Bin_Edges; UInt nBin; Vector_1D pdf; // represents dN [-/cm3] / d(ln r) }; @@ -128,8 +132,14 @@ class AIM::Grid_Aerosol inline void updateNy(int ny_new) { Ny = ny_new; }; /* Moments */ + template + Vector_2D Moment() const; Vector_2D Moment( UInt n ) const; + template + double Moment( const Vector_1D& PDF ) const; double Moment( UInt n, const Vector_1D& PDF ) const; + template + double Moment( UInt iNx, UInt jNy ) const; double Moment( UInt n, UInt iNx, UInt jNy ) const; /* Extra utils */ @@ -181,6 +191,7 @@ class AIM::Grid_Aerosol inline const Vector_3D& getBinVCenters() const { return bin_VCenters; }; inline const Vector_1D& getBinEdges() const { return bin_Edges; }; inline const Vector_1D& getBinSizes() const { return bin_Sizes; }; + inline const Vector_1D& getLogBinEdges() const { return log_Bin_Edges; }; inline UInt getNBin() const { return nBin; }; inline const Vector_3D& getPDF() const { return pdf; }; inline Vector_3D& getPDF() { return pdf; }; @@ -193,6 +204,7 @@ class AIM::Grid_Aerosol Vector_1D bin_Edges; Vector_1D bin_VEdges; Vector_1D bin_Sizes; + Vector_1D log_Bin_Edges; UInt nBin; unsigned int Nx, Ny; Vector_3D pdf; //Everything with the pdf is implicitly in [ / cm3] diff --git a/Code.v05-00/src/AIM/Aerosol.cpp b/Code.v05-00/src/AIM/Aerosol.cpp index d6aa6706..8368985f 100644 --- a/Code.v05-00/src/AIM/Aerosol.cpp +++ b/Code.v05-00/src/AIM/Aerosol.cpp @@ -31,6 +31,7 @@ namespace AIM bin_VCenters(bin_Centers_.size()), bin_Edges(bin_Edges_), bin_Sizes(bin_Centers_.size()), + log_Bin_Edges(bin_Edges_.size()), nBin(bin_Centers_.size()), pdf(bin_Centers_.size()) { /* In the following cases, the aerosol pdf represents the following quantity: dn/d(ln(r)) @@ -49,7 +50,9 @@ namespace AIM for (UInt iBin = 0; iBin < nBin; iBin++) { bin_Sizes[iBin] = bin_Edges[iBin + 1] - bin_Edges[iBin]; bin_VCenters[iBin] = 4.0 / 3.0 * PI * (bin_Edges_[iBin] * bin_Edges_[iBin] * bin_Edges_[iBin] + bin_Edges_[iBin + 1] * bin_Edges_[iBin + 1] * bin_Edges_[iBin + 1]) * 0.5; + log_Bin_Edges[iBin] = log(bin_Edges[iBin]); } + log_Bin_Edges[nBin] = log(bin_Edges[nBin]); /* Check mean and standard deviation */ if (mu <= 0) { @@ -115,12 +118,15 @@ namespace AIM bin_VCenters(bin_Centers.size()), bin_Edges(bin_Edges), bin_Sizes(bin_Centers.size()), + log_Bin_Edges(bin_Edges.size()), nBin(bin_Centers.size()), pdf(pdf) { for (UInt iBin = 0; iBin < nBin; iBin++) { bin_Sizes[iBin] = bin_Edges[iBin + 1] - bin_Edges[iBin]; bin_VCenters[iBin] = 4.0 / 3.0 * PI * (bin_Edges[iBin] * bin_Edges[iBin] * bin_Edges[iBin] + bin_Edges[iBin + 1] * bin_Edges[iBin + 1] * bin_Edges[iBin + 1]) * 0.5; + log_Bin_Edges[iBin] = log(bin_Edges[iBin]); } + log_Bin_Edges[nBin] = log(bin_Edges[nBin]); } @@ -302,41 +308,43 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { } /* End of Aerosol::Coagulate */ - double Aerosol::Moment(UInt n) const + template + double Aerosol::Moment() const { double moment = 0; - double pow_value; - double log_current = log(bin_Edges[0]); - double log_next; for (UInt iBin = 0; iBin < nBin; iBin++) { - log_next = log(bin_Edges[iBin + 1]); - switch(n) { - case 0: - pow_value = 1; - break; - case 1: - pow_value = bin_Centers[iBin]; - break; - case 2: - pow_value = bin_Centers[iBin] * bin_Centers[iBin]; - break; - case 3: - pow_value = bin_Centers[iBin] * bin_Centers[iBin] * bin_Centers[iBin]; - break; - default: - pow_value = pow(bin_Centers[iBin], n); // fallback for safety - break; - } - moment += (log_next - log_current) * pow_value * pdf[iBin]; - log_current = log_next; // reuse for next iteration + double pow_value; + if constexpr (N == 0) + pow_value = 1; + else if constexpr (N == 1) + pow_value = bin_Centers[iBin]; + else if constexpr (N == 2) + pow_value = bin_Centers[iBin] * bin_Centers[iBin]; + else if constexpr (N == 3) + pow_value = bin_Centers[iBin] * bin_Centers[iBin] * bin_Centers[iBin]; + else + pow_value = pow(bin_Centers[iBin], N); // fallback for safety + + moment += (log_Bin_Edges[iBin+1] - log_Bin_Edges[iBin]) * pow_value * pdf[iBin]; } return moment; - } /* End of Aerosol::Moment */ + } + double Aerosol::Moment(UInt n) const + { + switch(n) { + case 0: return Moment<0>(); + case 1: return Moment<1>(); + case 2: return Moment<2>(); + case 3: return Moment<3>(); + default: return Moment<4>(); + } + } + /* End of Aerosol::Moment */ double Aerosol::Radius() const { @@ -426,6 +434,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { bin_Edges(bin_Edges_), bin_VEdges(bin_Centers_.size() + 1), bin_Sizes(bin_Centers_.size()), + log_Bin_Edges(bin_Edges_.size()), nBin(bin_Centers.size()), Nx(Nx_), Ny(Ny_) { /* In the following cases, the aerosol pdf represents the following quantity: dn/d(ln(r)) @@ -444,7 +453,9 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { { bin_VEdges[iBin] = 4.0 / 3.0 * PI * bin_Edges[iBin] * bin_Edges[iBin] * bin_Edges[iBin]; bin_Sizes[iBin] = bin_Edges[iBin + 1] - bin_Edges[iBin]; + log_Bin_Edges[iBin] = log(bin_Edges[iBin]); } + log_Bin_Edges[nBin] = log(bin_Edges[nBin]); bin_VEdges[nBin] = 4.0 / 3.0 * PI * bin_Edges[nBin] * bin_Edges[nBin] * bin_Edges[nBin]; bin_VCenters.resize(nBin, Vector_2D(Ny, Vector_1D(Nx, 0.0E+00))); @@ -661,7 +672,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { for (jBin = 0; jBin < nBin; jBin++) { - nPart = pdf[jBin][jNy][iNx] * log(bin_Edges[jBin + 1] / bin_Edges[jBin]); + nPart = pdf[jBin][jNy][iNx] * (log_Bin_Edges[jBin + 1] - log_Bin_Edges[jBin]); if (jBin <= iBin) { @@ -974,7 +985,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { { // Bin is not empty. Compute particle volume, and clip it between min and max volume allowed. bin_VCenters[iBin][y_index][x_index] = std::max(std::min(iceVol_ / icePart_, bin_VEdges[iBin + 1]), bin_VEdges[iBin]); - pdf[iBin][y_index][x_index] = icePart_ / (log(bin_Edges[iBin + 1] / bin_Edges[iBin])); + pdf[iBin][y_index][x_index] = icePart_ / (log_Bin_Edges[iBin + 1] - log_Bin_Edges[iBin]); } else { @@ -1125,7 +1136,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { { //Must resize the bin_VCenters to avoid indexing errors bin_VCenters[iBin] = Vector_2D(Ny, Vector_1D(Nx)); - double ratio = log(bin_Edges[iBin + 1] / bin_Edges[iBin]); + double ratio = log_Bin_Edges[iBin + 1] - log_Bin_Edges[iBin]; for (UInt jNy = 0; jNy < Ny; jNy++) { for (UInt iNx = 0; iNx < Nx; iNx++) @@ -1145,7 +1156,8 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { } /* End of Grid_Aerosol::UpdateCenters */ - Vector_2D Grid_Aerosol::Moment(UInt n) const + template + Vector_2D Grid_Aerosol::Moment() const { UInt jNy = 0; @@ -1155,44 +1167,42 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { Vector_2D moment(Ny, Vector_1D(Nx, 0.0E+00)); const double FACTOR = 3.0 / double(4.0 * PI); - double log_current = log(bin_Edges[0]); - double log_next; - double pow_value; - #pragma omp parallel for default(shared) private(iNx, jNy, iBin) \ schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++) { - log_next = log(bin_Edges[iBin + 1]); + double pow_value; for (jNy = 0; jNy < Ny; jNy++) { for (iNx = 0; iNx < Nx; iNx++) { - switch(n) { - case 0: - pow_value = 1; - break; - case 1: - pow_value = cbrt(FACTOR * bin_VCenters[iBin][jNy][iNx]); - break; - case 2: - pow_value = cbrt(FACTOR * FACTOR * bin_VCenters[iBin][jNy][iNx] * bin_VCenters[iBin][jNy][iNx]); - break; - case 3: - pow_value = FACTOR * bin_VCenters[iBin][jNy][iNx]; - break; - default: - pow_value = pow(FACTOR * bin_VCenters[iBin][jNy][iNx], n / 3.0); // fallback for safety - break; - } - moment[jNy][iNx] += (log_next - log_current) * pow_value * pdf[iBin][jNy][iNx]; + if constexpr (N == 0) + pow_value = 1; + else if constexpr (N == 1) + pow_value = cbrt(FACTOR * bin_VCenters[iBin][jNy][iNx]); + else if constexpr (N == 2) + pow_value = cbrt(FACTOR * FACTOR * bin_VCenters[iBin][jNy][iNx] * bin_VCenters[iBin][jNy][iNx]); + else if constexpr (N == 3) + pow_value = FACTOR * bin_VCenters[iBin][jNy][iNx]; + else + pow_value = pow(FACTOR * bin_VCenters[iBin][jNy][iNx], N / 3.0); // fallback for safety + moment[jNy][iNx] += (log_Bin_Edges[iBin+1] - log_Bin_Edges[iBin]) * pow_value * pdf[iBin][jNy][iNx]; } } - log_current = log_next; // reuse for next iteration } return moment; + } + Vector_2D Grid_Aerosol::Moment(UInt n) const + { + switch(n) { + case 0: return Moment<0>(); + case 1: return Moment<1>(); + case 2: return Moment<2>(); + case 3: return Moment<3>(); + default: return Moment<4>(); + } } /* End of Grid_Aerosol::Moment */ Vector_3D Grid_Aerosol::Number() const @@ -1209,7 +1219,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++) { - ratio = log(bin_Edges[iBin + 1] / bin_Edges[iBin]); + ratio = log_Bin_Edges[iBin + 1] - log_Bin_Edges[iBin]; for (jNy = 0; jNy < Ny; jNy++) { for (iNx = 0; iNx < Nx; iNx++) @@ -1293,7 +1303,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++) { - ratio = log(bin_Edges[iBin + 1] / bin_Edges[iBin]); + ratio = log_Bin_Edges[iBin + 1] - log_Bin_Edges[iBin]; for (jNy = 0; jNy < Ny; jNy++) { for (iNx = 0; iNx < Nx; iNx++) @@ -1662,52 +1672,51 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { } /* End of Grid_Aerosol::StdDev */ - double Grid_Aerosol::Moment(UInt n, const Vector_1D& PDF) const + template + double Grid_Aerosol::Moment(const Vector_1D& PDF) const { UInt iBin = 0; double moment = 0.0E+00; - double log_current = log(bin_Edges[0]); - double log_next; - double pow_value; - #pragma omp parallel for default(shared) private(iBin) \ reduction(+ \ : moment) \ schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++) { - log_next = log(bin_Edges[iBin+1]); - - switch(n) { - case 0: - pow_value = 1; - break; - case 1: - pow_value = bin_Centers[iBin]; - break; - case 2: - pow_value = bin_Centers[iBin] * bin_Centers[iBin]; - break; - case 3: - pow_value = bin_Centers[iBin] * bin_Centers[iBin] * bin_Centers[iBin]; - break; - default: - pow_value = pow(bin_Centers[iBin], n); // fallback for safety - break; - } + double pow_value; + if constexpr (N == 0) + pow_value = 1; + else if constexpr (N == 1) + pow_value = bin_Centers[iBin]; + else if constexpr (N == 2) + pow_value = bin_Centers[iBin] * bin_Centers[iBin]; + else if constexpr (N == 3) + pow_value = bin_Centers[iBin] * bin_Centers[iBin] * bin_Centers[iBin]; + else + pow_value = pow(bin_Centers[iBin], N); // fallback for safety - moment += (log_next - log_current) * pow_value * PDF[iBin]; - log_current = log_next; + moment += (log_Bin_Edges[iBin+1] - log_Bin_Edges[iBin]) * pow_value * PDF[iBin]; } return moment; - } /* End of Grid_Aerosol::Moment */ + } + double Grid_Aerosol::Moment(UInt n, const Vector_1D& PDF) const + { + switch(n) { + case 0: return Moment<0>(PDF); + case 1: return Moment<1>(PDF); + case 2: return Moment<2>(PDF); + case 3: return Moment<3>(PDF); + default: return Moment<4>(PDF); + } + }/* End of Grid_Aerosol::Moment */ - double Grid_Aerosol::Moment(UInt n, UInt jNy, UInt iNx) const + template + double Grid_Aerosol::Moment(UInt jNy, UInt iNx) const { UInt iBin = 0; @@ -1715,40 +1724,40 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { double moment = 0.0E+00; const double FACTOR = 3.0 / double(4.0 * PI); - double pow_value; - double log_current = log(bin_Edges[0]); - double log_next; - #pragma omp parallel for default(shared) private(iBin) \ reduction(+ \ : moment) \ schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++){ - log_next = log(bin_Edges[iBin+1]); - switch(n) { - case 0: - pow_value = 1; - break; - case 1: - pow_value = cbrt(FACTOR * bin_VCenters[iBin][jNy][iNx]); - break; - case 2: - pow_value = cbrt(FACTOR * FACTOR * bin_VCenters[iBin][jNy][iNx] * bin_VCenters[iBin][jNy][iNx]); - break; - case 3: - pow_value = FACTOR * bin_VCenters[iBin][jNy][iNx]; - break; - default: - pow_value = pow(FACTOR * bin_VCenters[iBin][jNy][iNx], n / 3.0); // fallback for safety - break; - } - moment += (log_next - log_current) * pow_value * pdf[iBin][jNy][iNx]; - log_current = log_next; + double pow_value; + if constexpr (N == 0) + pow_value = 1; + else if constexpr (N == 1) + pow_value = cbrt(FACTOR * bin_VCenters[iBin][jNy][iNx]); + else if constexpr (N == 2) + pow_value = cbrt(FACTOR * FACTOR * bin_VCenters[iBin][jNy][iNx] * bin_VCenters[iBin][jNy][iNx]); + else if constexpr (N == 3) + pow_value = FACTOR * bin_VCenters[iBin][jNy][iNx]; + else + pow_value = pow(FACTOR * bin_VCenters[iBin][jNy][iNx], N / 3.0); // fallback for safety + + moment += (log_Bin_Edges[iBin+1] - log_Bin_Edges[iBin]) * pow_value * pdf[iBin][jNy][iNx]; } return moment; - } /* End of Grid_Aerosol::Moment */ + } + double Grid_Aerosol::Moment(UInt n, UInt jNy, UInt iNx) const + { + switch(n) { + case 0: return Moment<0>(jNy, iNx); + case 1: return Moment<1>(jNy, iNx); + case 2: return Moment<2>(jNy, iNx); + case 3: return Moment<3>(jNy, iNx); + default: return Moment<4>(jNy, iNx); + } + } + /* End of Grid_Aerosol::Moment */ double Grid_Aerosol::Radius(UInt jNy, UInt iNx) const { From 0238044a39e049e863217ca0d7de72d9cd79b95b Mon Sep 17 00:00:00 2001 From: coco-yeung Date: Mon, 30 Mar 2026 16:14:52 +0800 Subject: [PATCH 5/5] Added comments --- Code.v05-00/include/AIM/Aerosol.hpp | 7 ++++--- Code.v05-00/src/AIM/Aerosol.cpp | 8 ++++++++ 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/Code.v05-00/include/AIM/Aerosol.hpp b/Code.v05-00/include/AIM/Aerosol.hpp index b03ef4fa..68808bba 100644 --- a/Code.v05-00/include/AIM/Aerosol.hpp +++ b/Code.v05-00/include/AIM/Aerosol.hpp @@ -61,6 +61,7 @@ class AIM::Aerosol /* Moments */ //NOTE: This gives the moment in [- / cm3]. You need to multiply by factors to get it in [ / m] or something. + //NOTE: now returning binMoment without the log bin ratio, as the only instance calling this function already calculates it inline double binMoment(int iBin, int n = 0) const { if (n == 0){ return pdf[iBin]; // if n = 0, it can skip the pow() calculation @@ -69,7 +70,7 @@ class AIM::Aerosol return pow(bin_Centers[iBin], n) * pdf[iBin]; } //NOTE: This gives the moment in [- / cm3]. You need to multiply by factors to get it in [ / m] or something. - template + template //Template allows us to evaulate the power at compile time to decrease runtime double Moment() const; double Moment( UInt n = 0 ) const; double Radius( ) const; @@ -132,10 +133,10 @@ class AIM::Grid_Aerosol inline void updateNy(int ny_new) { Ny = ny_new; }; /* Moments */ - template + template //Templates allow us to evaulate the power at compile time to decrease runtime Vector_2D Moment() const; Vector_2D Moment( UInt n ) const; - template + template double Moment( const Vector_1D& PDF ) const; double Moment( UInt n, const Vector_1D& PDF ) const; template diff --git a/Code.v05-00/src/AIM/Aerosol.cpp b/Code.v05-00/src/AIM/Aerosol.cpp index 8368985f..fc9acc99 100644 --- a/Code.v05-00/src/AIM/Aerosol.cpp +++ b/Code.v05-00/src/AIM/Aerosol.cpp @@ -317,6 +317,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { for (UInt iBin = 0; iBin < nBin; iBin++) { double pow_value; + // Moment function only used for N = 0, 1, 2, 3 --> eliminate usage of pow() to save time if constexpr (N == 0) pow_value = 1; else if constexpr (N == 1) @@ -334,6 +335,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { return moment; } + // overload function to call and evaluate branch in template function for different powers N at runtime double Aerosol::Moment(UInt n) const { switch(n) { @@ -1176,6 +1178,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { { for (iNx = 0; iNx < Nx; iNx++) { + // Moment function only used for N = 0, 1, 2, 3 --> eliminate usage of pow() to save time if constexpr (N == 0) pow_value = 1; else if constexpr (N == 1) @@ -1194,6 +1197,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { return moment; } + // overload function to call and evaluate branch in template function for different powers N at runtime Vector_2D Grid_Aerosol::Moment(UInt n) const { switch(n) { @@ -1686,6 +1690,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++) { + // Moment function only used for N = 0, 1, 2, 3 --> eliminate usage of pow() to save time double pow_value; if constexpr (N == 0) pow_value = 1; @@ -1704,6 +1709,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { return moment; } + // overload function to call and evaluate branch in template function for different powers N at runtime double Grid_Aerosol::Moment(UInt n, const Vector_1D& PDF) const { switch(n) { @@ -1730,6 +1736,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++){ double pow_value; + // Moment function only used for N = 0, 1, 2, 3 --> eliminate usage of pow() to save time if constexpr (N == 0) pow_value = 1; else if constexpr (N == 1) @@ -1747,6 +1754,7 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { return moment; } + // overload function to call and evaluate branch in template function for different powers N at runtime double Grid_Aerosol::Moment(UInt n, UInt jNy, UInt iNx) const { switch(n) {