diff --git a/Code.v05-00/CMakeLists.txt b/Code.v05-00/CMakeLists.txt index 895dad2a8..8264b543a 100644 --- a/Code.v05-00/CMakeLists.txt +++ b/Code.v05-00/CMakeLists.txt @@ -19,6 +19,9 @@ if (CMAKE_BUILD_TYPE STREQUAL "Debug") #add_link_options(-fsanitize=bounds -fsanitize=address -fsanitize=undefined) add_compile_options(-fno-omit-frame-pointer -fsanitize=bounds -fsanitize=undefined) add_link_options(-fsanitize=bounds -fsanitize=undefined) +elseif (CMAKE_BUILD_TYPE STREQUAL "Profiling") + include(CheckCXXCompilerFlag) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g -O2 -fno-omit-frame-pointer") elseif (NOT CMAKE_BUILD_TYPE OR CMAKE_BUILD_TYPE STREQUAL "") set(CMAKE_BUILD_TYPE "Release" CACHE STRING "" FORCE) include(CheckCXXCompilerFlag) @@ -132,4 +135,4 @@ add_subdirectory(tests) if (DEBUG) message(STATUS "DEBUG mode is enabled") -endif() +endif() \ No newline at end of file diff --git a/Code.v05-00/include/AIM/Aerosol.hpp b/Code.v05-00/include/AIM/Aerosol.hpp index 77b5c56f7..c259e076c 100644 --- a/Code.v05-00/include/AIM/Aerosol.hpp +++ b/Code.v05-00/include/AIM/Aerosol.hpp @@ -61,10 +61,17 @@ 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 (log(bin_Edges[iBin + 1]) - log(bin_Edges[iBin])) * pdf[iBin]; // if n = 0, it can skip the pow() calculation + } + return (log(bin_Edges[iBin + 1]) - log(bin_Edges[iBin])) * 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 allows us to evaulate the power at compile time to decrease runtime + double Moment() const; double Moment( UInt n = 0 ) const; double Radius( ) const; double EffRadius( ) const; @@ -79,6 +86,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; }; @@ -88,6 +96,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) }; @@ -124,8 +133,14 @@ class AIM::Grid_Aerosol inline void updateNy(int ny_new) { Ny = ny_new; }; /* Moments */ + 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 + 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 */ @@ -177,6 +192,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; }; @@ -189,6 +205,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/include/Util/MetFunction.hpp b/Code.v05-00/include/Util/MetFunction.hpp index 0c6cc4e67..cd20dc4e1 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 79523ab33..aa1ca2a90 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)) @@ -48,8 +49,12 @@ 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; + double bin_edge_cubed = bin_Edges[iBin] * bin_Edges[iBin] * bin_Edges[iBin]; + double next_bin_edge_cubed = bin_Edges[iBin + 1] * bin_Edges[iBin + 1] * bin_Edges[iBin + 1]; + bin_VCenters[iBin] = 4.0 / 3.0 * PI * (bin_edge_cubed + next_bin_edge_cubed) * 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) { @@ -84,8 +89,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 +102,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"; @@ -113,12 +120,17 @@ 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 * (pow(bin_Edges[iBin], 3) + pow(bin_Edges[iBin + 1], 3)) * 0.5; + double bin_edge_cubed = bin_Edges[iBin] * bin_Edges[iBin] * bin_Edges[iBin]; + double next_bin_edge_cubed = bin_Edges[iBin + 1] * bin_Edges[iBin + 1] * bin_Edges[iBin + 1]; + bin_VCenters[iBin] = 4.0 / 3.0 * PI * (bin_edge_cubed + next_bin_edge_cubed) * 0.5; + log_Bin_Edges[iBin] = log(bin_Edges[iBin]); } + log_Bin_Edges[nBin] = log(bin_Edges[nBin]); } @@ -300,19 +312,56 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { } /* End of Aerosol::Coagulate */ - double Aerosol::Moment(UInt n) const + template + double Aerosol::Moment() const { double moment = 0; for (UInt iBin = 0; iBin < nBin; iBin++) { - moment += (log(bin_Edges[iBin + 1]) - log(bin_Edges[iBin])) * pow(bin_Centers[iBin], n) * pdf[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) + 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 */ + } + // overload function to call and evaluate branch in template function for different powers N at runtime + 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: + { + // For arbitrary N, fall back to pow() directly + double moment = 0; + for (UInt iBin = 0; iBin < nBin; iBin++) + { + moment += (log_Bin_Edges[iBin+1] - log_Bin_Edges[iBin]) + * pow(bin_Centers[iBin], n) + * pdf[iBin]; + } + return moment; + } + } + } + /* End of Aerosol::Moment */ double Aerosol::Radius() const { @@ -355,7 +404,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 { @@ -401,6 +451,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)) @@ -417,10 +468,12 @@ 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]; + log_Bin_Edges[iBin] = log(bin_Edges[iBin]); } - bin_VEdges[nBin] = 4.0 / 3.0 * PI * pow(bin_Edges[nBin], 3); + 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))); @@ -485,10 +538,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 +563,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++) @@ -634,7 +689,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) { @@ -947,7 +1002,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 { @@ -1098,7 +1153,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++) @@ -1118,7 +1173,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; @@ -1132,17 +1188,65 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { schedule(dynamic, 1) if (!PARALLEL_CASES) for (iBin = 0; iBin < nBin; iBin++) { + double pow_value; 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]; + // 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) + 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]; } } } 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) { + case 0: return Moment<0>(); + case 1: return Moment<1>(); + case 2: return Moment<2>(); + case 3: return Moment<3>(); + default: + { + UInt jNy = 0; + UInt iNx = 0; + UInt iBin = 0; + + Vector_2D moment(Ny, Vector_1D(Nx, 0.0E+00)); + const double FACTOR = 3.0 / double(4.0 * PI); + + #pragma omp parallel for default(shared) private(iNx, jNy, iBin) \ + schedule(dynamic, 1) if (!PARALLEL_CASES) + for (iBin = 0; iBin < nBin; iBin++) + { + double pow_value; + for (jNy = 0; jNy < Ny; jNy++) + { + for (iNx = 0; iNx < Nx; iNx++) + { + pow_value = pow(FACTOR * bin_VCenters[iBin][jNy][iNx], n / 3.0); + moment[jNy][iNx] += (log_Bin_Edges[iBin+1] - log_Bin_Edges[iBin]) * pow_value * pdf[iBin][jNy][iNx]; + } + } + } + + return moment; + } + } } /* End of Grid_Aerosol::Moment */ Vector_3D Grid_Aerosol::Number() const @@ -1159,7 +1263,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++) @@ -1243,7 +1347,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++) @@ -1612,7 +1716,8 @@ 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; @@ -1625,14 +1730,54 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { 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]; + // 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; + 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 Grid_Aerosol::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) { + 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: + { + UInt iBin = 0; + double moment = 0.0E+00; - double Grid_Aerosol::Moment(UInt n, UInt jNy, UInt iNx) const + #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] - log_Bin_Edges[iBin]) * pow(bin_Centers[iBin], n) * PDF[iBin]; + } + + return moment; + } + } + }/* End of Grid_Aerosol::Moment */ + + template + double Grid_Aerosol::Moment(UInt jNy, UInt iNx) const { UInt iBin = 0; @@ -1644,12 +1789,53 @@ void Aerosol::addAerosolToPDF( const Aerosol &rhs ) { 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++){ + 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) + 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 */ + } + // 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) { + 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: + { + UInt iBin = 0; + double moment = 0.0E+00; + const double FACTOR = 3.0 / double(4.0 * PI); + + #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] - log_Bin_Edges[iBin]) * pow(FACTOR * bin_VCenters[iBin][jNy][iNx], n / 3.0) * pdf[iBin][jNy][iNx]; + } + + return moment; + } + } + } + /* End of Grid_Aerosol::Moment */ double Grid_Aerosol::Radius(UInt jNy, UInt iNx) const { @@ -1692,7 +1878,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 e4b206f54..2fd3ccd64 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 8dee9e3b3..ffb6cdcfd 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 242f1cce3..b72444418 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 895cb6a27..f66313656 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(); diff --git a/Code.v05-00/src/Core/LiquidAer.cpp b/Code.v05-00/src/Core/LiquidAer.cpp index 852526721..95523becc 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 static double SLA_VA = 8.406E-08 * pow( 10.0, 12.0 * 0.751E+00 ); /* Factor to convert effective radius to liquid radius (unitless) */ - const double SLA_RR = exp( -0.173E+00 ); + const static double SLA_RR = exp( -0.173E+00 ); /* 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 static double SLA_VR = 0.357E-06 * pow( 10.0, 12.0 * 0.249E+00 ); /* 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 d50000b7a..9fb29b77a 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 b95ec4b18..8dbd4bf7f 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 inv_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 ) * inv_denominator ); } UInt iTime_3mins; diff --git a/Code.v05-00/src/EPM/Solution.cpp b/Code.v05-00/src/EPM/Solution.cpp index 2304c7499..b0cde6ecc 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/LAGRID/RemappingFunctions.cpp b/Code.v05-00/src/LAGRID/RemappingFunctions.cpp index 67a59575f..f1b3564aa 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 458932e9c..c308d6a3e 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,8 +344,11 @@ 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; + double length_scale = 2.0 * r + l; + double length_scale_cubed = length_scale * length_scale * length_scale; + double pow_var = 4.0 * r * r + l * l; + return ( length_scale_cubed - pow_var * sqrt(pow_var) ) / ( 6.0 * r * l ) - 2.0 * r; } /* End of delta_p */ double Reynolds_p( const double r, const double rho, \ @@ -468,7 +471,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 +482,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; } diff --git a/examples/issl_rhi140/input.yaml b/examples/issl_rhi140/input.yaml index 686704ebb..b874d3edf 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