From 37b9c073671515a5bcf87cef970ba0800e94fd03 Mon Sep 17 00:00:00 2001 From: Aditya Dendukuri Date: Wed, 17 Jul 2024 11:47:59 -0600 Subject: [PATCH] split solver code into multiple files --- "\\" | 332 ------------------ .../solvers/assemble_tridiagonal_system.hh | 112 ++++++ .../solvers/compute_radiation_field.hh | 62 ++++ .../solvers/delta_eddington.hh | 292 +-------------- .../solvers/delta_eddington.hpp | 50 +-- .../solvers/initialize_solver.hh | 129 +++++++ .../tuvx/radiative_transfer/solvers/solve.hh | 74 ++++ 7 files changed, 398 insertions(+), 653 deletions(-) delete mode 100644 "\\" create mode 100644 include/tuvx/radiative_transfer/solvers/assemble_tridiagonal_system.hh create mode 100644 include/tuvx/radiative_transfer/solvers/compute_radiation_field.hh create mode 100644 include/tuvx/radiative_transfer/solvers/initialize_solver.hh create mode 100644 include/tuvx/radiative_transfer/solvers/solve.hh diff --git "a/\\" "b/\\" deleted file mode 100644 index fef9fe4f..00000000 --- "a/\\" +++ /dev/null @@ -1,332 +0,0 @@ -// 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 - -#include - -namespace tuvx -{ - template - void DeltaEddingtonApproximation( - const RadiatorStatePolicy& accumulated_radiator_states, - const std::map> solution_parameters, - const std::vector solar_zenith_angles) - { - const std::size_t number_of_columns = solar_zenith_angles.size(); - ArrayPolicy& omega = accumulated_radiator_states.optical_depth_; - ArrayPolicy& g = accumulated_radiator_states.single_scattering_albedo_; - ArrayPolicy& tau = accumulated_radiator_states.assymetry_parameter_; - - // delta eddington parameters - std::vector& gamma1 = solution_parameters.at("gamma1"); - std::vector& gamma2 = solution_parameters.at("gamma2"); - std::vector& gamma3 = solution_parameters.at("gamma3"); - std::vector& gamma4 = solution_parameters.at("gamma4"); - std::vector& mu = solution_parameters.at("mu"); - std::vector& lambda = solution_parameters.at("mu"); - std::vector& Gamma = solution_parameters.at("mu"); - - // simulation parameters - T mu_0; - for (std::size_t i = 0; i < number_of_columns; i++) - { - // compute delta eddington parameters - mu_0 = std::acos(solar_zenith_angles[i]); - gamma1[i] = 7 - omega[i] * (4 + 3 * g[i]); - gamma2[i] = -(1 - omega[i] * (4 - 3 * g[i])) / 4; - gamma3[i] = (2 - 3 * g[i] * mu_0) / 4; - lambda[i] = std::sqrt(gamma1[i] * gamma1[i] - gamma2[i] * gamma2[i]); - Gamma[i] = std::sqrt(gamma1[i] * gamma1[i] - gamma2[i] * gamma2[i]); - mu = (T)0.5; - } - } - - template< - typename T, - typename GridPolicy, - typename ProfilePolicy, - typename RadiatorStatePolicy, - typename RadiationFieldPolicy, - typename ArrayPolicy> - inline void InitializeVariables( - const std::vector& solar_zenith_angles, - const std::map& grids, - const std::map& profiles, - const std::map& solver_parameters, - const std::map& solution_parameters, - const std::map& source_terms, - const RadiatorStatePolicy& accumulated_radiator_states) - { - // determine number of layers - const std::size_t number_of_columns = solar_zenith_angles.size(); - const auto& vertical_grid = grids.at("altitude [m]"); - const auto& wavelength_grid = grids.at("wavelength [m]"); - - // Check for consistency between the grids and profiles. - assert(vertical_grid.NumberOfColumns() == number_of_columns); - assert(wavelength_grid.NumberOfColumns() == 1); - - // Radiator state variables - ArrayPolicy& tau = accumulated_radiator_states.optical_depth_; - ArrayPolicy& omega = accumulated_radiator_states.single_scattering_albedo_; - ArrayPolicy& g = accumulated_radiator_states.assymetry_parameter_; - - // Delta scaling - T f; - for (std::size_t i = 0; i < number_of_columns; i++) - { - f = omega[i] * omega[i]; - omega[i] = (omega - f) / (1 - f); - g[i] = (1 - f) * g[i] / (1 - g[i] * f); - omega[i] = (1 - g[i] * f) * omega[i]; - } - - // TODO slant optical depth computation - for (auto& tau_n : tau) - { - tau_n = tau_n; - } - - { - // Source terms (C1 and C2 from the paper) - auto& C_upwelling = source_terms.at("C_upwelling"); - auto& C_downwelling = source_terms.at("C_downwelling"); - auto& S_sfc_i = solution_parameters.at("infrared source flux"); - auto& S_sfc_s = solution_parameters.at("solar source flux"); - - // other parameters - auto& lambda = solution_parameters.at("lambda"); - std::vector& gamma1 = solution_parameters.at("gamma1"); - std::vector& gamma2 = solution_parameters.at("gamma2"); - std::vector& gamma3 = solution_parameters.at("gamma3"); - std::vector& gamma4 = solution_parameters.at("gamma4"); - std::vector& mu = solution_parameters.at("mu"); - - auto& R_sfc = solution_parameters.at("source flux"); - - // temporary variables - T tau_cumulative = 0; - T exponential_term, denominator_term, mu_0; - - // source terms (C equations from 16, 290; eqns 23, 24) - for (std::size_t i = 0; i < number_of_columns; i++) - { - mu_0 = std::acos(solar_zenith_angles[i]); - exponential_term = omega * M_PI * R_sfc * std::exp(-(tau_cumulative - tau[i]) / mu_0); - denominator_term = (lambda * lambda - 1 / (mu_0 * mu_0)); - tau_cumulative += tau[i]; - - S_sfc_i[i] = R_sfc * mu_0 * std::exp(-tau_cumulative / mu_0); - S_sfc_s[i] = M_PI * R_sfc; - C_downwelling[i] = exponential_term * (((gamma1[i] + 1) / mu_0) * gamma4[i] + gamma2[i] * gamma3[i]); - C_upwelling[i] = exponential_term * (((gamma1[i] - 1) / mu_0) * gamma3[i] + gamma4[i] * gamma2[i]); - } - } - } - - template - void AssembleTridiagonalMatrix( - std::size_t number_of_layers, - const std::map> solution_parameters, - const std::map solver_parameters, - const TridiagonalMatrix& coeffcient_matrix) - { - // get linear system size - std::size_t matrix_size = 2 * number_of_layers; - { - // LEFT HAND SIDE 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_; - - // extract internal variables to build the matrix - const std::vector& e1 = solution_parameters.at("e1"); - const std::vector& e2 = solution_parameters.at("e2"); - const std::vector& e3 = solution_parameters.at("e3"); - const std::vector& e4 = solution_parameters.at("e4"); - - // extract surface reflectivity - const T& R_sfc = solver_parameters.at("Surface Reflectivity"); - - // first row - upper_diagonal.front() = 0; - main_diagonal.front() = e1.front(); - lower_diagonal.front() = -e2.front(); - - // odd rows - for (std::size_t n = 1; n < matrix_size - 1; n += 2) - { - upper_diagonal[n] = e2[n + 1] * e1[n] - e3[n] * e4[n + 1]; - main_diagonal[n] = e2[n] * e2[n + 1] - e3[n] * e4[n + 1]; - lower_diagonal[n] = e3[n] * e4[n + 1] - e1[n + 1] * e2[n + 1]; - } - - // even rows - for (std::size_t n = 2; n < matrix_size - 2; n += 2) - { - 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 - lower_diagonal.back() = e1.back() - R_sfc * e3.back(); - main_diagonal.back() = e2.back() - R_sfc * e4.back(); - upper_diagonal.back() = 0; - } - } - - template - void AssembleCoeffcientVector( - std::size_t number_of_layers, - const std::map> solution_parameters, - const std::map> source_terms, - const std::map solver_parameters, - std::vector& coeffcient_vector) - { - // get linear system size - std::size_t matrix_size = 2 * number_of_layers; - { - // extract internal variables to build the matrix - const std::vector& e1 = solution_parameters.at("e1"); - const std::vector& e2 = solution_parameters.at("e2"); - const std::vector& e3 = solution_parameters.at("e3"); - const std::vector& e4 = solution_parameters.at("e4"); - - // extract surface reflectivity and flux source - const T& R_sfc = solver_parameters.at("Surface Reflectivity"); - const T& f_0 = solver_parameters.at("source flux"); - - // extract source terms - const auto& C_upwelling = source_terms.at("C_upwelling"); - const auto& C_downwelling = source_terms.at("C_downwelling"); - - // first row - coeffcient_vector.front() = f_0 - C_downwelling.front(); - - // odd rows - for (std::size_t n = 1; n < matrix_size - 1; n += 2) - { - coeffcient_vector[n] = - e3[n] * (C_upwelling.start() - C_upwelling[n]) + e1[n] * (C_downwelling[n] - C_downwelling[0]); - } - - // even rows - for (std::size_t n = 2; n < matrix_size - 2; n += 2) - { - coeffcient_vector[n] = - e2[n + 1] * (C_upwelling.start() - C_upwelling[n]) + e4[n + 1] * (C_downwelling.start() - C_downwelling[n]); - } - - // last row - } - } - - 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 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 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) - { - // Solve the radiative transfer equation. - // - // [DEV NOTES] This is a placeholder for the actual implementation. - // The spherical geometry argument of the original solver was left out - // until we determine whether it needs to be an object or just a set of functions. - // - // Things that will change from the original solver: - // 1. All variables will be in SI units. Some of the original solver's - // variables were in non-SI units. - // 2. We will be solving for collections of columns. The original solver - // was for a single column. - // 3. The variable naming and source-code documentation will be improved. - const std::size_t number_of_columns = solar_zenith_angles.size(); - const auto& vertical_grid = grids.at("altitude [m]"); - const auto& wavelength_grid = grids.at("wavelength [m]"); - - // Check for consistency between the grids and profiles. - assert(vertical_grid.NumberOfColumns() == number_of_columns); - assert(wavelength_grid.NumberOfColumns() == 1); - - // internal solver variables - Array2D solution_parameters; - Array2D simulation_parameters; - - // tridiagonal system variables - TridiagonalMatrix coeffcient_matrix; - std::vector coeffcient_vector; - - tuvx::InitializeVariables( - solar_zenith_angles, grids, profiles, accumulated_radiator_state); - - ApproximationFunction(accumulated_radiator_state, solar_zenith_angles, simulation_parameters); - - tuvx::AssembleTridiagonalSystem( - solar_zenith_angles, grids, profiles, solution_parameters, coeffcient_matrix, coeffcient_vector); - - tuvx::Solve(coeffcient_matrix, coeffcient_vector); - - tuvx::ComputeRadiationField( - solar_zenith_angles, grids, profiles, solution_parameters, coeffcient_matrix, coeffcient_vector, radiation_field); - } - -} // namespace tuvx diff --git a/include/tuvx/radiative_transfer/solvers/assemble_tridiagonal_system.hh b/include/tuvx/radiative_transfer/solvers/assemble_tridiagonal_system.hh new file mode 100644 index 00000000..d11e5d0c --- /dev/null +++ b/include/tuvx/radiative_transfer/solvers/assemble_tridiagonal_system.hh @@ -0,0 +1,112 @@ + +#include "tuvx/linear_algebra/linear_algebra.hpp" + +#include + +#include + +namespace tuvx +{ + + template + void AssembleTridiagonalMatrix( + std::size_t number_of_layers, + const std::map> solution_parameters, + const std::map solver_parameters, + const TridiagonalMatrix& coeffcient_matrix) + { + // get linear system size + std::size_t matrix_size = 2 * number_of_layers; + { + // LEFT HAND SIDE 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_; + + // extract internal variables to build the matrix + const std::vector& e1 = solution_parameters.at("e1"); + const std::vector& e2 = solution_parameters.at("e2"); + const std::vector& e3 = solution_parameters.at("e3"); + const std::vector& e4 = solution_parameters.at("e4"); + + // extract surface reflectivity + const T& R_sfc = solver_parameters.at("Surface Reflectivity"); + + // first row + upper_diagonal.front() = 0; + main_diagonal.front() = e1.front(); + lower_diagonal.front() = -e2.front(); + + // odd rows + for (std::size_t n = 1; n < matrix_size - 1; n += 2) + { + upper_diagonal[n] = e2[n + 1] * e1[n] - e3[n] * e4[n + 1]; + main_diagonal[n] = e2[n] * e2[n + 1] - e3[n] * e4[n + 1]; + lower_diagonal[n] = e3[n] * e4[n + 1] - e1[n + 1] * e2[n + 1]; + } + + // even rows + for (std::size_t n = 2; n < matrix_size - 2; n += 2) + { + 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 + lower_diagonal.back() = e1.back() - R_sfc * e3.back(); + main_diagonal.back() = e2.back() - R_sfc * e4.back(); + upper_diagonal.back() = 0; + } + } + + template + void AssembleCoeffcientVector( + std::size_t number_of_layers, + const std::map> solution_parameters, + const std::map> source_terms, + const std::map solver_parameters, + std::vector& coeffcient_vector) + { + // set up tridiagonal matrix + std::size_t matrix_size = 2 * number_of_layers; + { + // extract internal variables to build the matrix + const std::vector& e1 = solution_parameters.at("e1"); + const std::vector& e2 = solution_parameters.at("e2"); + const std::vector& e3 = solution_parameters.at("e3"); + const std::vector& e4 = solution_parameters.at("e4"); + + // extract surface reflectivity and flux source + const std::vector& R_sfc = solver_parameters.at("Surface Reflectivity"); + const std::vector& f_0 = solver_parameters.at("Initial Source Flux"); + const std::vector& tau = solver_parameters.at("Optical Depth"); + const std::vector& S_sfc = solver_parameters.at("Solar Reflectivity"); + + // extract source functions + // this is a list of std::functions + const auto& C_upwelling = source_terms.at("C_upwelling"); + const auto& C_downwelling = source_terms.at("C_downwelling"); + + // first row + coeffcient_vector.front() = f_0 - C_downwelling[0](tau[0]); + + // odd rows + for (std::size_t n = 1; n < matrix_size - 1; n += 2) + { + coeffcient_vector[n] = e3[n] * (C_upwelling[n + 1](0) - C_upwelling[n](tau[n])) + + e1[n] * (C_downwelling[n](tau[n]) - C_downwelling[n](0)); + } + + // even rows + for (std::size_t n = 2; n < matrix_size - 2; n += 2) + { + coeffcient_vector[n] = e2[n + 1] * (C_upwelling[n + 1](0) - C_upwelling[n](tau[n])) + + e1[4 * n + 1] * (C_downwelling[n + 1](0) - C_downwelling[n + 1](tau[n])); + } + + // last row + coeffcient_vector.back() = S_sfc - C_upwelling.back()(tau.back()) + R_sfc * C_downwelling.back()(tau.back()); + } + } +} // namespace tuvx diff --git a/include/tuvx/radiative_transfer/solvers/compute_radiation_field.hh b/include/tuvx/radiative_transfer/solvers/compute_radiation_field.hh new file mode 100644 index 00000000..85928fb7 --- /dev/null +++ b/include/tuvx/radiative_transfer/solvers/compute_radiation_field.hh @@ -0,0 +1,62 @@ + + +// 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 + +#include + +namespace tuvx +{ + + 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 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++; + } + } + +} // namespace tuvx diff --git a/include/tuvx/radiative_transfer/solvers/delta_eddington.hh b/include/tuvx/radiative_transfer/solvers/delta_eddington.hh index 7898d1b7..c5dc8b91 100644 --- a/include/tuvx/radiative_transfer/solvers/delta_eddington.hh +++ b/include/tuvx/radiative_transfer/solvers/delta_eddington.hh @@ -26,8 +26,8 @@ namespace tuvx std::vector& gamma3 = solution_parameters.at("gamma3"); std::vector& gamma4 = solution_parameters.at("gamma4"); std::vector& mu = solution_parameters.at("mu"); - std::vector& lambda = solution_parameters.at("mu"); - std::vector& Gamma = solution_parameters.at("mu"); + std::vector& lambda = solution_parameters.at("lambda"); + std::vector& gamma = solution_parameters.at("Gamma"); // simulation parameters T mu_0; @@ -39,294 +39,8 @@ namespace tuvx gamma2[i] = -(1 - omega[i] * (4 - 3 * g[i])) / 4; gamma3[i] = (2 - 3 * g[i] * mu_0) / 4; lambda[i] = std::sqrt(gamma1[i] * gamma1[i] - gamma2[i] * gamma2[i]); - Gamma[i] = std::sqrt(gamma1[i] * gamma1[i] - gamma2[i] * gamma2[i]); + gamma[i] = (gamma1[i] - lambda[i]) / gamma2[i]; mu = (T)0.5; } } - - template< - typename T, - typename GridPolicy, - typename ProfilePolicy, - typename RadiatorStatePolicy, - typename RadiationFieldPolicy, - typename ArrayPolicy> - inline void InitializeVariables( - const std::vector& solar_zenith_angles, - const std::map& grids, - const std::map& profiles, - const std::map& solver_parameters, - const std::map& solution_parameters, - const std::map& source_terms, - const RadiatorStatePolicy& accumulated_radiator_states) - { - // determine number of layers - const std::size_t number_of_columns = solar_zenith_angles.size(); - const auto& vertical_grid = grids.at("altitude [m]"); - const auto& wavelength_grid = grids.at("wavelength [m]"); - - // Check for consistency between the grids and profiles. - assert(vertical_grid.NumberOfColumns() == number_of_columns); - assert(wavelength_grid.NumberOfColumns() == 1); - - // Radiator state variables - ArrayPolicy& tau = accumulated_radiator_states.optical_depth_; - ArrayPolicy& omega = accumulated_radiator_states.single_scattering_albedo_; - ArrayPolicy& g = accumulated_radiator_states.assymetry_parameter_; - - // Delta scaling - T f; - for (std::size_t i = 0; i < number_of_columns; i++) - { - f = omega[i] * omega[i]; - omega[i] = (omega - f) / (1 - f); - g[i] = (1 - f) * g[i] / (1 - g[i] * f); - omega[i] = (1 - g[i] * f) * omega[i]; - } - - // TODO slant optical depth computation - for (auto& tau_n : tau) - { - tau_n = tau_n; - } - - { - // Source terms (C1 and C2 from the paper) - auto& C_upwelling = source_terms.at("C_upwelling"); - auto& C_downwelling = source_terms.at("C_downwelling"); - auto& S_sfc_i = solution_parameters.at("infrared source flux"); - auto& S_sfc_s = solution_parameters.at("solar source flux"); - - // other parameters - auto& lambda = solution_parameters.at("lambda"); - std::vector& gamma1 = solution_parameters.at("gamma1"); - std::vector& gamma2 = solution_parameters.at("gamma2"); - std::vector& gamma3 = solution_parameters.at("gamma3"); - std::vector& gamma4 = solution_parameters.at("gamma4"); - std::vector& mu = solution_parameters.at("mu"); - - auto& R_sfc = solution_parameters.at("source flux"); - - // temporary variables - T tau_cumulative = 0; - T exponential_term, denominator_term, mu_0; - - // source terms (C equations from 16, 290; eqns 23, 24) - for (std::size_t i = 0; i < number_of_columns; i++) - { - mu_0 = std::acos(solar_zenith_angles[i]); - exponential_term = omega * M_PI * R_sfc * std::exp(-(tau_cumulative - tau[i]) / mu_0); - denominator_term = (lambda * lambda - 1 / (mu_0 * mu_0)); - tau_cumulative += tau[i]; - - S_sfc_i[i] = R_sfc * mu_0 * std::exp(-tau_cumulative / mu_0); - S_sfc_s[i] = M_PI * R_sfc; - C_downwelling[i] = exponential_term * (((gamma1[i] + 1) / mu_0) * gamma4[i] + gamma2[i] * gamma3[i]); - C_upwelling[i] = exponential_term * (((gamma1[i] - 1) / mu_0) * gamma3[i] + gamma4[i] * gamma2[i]); - } - } - } - - template - void AssembleTridiagonalMatrix( - std::size_t number_of_layers, - const std::map> solution_parameters, - const std::map solver_parameters, - const TridiagonalMatrix& coeffcient_matrix) - { - // get linear system size - std::size_t matrix_size = 2 * number_of_layers; - { - // LEFT HAND SIDE 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_; - - // extract internal variables to build the matrix - const std::vector& e1 = solution_parameters.at("e1"); - const std::vector& e2 = solution_parameters.at("e2"); - const std::vector& e3 = solution_parameters.at("e3"); - const std::vector& e4 = solution_parameters.at("e4"); - - // extract surface reflectivity - const T& R_sfc = solver_parameters.at("Surface Reflectivity"); - - // first row - upper_diagonal.front() = 0; - main_diagonal.front() = e1.front(); - lower_diagonal.front() = -e2.front(); - - // odd rows - for (std::size_t n = 1; n < matrix_size - 1; n += 2) - { - upper_diagonal[n] = e2[n + 1] * e1[n] - e3[n] * e4[n + 1]; - main_diagonal[n] = e2[n] * e2[n + 1] - e3[n] * e4[n + 1]; - lower_diagonal[n] = e3[n] * e4[n + 1] - e1[n + 1] * e2[n + 1]; - } - - // even rows - for (std::size_t n = 2; n < matrix_size - 2; n += 2) - { - 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 - lower_diagonal.back() = e1.back() - R_sfc * e3.back(); - main_diagonal.back() = e2.back() - R_sfc * e4.back(); - upper_diagonal.back() = 0; - } - } - - template - void AssembleCoeffcientVector( - std::size_t number_of_layers, - const std::map> solution_parameters, - const std::map> source_terms, - const std::map solver_parameters, - std::vector& coeffcient_vector) - { - // get linear system size - std::size_t matrix_size = 2 * number_of_layers; - { - // extract internal variables to build the matrix - const std::vector& e1 = solution_parameters.at("e1"); - const std::vector& e2 = solution_parameters.at("e2"); - const std::vector& e3 = solution_parameters.at("e3"); - const std::vector& e4 = solution_parameters.at("e4"); - - // extract surface reflectivity and flux source - const T& R_sfc = solver_parameters.at("Surface Reflectivity"); - const T& f_0 = solver_parameters.at("source flux"); - - // extract source terms - const auto& C_upwelling = source_terms.at("C_upwelling"); - const auto& C_downwelling = source_terms.at("C_downwelling"); - - // first row - coeffcient_vector.front() = f_0 - C_downwelling.front(); - - // odd rows - for (std::size_t n = 1; n < matrix_size - 1; n += 2) - { - coeffcient_vector[n] = - e3[n] * (C_upwelling.start() - C_upwelling[n]) + e1[n] * (C_downwelling[n] - C_downwelling.start()); - } - - // even rows - for (std::size_t n = 2; n < matrix_size - 2; n += 2) - { - coeffcient_vector[n] = - e2[n + 1] * (C_upwelling.start() - C_upwelling[n]) + e4[n + 1] * (C_downwelling.start() - C_downwelling[n]); - } - - // last row - } - } - - 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 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 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) - { - // Solve the radiative transfer equation. - // - // [DEV NOTES] This is a placeholder for the actual implementation. - // The spherical geometry argument of the original solver was left out - // until we determine whether it needs to be an object or just a set of functions. - // - // Things that will change from the original solver: - // 1. All variables will be in SI units. Some of the original solver's - // variables were in non-SI units. - // 2. We will be solving for collections of columns. The original solver - // was for a single column. - // 3. The variable naming and source-code documentation will be improved. - const std::size_t number_of_columns = solar_zenith_angles.size(); - const auto& vertical_grid = grids.at("altitude [m]"); - const auto& wavelength_grid = grids.at("wavelength [m]"); - - // Check for consistency between the grids and profiles. - assert(vertical_grid.NumberOfColumns() == number_of_columns); - assert(wavelength_grid.NumberOfColumns() == 1); - - // internal solver variables - Array2D solution_parameters; - Array2D simulation_parameters; - - // tridiagonal system variables - TridiagonalMatrix coeffcient_matrix; - std::vector coeffcient_vector; - - tuvx::InitializeVariables( - solar_zenith_angles, grids, profiles, accumulated_radiator_state); - - ApproximationFunction(accumulated_radiator_state, solar_zenith_angles, simulation_parameters); - - tuvx::AssembleTridiagonalSystem( - solar_zenith_angles, grids, profiles, solution_parameters, coeffcient_matrix, coeffcient_vector); - - tuvx::Solve(coeffcient_matrix, coeffcient_vector); - - tuvx::ComputeRadiationField( - solar_zenith_angles, grids, profiles, solution_parameters, coeffcient_matrix, coeffcient_vector, radiation_field); - } - } // namespace tuvx diff --git a/include/tuvx/radiative_transfer/solvers/delta_eddington.hpp b/include/tuvx/radiative_transfer/solvers/delta_eddington.hpp index 35950e1f..4d306f23 100644 --- a/include/tuvx/radiative_transfer/solvers/delta_eddington.hpp +++ b/include/tuvx/radiative_transfer/solvers/delta_eddington.hpp @@ -53,13 +53,6 @@ namespace tuvx 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, @@ -72,31 +65,21 @@ namespace tuvx 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); + template + void AssembleTridiagonalMatrix( + std::size_t number_of_layers, + const std::map> solution_parameters, + const std::map solver_parameters, + const TridiagonalMatrix& coeffcient_matrix); + + template + void AssembleCoeffcientVector( + std::size_t number_of_layers, + const std::map> solution_parameters, + const std::map> source_terms, + const std::map solver_parameters, + 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, @@ -114,5 +97,8 @@ namespace tuvx /// @brief Compute the Eddington parameters template - void DeltaEddingtonApproximation(const RadiatorStatePolicy& radiator_state, const ArrayPolicy& parameters); + void DeltaEddingtonApproximation(const RadiatorStatePolicy& radiator_state, const ArrayPolicy& parotherameters); + + /// @brief Load system + }; // namespace tuvx diff --git a/include/tuvx/radiative_transfer/solvers/initialize_solver.hh b/include/tuvx/radiative_transfer/solvers/initialize_solver.hh new file mode 100644 index 00000000..1f3ea750 --- /dev/null +++ b/include/tuvx/radiative_transfer/solvers/initialize_solver.hh @@ -0,0 +1,129 @@ + +// 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 + +#include +#include + +namespace tuvx +{ + + template< + typename T, + typename GridPolicy, + typename ProfilePolicy, + typename RadiatorStatePolicy, + typename RadiationFieldPolicy, + typename SourceFunctionPolicy, + typename ArrayPolicy> + inline void InitializeVariables( + const std::vector& solar_zenith_angles, + const std::map& grids, + const std::map& profiles, + const std::map& solver_parameters, + const std::map& solution_parameters, + std::map& source_functions, + const RadiatorStatePolicy& accumulated_radiator_states) + { + // determine number of layers + const std::size_t number_of_columns = solar_zenith_angles.size(); + const auto& vertical_grid = grids.at("altitude [m]"); + const auto& wavelength_grid = grids.at("wavelength [m]"); + + // Check for consistency between the grids and profiles. + assert(vertical_grid.NumberOfColumns() == number_of_columns); + assert(wavelength_grid.NumberOfColumns() == 1); + + // Radiator state variables + ArrayPolicy& tau = accumulated_radiator_states.optical_depth_; + ArrayPolicy& omega = accumulated_radiator_states.single_scattering_albedo_; + ArrayPolicy& g = accumulated_radiator_states.assymetry_parameter_; + } + + template + void ScaleVariables(std::map>& solver_parameters) + { + // solver parameters (environment related variables) + std::vector& solar_zenith_angles = solver_parameters.at("Solar Zenith Angles"); + std::vector& tau = solver_parameters.at("Optical Depth"); + std::vector& g = solver_parameters.at("Assymetry Parameter"); + std::vector& omega = solver_parameters.at("Single Scattering Albedo"); + std::size_t number_of_columns = solar_zenith_angles.size(); + + // Delta scaling + T f; + for (std::size_t i = 0; i < number_of_columns; i++) + { + f = omega[i] * omega[i]; + g[i] = (g[i] - f) / (1 - f); + omega[i] = (1 - f) * omega[i] / (1 - omega[i] * f); + tau[i] = (1 - omega[i] * f) * tau[i]; + } + + // TODO slant optical depth computation + for (auto& tau_n : tau) + { + tau_n = tau_n; + } + } + + template + void BuildSourceFunctions( + std::map> solver_parameters, + std::map> solution_parameters, + std::map> source_functions) + { + // Source terms (C1 and C2 from the paper) + auto& C_upwelling = source_functions.at("C_upwelling"); + auto& C_downwelling = source_functions.at("C_downwelling"); + + // solver parameters (environment related variables) + std::vector& solar_zenith_angles = solver_parameters.at("Solar Zenith Angles"); + std::vector& tau = solver_parameters.at("Optical Depth"); + std::vector& omega = solver_parameters.at("Assymetry Parameter"); + std::size_t number_of_columns = solar_zenith_angles.size(); + + // parameters used to compute the soulution + std::vector& lambda = solution_parameters.at("lambda"); + std::vector& gamma1 = solution_parameters.at("gamma1"); + std::vector& gamma2 = solution_parameters.at("gamma2"); + std::vector& gamma3 = solution_parameters.at("gamma3"); + std::vector& gamma4 = solution_parameters.at("gamma4"); + std::vector& mu = solution_parameters.at("mu"); + + // solution parameters + auto& S_sfc_i = solution_parameters.at("Infrared Source Flux"); + auto& S_sfc_s = solution_parameters.at("Solar Source Flux"); + auto& R_sfc = solution_parameters.at("source flux"); + + // temporary variables + T tau_cumulative = 0; + T exponential_term, denominator_term, mu_0; + + // source terms (C equations from 16, 290; eqns 23, 24) + for (std::size_t i = 0; i < number_of_columns; i++) + { + mu_0 = std::acos(solar_zenith_angles[i]); + denominator_term = (lambda * lambda - 1 / (mu_0 * mu_0)); + tau_cumulative += tau[i]; + + S_sfc_i[i] = R_sfc * mu_0 * std::exp(-tau_cumulative / mu_0); + S_sfc_s[i] = M_PI * R_sfc; + + C_downwelling[i] = [&i, &mu_0](T tau) -> T + { + T exponential_term = omega * M_PI * R_sfc * std::exp(-(tau_cumulative - tau) / mu_0); + return exponential_term * (((gamma1[i] + 1) / mu_0) * gamma4[i] + gamma2[i] * gamma3[i]); + }; + C_upwelling[i] = [&i, &mu_0](T tau) -> T + { + T exponential_term = omega * M_PI * R_sfc * std::exp(-(tau_cumulative - tau) / mu_0); + return exponential_term * (((gamma1[i] - 1) / mu_0) * gamma3[i] + gamma4[i] * gamma2[i]); + }; + } + } +} // namespace tuvx diff --git a/include/tuvx/radiative_transfer/solvers/solve.hh b/include/tuvx/radiative_transfer/solvers/solve.hh new file mode 100644 index 00000000..205fcfd8 --- /dev/null +++ b/include/tuvx/radiative_transfer/solvers/solve.hh @@ -0,0 +1,74 @@ + + +// 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 + +#include + +namespace tuvx +{ + template< + typename T, + typename ArrayPolicy, + typename GridPolicy, + typename ProfilePolicy, + typename RadiatorStatePolicy, + typename RadiationFieldPolicy> + 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) + { + // Solve the radiative transfer equation. + // + // Things that will change from the original solver: + // 1. All variables will be in SI units. Some of the original solver's + // variables were in non-SI units. + // 2. We will be solving for collections of columns. The original solver + // was for a single column. + // 3. The variable naming and source-code documentation will be improved. + const std::size_t number_of_columns = solar_zenith_angles.size(); + const auto& vertical_grid = grids.at("altitude [m]"); + const auto& wavelength_grid = grids.at("wavelength [m]"); + + // Check for consistency between the grids and profiles. + assert(vertical_grid.NumberOfColumns() == number_of_columns); + assert(wavelength_grid.NumberOfColumns() == 1); + + // tridiagonal system variables + TridiagonalMatrix coeffcient_matrix(number_of_columns, 0); + std::vector coeffcient_vector(number_of_columns, 0); + + // internal solver variables + std::map> solution_parameters; + std::map> simulation_parameters; + std::map> source_functions; + for (std::size_t i = 0; i < wavelength_grid.NumberOfColumns(); i++) + { + solution_parameters.at("Single Scattering Albedo") = accumulated_radiator_state.single_scattering_albedo_[0][i]; + solution_parameters.at("Optical Depth") = accumulated_radiator_state.optical_depth_[0][i]; + solution_parameters.at("Assymetry Parameter") = accumulated_radiator_state.assymetry_parameter_[0][i]; + } + + tuvx::InitializeVariables( + solar_zenith_angles, grids, profiles, accumulated_radiator_state); + + // delta eddington approximation or quadrature + ApproximationFunction(accumulated_radiator_state, solar_zenith_angles, simulation_parameters); + + tuvx::AssembleTridiagonalMatrix( + solar_zenith_angles, grids, profiles, solution_parameters, coeffcient_matrix, coeffcient_vector); + + tuvx::Solve(coeffcient_matrix, coeffcient_vector); + + tuvx::ComputeRadiationField( + solar_zenith_angles, grids, profiles, solution_parameters, coeffcient_matrix, coeffcient_vector, radiation_field); + } +} // namespace tuvx