From e7514feb5a997be251abe9d0ae075798d63ca02f Mon Sep 17 00:00:00 2001 From: Aditya Dendukuri Date: Mon, 15 Jul 2024 20:47:04 -0600 Subject: [PATCH] assemble tridiagonal system --- etc/derecho/job.pbs | 24 ++ .../solvers/delta_eddington.hh | 221 ++++++++++++------ .../solvers/delta_eddington.hpp | 180 +++++++------- 3 files changed, 269 insertions(+), 156 deletions(-) create mode 100644 etc/derecho/job.pbs diff --git a/etc/derecho/job.pbs b/etc/derecho/job.pbs new file mode 100644 index 00000000..dfb1fac9 --- /dev/null +++ b/etc/derecho/job.pbs @@ -0,0 +1,24 @@ +#PBS -A NTDD0005 +#PBS -N hybrid_job +#PBS -q main +#PBS -l walltime=00:05:00 +#PBS -l select=1:ncpus=128:mpiprocs=1:ompthreads=1 + +# Load modules to match compile-time environment +module purge +module load ncarenv/23.09 intel-oneapi/2023.2.1 craype/2.7.23 cray-mpich/8.1.27 + +# move to the build directory +cd ../../ +rm -rf build +mkdir build +cd build + +# cmake command to generate build scrips +cmake -D CMAKE_BUILD_TYPE=release -D TUVX_ENABLE_MEMCHECK=OFF -D TUVX_ENABLE_REGRESSION_TESTS=OFF -D CMAKE_EXPORT_COMPILE_COMMANDS=1 -D TUVX_ENABLE_BENCHMARK=ON -D TUVX_ENABLE_TESTS=ON .. + +# build the benchmarking code target +make benchmark_tridiagonal_solver -j 8 + +# run the executable +./benchmark_tridiagonal_solver diff --git a/include/tuvx/radiative_transfer/solvers/delta_eddington.hh b/include/tuvx/radiative_transfer/solvers/delta_eddington.hh index 29803d7a..3c8767f2 100644 --- a/include/tuvx/radiative_transfer/solvers/delta_eddington.hh +++ b/include/tuvx/radiative_transfer/solvers/delta_eddington.hh @@ -1,8 +1,43 @@ // Copyright (C) 2023-2024 National Center for Atmospheric Research // SPDX-License-Identifier: Apache-2.0 #include "tuvx/linear_algebra/linear_algebra.hpp" +#include "tuvx/radiative_transfer/radiator.hpp" #include + +template +void DeltaEddingtonApproximation( + const RadiatorStatePolicy& accumulated_radiator_states, + const ArrayPolicy& parameters, + const std::vector solar_zenith_angles) +{ + const std::size_t number_of_columns = solar_zenith_angles.size(); + ArrayPolicy& optical_depth = accumulated_radiator_states.optical_depth_; + ArrayPolicy& single_scattering_albedo = accumulated_radiator_states.single_scattering_albedo_; + ArrayPolicy& assymetry_parameter = accumulated_radiator_states.assymetry_parameter_; + for (std::size_t i = 0; i < number_of_columns; i++) + { + // delta eddington parameters + auto& gamma1 = parameters[i][0]; + auto& gamma2 = parameters[i][1]; + auto& gamma3 = parameters[i][2]; + auto& gamma4 = parameters[i][3]; + auto& mu = parameters[i][4]; + + // radiator state variables + auto mu_0 = std::acos(solar_zenith_angles[i]); + auto& omega = optical_depth[i]; + auto& g = single_scattering_albedo[i]; + auto& tau = optical_depth[i]; + + // compute delta eddington parameters + gamma1 = 7 - omega * (4 + 3 * g); + gamma2 = -(1 - omega * (4 - 3 * g)) / 4; + gamma3 = (2 - 3 * g * mu_0) / 4; + mu = (T)0.5; + } +} + namespace tuvx { @@ -24,6 +59,7 @@ namespace tuvx const std::vector& solar_zenith_angles, const std::map& grids, const std::map& profiles, + const RadiatorStatePolicy& accumulated_radiator_states) { // determine number of layers @@ -56,46 +92,136 @@ namespace tuvx // TODO slant optical depth computation tau = tau; } + } + + template + void AssembleTridiagonalSystem( + const std::vector& solar_zenith_angles, + const std::map& grids, + const std::map& profiles, + const std::vector solution_parameters, + const TridiagonalMatrix& coeffcient_matrix, + const std::vector& coeffcient_vector) + { + // get system size + std::size_t system_size = coeffcient_matrix.main_diagonal_.size() / 2; + + // coeffcient matrix diagonals + std::vector& upper_diagonal = coeffcient_matrix.upper_diagonal_; + std::vector& main_diagonal = coeffcient_matrix.main_diagonal_; + std::vector& lower_diagonal = coeffcient_matrix.lower_diagonal; - // delta eddington parameters (little gammas 1, 2, 3 and mu from the paper) - ArrayPolicy delta_eddington_parameters(number_of_columns, 5); { - for (std::size_t i = 0; i < number_of_columns; i++) + // first row + lower_diagonal[0] = 0; + main_diagonal[0] = solution_parameters[0][0]; + upper_diagonal[0] = -solution_parameters[1][0]; + + // even rows + for (std::size_t n = 2; n < 2 * system_size / 2; n += 2) { - // delta eddington parameters - auto& gamma1 = delta_eddington_parameters[i][0]; - auto& gamma2 = delta_eddington_parameters[i][1]; - auto& gamma3 = delta_eddington_parameters[i][2]; - auto& gamma4 = delta_eddington_parameters[i][3]; - auto& mu = delta_eddington_parameters[i][4]; - - // radiator state variables - auto mu_0 = std::acos(solar_zenith_angles[i]); - auto& omega = optical_depth[i]; - auto& g = single_scattering_albedo[i]; - auto& tau = optical_depth[i]; - - // compute delta eddington parameters - gamma1 = 7 - omega * (4 + 3 * g); - gamma2 = -(1 - omega * (4 - 3 * g)) / 4; - gamma3 = (2 - 3 * g * mu_0) / 4; - mu = (T)0.5; + const auto& e1 = solution_parameters[n][0]; + const auto& e2 = solution_parameters[n][1]; + const auto& e3 = solution_parameters[n][2]; + const auto& e4 = solution_parameters[n][3]; + upper_diagonal[n] = e2[n + 1] * e1[n] - e3[n] * e4[n + 1]; + main_diagonal[n] = e2[n] * e2[n + 1] - e4[n] * e4[n + 1]; + lower_diagonal[n] = e1[n] * e4[n + 1] - e2[n + 1] * e3[n + 1]; } + + // odd rows + for (std::size_t n = 1; n < 2 * system_size / 2 - 1; n += 2) + { + const auto& e1 = solution_parameters[n][0]; + const auto& e2 = solution_parameters[n][1]; + const auto& e3 = solution_parameters[n][2]; + const auto& e4 = solution_parameters[n][3]; + upper_diagonal[n] = e2[n] * e3[n] - e4[n] * e1[n]; + main_diagonal[n] = e1[n] * e1[n + 1] - e3[n] * e3[n + 1]; + lower_diagonal[n] = e3[n] * e4[n + 1] - e1[n + 1] * e2[n + 1]; + } + + // last row + const auto& e1 = solution_parameters[n][0]; + const auto& e2 = solution_parameters[n][1]; + const auto& e3 = solution_parameters[n][2]; + const auto& e4 = solution_parameters[n][3]; + std::size_t N = 2 * system_size - 1; + upper_diagonal[N] = e1[N]; } } + /// @brief Initilize the variables required for the delta eddington approximation + /// @param solar zenith angles Solar zenith angles for each column [radians] + /// @param grids Grids available for the radiative transfer calculation + /// @param profiles Profiles available for the radiative transfer calculation + /// @param radiation_field The calculated radiation field + /// + /// Delta scaling of the inputs, compuation of reflectivity, diffuse flux and + /// delta eddington coeffcients. template< typename T, typename GridPolicy, typename ProfilePolicy, + typename RadiatorStatePolicytypename, + typename RadiationFieldPolicy> + void ComputeRadiationField( + const std::vector& solar_zenith_angles, + const std::map& grids, + const std::map& profiles, + const Array2D solution_parameters, + const TridiagonalMatrix& coeffcient_matrix, + const std::vector& coeffcient_vector, + const RadiationFieldPolicy& radiation_field) + { + // [DEV NOTES] Temporarily return predictable values for the radiation field. + // This will be replaced with the actual results once the solver is implemented. + int offset = 42; + for (auto& elem : radiation_field.spectral_irradiance_.direct_) + { + elem = offset++; + } + offset = 93; + for (auto& elem : radiation_field.spectral_irradiance_.upwelling_) + { + elem = offset++; + } + offset = 52; + for (auto& elem : radiation_field.spectral_irradiance_.downwelling_) + { + elem = offset++; + } + offset = 5; + for (auto& elem : radiation_field.actinic_flux_.direct_) + { + elem = offset++; + } + offset = 24; + for (auto& elem : radiation_field.actinic_flux_.upwelling_) + { + elem = offset++; + } + offset = 97; + for (auto& elem : radiation_field.actinic_flux_.downwelling_) + { + elem = offset++; + } + } + + template< + typename T, + typename ArrayPolicy, + typename GridPolicy, + typename ProfilePolicy, typename RadiatorStatePolicy, typename RadiationFieldPolicy> - inline void DeltaEddington::Solve( + inline void Solve( const std::vector& solar_zenith_angles, const std::map& grids, const std::map& profiles, + const std::function)> ApproximationFunction, const RadiatorStatePolicy& accumulated_radiator_state, - RadiationFieldPolicy& radiation_field) const + RadiationFieldPolicy& radiation_field) { // Solve the radiative transfer equation. // @@ -119,59 +245,24 @@ namespace tuvx // internal solver variables Array2D solution_parameters; - Array2D delta_eddington_parameters; + Array2D simulation_parameters; // tridiagonal system variables TridiagonalMatrix coeffcient_matrix; std::vector coeffcient_vector; - // Initialize variables - this->InitializeVariables( + tuvx::InitializeVariables( solar_zenith_angles, grids, profiles, accumulated_radiator_state); - // Assemble the tridiagonal system - this->AssembleTridiagonalSystem( + ApproximationFunction(accumulated_radiator_state, solar_zenith_angles, simulation_parameters); + + tuvx::AssembleTridiagonalSystem( solar_zenith_angles, grids, profiles, solution_parameters, coeffcient_matrix, coeffcient_vector); - // solve tridiagonal system - Solve(coeffcient_matrix, coeffcient_vector); + tuvx::Solve(coeffcient_matrix, coeffcient_vector); - // Update radiation field - this->ComputeRadiationField( + tuvx::ComputeRadiationField( solar_zenith_angles, grids, profiles, solution_parameters, coeffcient_matrix, coeffcient_vector, radiation_field); - - // [DEV NOTES] Temporarily return predictable values for the radiation field. - // This will be replaced with the actual results once the solver is implemented. - int offset = 42; - for (auto& elem : radiation_field.spectral_irradiance_.direct_) - { - elem = offset++; - } - offset = 93; - for (auto& elem : radiation_field.spectral_irradiance_.upwelling_) - { - elem = offset++; - } - offset = 52; - for (auto& elem : radiation_field.spectral_irradiance_.downwelling_) - { - elem = offset++; - } - offset = 5; - for (auto& elem : radiation_field.actinic_flux_.direct_) - { - elem = offset++; - } - offset = 24; - for (auto& elem : radiation_field.actinic_flux_.upwelling_) - { - elem = offset++; - } - offset = 97; - for (auto& elem : radiation_field.actinic_flux_.downwelling_) - { - elem = offset++; - } } } // namespace tuvx diff --git a/include/tuvx/radiative_transfer/solvers/delta_eddington.hpp b/include/tuvx/radiative_transfer/solvers/delta_eddington.hpp index dcbd2526..35950e1f 100644 --- a/include/tuvx/radiative_transfer/solvers/delta_eddington.hpp +++ b/include/tuvx/radiative_transfer/solvers/delta_eddington.hpp @@ -22,99 +22,97 @@ namespace tuvx /// @brief Radiative flux calculator that applies the delta-Eddington Approximation. /// /// [DEV NOTES] We can determine whether this should be a class or a set of functions - template> - class DeltaEddington - { - public: - /// Construct a Delta-Eddington solver. - DeltaEddington() = default; - /// @brief Solve the radiative transfer equation for a collection of columns - /// @param solar_zenith_angles Solar zenith angles for each column [radians]. - /// @param grids Grids available for the radiative transfer calculation. - /// @param profiles Profiles available for the radiative transfer calculation. - /// @param radiation_field The calculated radiation field. - /// - /// Solves two-stream equations for multiple layers. These routines are based - /// on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989. - /// DOI: https://doi.org/10.1029/JD094iD13p16287 - /// It contains 9 two-stream methods to choose from. A pseudo-spherical - /// correction has also been added. - /// - /// The original delta-Eddington paper is: - /// Joseph and Wiscombe, J. Atmos. Sci., 33, 2453-2459, 1976 - /// DOI: https://doi.org/10.1175/1520-0469(1976)033%3C2452:TDEAFR%3E2.0.CO;2 - template< - typename T, - typename GridPolicy, - typename ProfilePolicy, - typename RadiatorStatePolicy, - typename RadiationFieldPolicy> - void Solve( - const std::vector& solar_zenith_angles, - const std::map& grids, - const std::map& profiles, - const RadiatorStatePolicy& accumulated_radiator_states, - RadiationFieldPolicy& radiation_field) const; + /// @brief Solve the radiative transfer equation for a collection of columns + /// @param solar_zenith_angles Solar zenith angles for each column [radians]. + /// @param grids Grids available for the radiative transfer calculation. + /// @param profiles Profiles available for the radiative transfer calculation. + /// @param radiation_field The calculated radiation field. + /// + /// Solves two-stream equations for multiple layers. These routines are based + /// on equations from: Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989. + /// DOI: https://doi.org/10.1029/JD094iD13p16287 + /// It contains 9 two-stream methods to choose from. A pseudo-spherical + /// correction has also been added. + /// + /// The original delta-Eddington paper is: + /// Joseph and Wiscombe, J. Atmos. Sci., 33, 2453-2459, 1976 + /// DOI: https://doi.org/10.1175/1520-0469(1976)033%3C2452:TDEAFR%3E2.0.CO;2 + template< + typename T, + typename ArrayPolicy, + typename GridPolicy, + typename ProfilePolicy, + typename RadiatorStatePolicy, + typename RadiationFieldPolicy> + void Solve( + const std::vector& solar_zenith_angles, + const std::map& grids, + const std::map& profiles, + const RadiatorStatePolicy& accumulated_radiator_states, + const std::function approximation_function, + RadiationFieldPolicy& radiation_field); - /// @brief Initilize the variables required for the delta eddington approximation - /// @param solar zenith angles Solar zenith angles for each column [radians] - /// @param grids Grids available for the radiative transfer calculation - /// @param profiles Profiles available for the radiative transfer calculation - /// @param radiation_field The calculated radiation field - /// - /// Delta scaling of the inputs, compuation of reflectivity, diffuse flux - template< - typename T, - typename GridPolicy, - typename ProfilePolicy, - typename RadiatorStatePolicy, - typename RadiationFieldPolicy> - void InitializeVariables( - const std::vector& solar_zenith_angles, - const std::map& grids, - const std::map& profiles, - const RadiatorStatePolicy& accumulated_radiator_states); + /// @brief Initilize the variables required for the delta eddington approximation + /// @param solar zenith angles Solar zenith angles for each column [radians] + /// @param grids Grids available for the radiative transfer calculation + /// @param profiles Profiles available for the radiative transfer calculation + /// @param radiation_field The calculated radiation field + /// + /// Delta scaling of the inputs, compuation of reflectivity, diffuse flux + template< + typename T, + typename GridPolicy, + typename ProfilePolicy, + typename RadiatorStatePolicy, + typename RadiationFieldPolicy> + void InitializeVariables( + const std::vector& solar_zenith_angles, + const std::map& grids, + const std::map& profiles, + const RadiatorStatePolicy& accumulated_radiator_states); - /// @brief Initilize the variables required for the delta eddington approximation - /// @param solar zenith angles Solar zenith angles for each column [radians] - /// @param grids Grids available for the radiative transfer calculation - /// @param profiles Profiles available for the radiative transfer calculation - /// @param radiation_field The calculated radiation field - /// - /// Delta scaling of the inputs, compuation of reflectivity, diffuse flux and - /// delta eddington coeffcients. - template - void AssembleTridiagonalSystem( - const std::vector& solar_zenith_angles, - const std::map& grids, - const std::map& profiles, - const Array2D solution_parameters, - const TridiagonalMatrix& coeffcient_matrix, - const std::vector& coeffcient_vector); + /// @brief Initilize the variables required for the delta eddington approximation + /// @param solar zenith angles Solar zenith angles for each column [radians] + /// @param grids Grids available for the radiative transfer calculation + /// @param profiles Profiles available for the radiative transfer calculation + /// @param radiation_field The calculated radiation field + /// + /// Delta scaling of the inputs, compuation of reflectivity, diffuse flux and + /// delta eddington coeffcients. + template + void AssembleTridiagonalSystem( + const std::vector& solar_zenith_angles, + const std::map& grids, + const std::map& profiles, + const Array2D solution_parameters, + const TridiagonalMatrix& coeffcient_matrix, + const std::vector& coeffcient_vector); - /// @brief Initilize the variables required for the delta eddington approximation - /// @param solar zenith angles Solar zenith angles for each column [radians] - /// @param grids Grids available for the radiative transfer calculation - /// @param profiles Profiles available for the radiative transfer calculation - /// @param radiation_field The calculated radiation field - /// - /// Delta scaling of the inputs, compuation of reflectivity, diffuse flux and - /// delta eddington coeffcients. - template< - typename T, - typename GridPolicy, - typename ProfilePolicy, - typename RadiatorStatePolicytypename, - typename RadiationFieldPolicy> - void ComputeRadiationField( - const std::vector& solar_zenith_angles, - const std::map& grids, - const std::map& profiles, - const Array2D solution_parameters, - const TridiagonalMatrix& coeffcient_matrix, - const std::vector& coeffcient_vector, - const RadiationFieldPolicy& radiation_field); - }; + /// @brief Initilize the variables required for the delta eddington approximation + /// @param solar zenith angles Solar zenith angles for each column [radians] + /// @param grids Grids available for the radiative transfer calculation + /// @param profiles Profiles available for the radiative transfer calculation + /// @param radiation_field The calculated radiation field + /// + /// Delta scaling of the inputs, compuation of reflectivity, diffuse flux and + /// delta eddington coeffcients. + template< + typename T, + typename GridPolicy, + typename ProfilePolicy, + typename RadiatorStatePolicytypename, + typename RadiationFieldPolicy> + void ComputeRadiationField( + const std::vector& solar_zenith_angles, + const std::map& grids, + const std::map& profiles, + const Array2D solution_parameters, + const TridiagonalMatrix& coeffcient_matrix, + const std::vector& coeffcient_vector, + const RadiationFieldPolicy& radiation_field); -} // namespace tuvx + /// @brief Compute the Eddington parameters + template + void DeltaEddingtonApproximation(const RadiatorStatePolicy& radiator_state, const ArrayPolicy& parameters); +}; // namespace tuvx