Skip to content

Commit

Permalink
generating source functions
Browse files Browse the repository at this point in the history
  • Loading branch information
AdityaDendukuri committed Sep 9, 2024
1 parent bb939ff commit 11a904d
Show file tree
Hide file tree
Showing 4 changed files with 253 additions and 153 deletions.
93 changes: 54 additions & 39 deletions include/tuvx/radiative_transfer/solvers/delta_eddington.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,28 @@
namespace tuvx
{

/// @brief Data structure holding approximation variables for the two stream approximation. These variables are used for
/// computing the source functions and assembling the tridiagonal system. The different types of approximations can be
/// found in Table 1 of the paper. For now, only the Eddington approximation is implemented see
template<typename ArrayPolicy>
struct ApproximationVariables
{
// Two stream approximation coeffcients
ArrayPolicy gamma1_;
ArrayPolicy gamma2_;
ArrayPolicy gamma3_;
ArrayPolicy gamma4_;
ArrayPolicy mu_;
ArrayPolicy lambda_;
ArrayPolicy gamma_;
};

/// @brief Solve Function that solves the radiative transfer equation using two stream approximation
/// The function first initializes the solver variables (delta scaling, source function definitions), assembles and solves
/// the tridiagonal linear system and computes the radiation field.
///
/// @param solar_zenith_angles Zenith angles for each column
/// @param grids Domain grids - wavelength and altitude
template<
typename T,
typename ArrayPolicy,
Expand All @@ -34,56 +56,49 @@ namespace tuvx
const RadiatorStatePolicy& accumulated_radiator_states,
RadiationFieldPolicy& radiation_field);

/// @brief Compute the Eddington parameters
template<typename T, typename RadiatorStatePolicy, typename ArrayPolicy>
void DeltaEddingtonApproximation(const RadiatorStatePolicy& radiator_state, const ArrayPolicy& parameters);

template<typename T, typename GridPolicy, typename ArrayPolicy, typename RadiatorStatePolicy>
/// @brief Initialize the solver variables
template<typename T, typename GridPolicy, typename RadiatorStatePolicy, typename SourceFunctionPolicy>
void InitializeSolver(
const std::vector<T>& solar_zenith_angles,
const std::map<std::string, GridPolicy>& grids,
const std::map<std::string, std::vector<T>>& solver_variables,
const RadiatorStatePolicy& accumulated_radiator_states);
RadiatorStatePolicy& accumulated_radiator_states,
SourceFunctionPolicy& source_functions);

template<typename T, typename GridPolicy, typename ArrayPolicy, typename SourceFunctionPolicy>
void BuildSourceFunctions(
GridPolicy grids,
/// @brief Compute the Eddington parameters
template<typename T, typename ArrayPolicy, typename RadiatorStatePolicy>
void EddingtonApproximation(
const RadiatorStatePolicy& accumulated_radiator_states,
const std::vector<T>& solar_zenith_angles,
const std::map<std::string, ArrayPolicy>& solver_variables,
std::map<std::string, SourceFunctionPolicy> source_functions);

template<typename T, typename ArrayPolicy, typename GridPolicy>
void
ScaleVariables(GridPolicy grids, std::vector<T> solar_zenith_angles, std::map<std::string, ArrayPolicy>& solver_variables);
ApproximationVariables<ArrayPolicy>& approximation_variables);

/*
template<typename T>
void AssembleTridiagonalMatrix(
std::size_t number_of_layers,
const std::map<std::string, T> solver_variables,
TridiagonalMatrix<T>& coeffcient_matrix);
template<typename T>
void AssembleCoeffcientVector(
std::size_t number_of_layers,
const std::map<std::string, std::vector<T>> solver_variables,
std::map<std::string, std::function<T(T)>> source_functions,
std::vector<T>& coeffcient_vector);
template<typename T, typename GridPolicy, typename RadiatorStatePolicy>
void ScaleVariables(
const std::map<std::string, GridPolicy>& grids,
const std::vector<T>& solar_zenith_angles,
RadiatorStatePolicy& radiator_state);

template<typename T, typename RadiationFieldPolicy>
void ComputeRadiationField(
/// @brief BuildSourceFunctions Source functions $C_1$ and $C_2$ functions from the paper
template<
typename T,
typename ArrayPolicy,
typename GridPolicy,
typename RadiatorStatePolicy,
typename SourceFunctionPolicy>
void BuildSourceFunctions(
const std::size_t& number_of_columns,
const std::size_t& number_of_wavelengths,
const std::size_t& number_of_layers,
const std::vector<T>& solar_zenith_angles,
std::map<std::string, std::vector<T>> solver_variables,
const TridiagonalMatrix<T>& coeffcient_matrix,
const std::vector<T>& coeffcient_vector,
const RadiationFieldPolicy& radiation_field);
*/
const RadiatorStatePolicy& accumulated_radiator_states,
const ApproximationVariables<ArrayPolicy> approximation_variables,
const ArrayPolicy& surface_reflectivity,
std::map<std::string, SourceFunctionPolicy> source_functions);

}; // namespace tuvx
#include "delta_eddington.inl";
#include "initialize_solver.inl";

#include "delta_eddington.inl"
#include "initialize_solver.inl"
/*
include "delta_eddington.inl";
include "initialize_solver.inl";
include "assemble_tridiagonal_system.inl";
include "compute_radiation_field.inl";
Expand Down
69 changes: 41 additions & 28 deletions include/tuvx/radiative_transfer/solvers/delta_eddington.inl
Original file line number Diff line number Diff line change
@@ -1,43 +1,56 @@
// Copyright (C) 2023-2024 National Center for Atmospheric Research
// SPDX-License-Identifier: Apache-2.0

namespace tuvx
{
template<typename T, typename ArrayPolicy, typename RadiatorStatePolicy>
void DeltaEddingtonApproximation(
inline void EddingtonApproximation(
const std::size_t& number_of_columns,
const std::size_t& number_of_wavelengths,
const std::size_t& number_of_layers,
const RadiatorStatePolicy& accumulated_radiator_states,
const std::map<std::string, std::vector<T>> solution_parameters,
const std::vector<T> solar_zenith_angles)
const std::vector<T>& solar_zenith_angles,
ApproximationVariables<ArrayPolicy>& approximation_variables)
{
// allocate memory for the approximation parameters
approximation_variables.gamma1_ = tuvx::Array3D<double>(number_of_layers, number_of_wavelengths, number_of_columns);
approximation_variables.gamma2_ = tuvx::Array3D<double>(number_of_layers, number_of_wavelengths, number_of_columns);
approximation_variables.gamma3_ = tuvx::Array3D<double>(number_of_layers, number_of_wavelengths, number_of_columns);
approximation_variables.gamma4_ = tuvx::Array3D<double>(number_of_layers, number_of_wavelengths, number_of_columns);
approximation_variables.mu_ = tuvx::Array3D<double>(number_of_layers, number_of_wavelengths, number_of_columns);
approximation_variables.lambda_ = tuvx::Array3D<double>(number_of_layers, number_of_wavelengths, number_of_columns);
approximation_variables.gamma_ = tuvx::Array3D<double>(number_of_layers, number_of_wavelengths, number_of_columns);

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_;
// unpack optical properties from radiator
auto& omega = accumulated_radiator_states.optical_depth_;
auto& g = accumulated_radiator_states.single_scattering_albedo_;
auto& tau = accumulated_radiator_states.asymmetry_parameter_;

// delta eddington parameters
std::vector<T>& gamma1 = solution_parameters.at("gamma1");
std::vector<T>& gamma2 = solution_parameters.at("gamma2");
std::vector<T>& gamma3 = solution_parameters.at("gamma3");
std::vector<T>& gamma4 = solution_parameters.at("gamma4");
std::vector<T>& mu = solution_parameters.at("mu");
std::vector<T>& lambda = solution_parameters.at("lambda");
std::vector<T>& gamma = solution_parameters.at("Gamma");
// unpack delta eddington variables
auto& gamma1 = approximation_variables.gamma1_;
auto& gamma2 = approximation_variables.gamma2_;
auto& gamma3 = approximation_variables.gamma3_;
auto& gamma4 = approximation_variables.gamma4_;
auto& lambda = approximation_variables.lambda_;
auto& gamma = approximation_variables.gamma_;
auto& mu = approximation_variables.mu_;

// simulation parameters
T mu_0;
for (std::size_t i = 0; i < number_of_columns; i++)
for (std::size_t i = 0; i < number_of_layers; 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;
gamma4[i] = 1 - gamma3[i];
lambda[i] = std::sqrt(gamma1[i] * gamma1[i] - gamma2[i] * gamma2[i]);
gamma[i] = (gamma1[i] - lambda[i]) / gamma2[i];
mu = (T)0.5;
for (std::size_t j = 0; j < number_of_wavelengths; j++)
{
for (std::size_t k = 0; k < number_of_columns; k++)
{
mu_0 = std::acos(solar_zenith_angles[i]);
gamma1(i, j, k) = 7 - omega(i, j, k) * (4 + 3 * g(i, j, k));
gamma2(i, j, k) = -(1 - omega(i, j, k) * (4 - 3 * g(i, j, k))) / 4;
gamma3(i, j, k) = (2 - 3 * g(i, j, k) * mu_0) / 4;
gamma4(i, j, k) = 1 - gamma3(i, j, k);
lambda(i, j, k) = std::sqrt(gamma1(i, j, k) * gamma1(i, j, k) - gamma2(i, j, k) * gamma2(i, j, k));
gamma(i, j, k) = (gamma1(i, j, k) - lambda(i, j, k)) / gamma2(i, j, k);
mu(i, j, k) = 0.5;
}
}
}

}
} // namespace tuvx
125 changes: 52 additions & 73 deletions include/tuvx/radiative_transfer/solvers/initialize_solver.inl
Original file line number Diff line number Diff line change
Expand Up @@ -5,18 +5,11 @@
namespace tuvx
{

template<
typename T,
typename GridPolicy,
typename RadiatorStatePolicy,
typename RadiationFieldPolicy,
typename SourceFunctionPolicy,
typename ArrayPolicy>
inline void InitializeVariables(
template<typename T, typename GridPolicy, typename RadiatorStatePolicy, typename SourceFunctionPolicy>
inline void InitializeSolver(
const std::vector<T>& solar_zenith_angles,
const std::map<std::string, GridPolicy>& grids,
const RadiatorStatePolicy& accumulated_radiator_states,
std::map<std::string, ArrayPolicy>& solver_variables,
RadiatorStatePolicy& accumulated_radiator_states,
std::map<std::string, SourceFunctionPolicy>& source_functions)
{
// determine number of layers
Expand All @@ -26,34 +19,33 @@ namespace tuvx

// Check for consistency between the grids and profiles.
assert(vertical_grid.NumberOfColumns() == number_of_columns);
assert(wavelength_grid.NumberOfColumns() == 1);

// Scale variables for the Delta-Eddington approximation
ScaleVariables(grids, solar_zenith_angles, solver_variables);
ScaleVariables(grids, solar_zenith_angles, accumulated_radiator_states);

// Generate functions that define the source variables (C functions from the paper)
BuildSourceFunctions(grids, solar_zenith_angles, solver_variables, source_functions);
// BuildSourceFunctions(grids, solar_zenith_angles, solver_variables, source_functions);
}

template<typename T, typename ArrayPolicy, typename GridPolicy>
template<typename T, typename GridPolicy, typename RadiatorStatePolicy>
inline void ScaleVariables(
GridPolicy grids,
std::vector<T> solar_zenith_angles,
std::map<std::string, std::vector<T>>& solver_variables)
const std::map<std::string, GridPolicy>& grids,
const std::vector<T>& solar_zenith_angles,
RadiatorStatePolicy& radiator_state)
{
// solver parameters (environment related variables)
std::vector<T>& tau = solver_variables.at("Optical Depth");
std::vector<T>& g = solver_variables.at("Assymetry Parameter");
std::vector<T>& omega = solver_variables.at("Single Scattering Albedo");
auto& tau = radiator_state.optical_depth_;
auto& g = radiator_state.asymmetry_parameter_;
auto& omega = radiator_state.single_scattering_albedo_;

// grid and dimensions
// grids
const auto& vertical_grid = grids.at("altitude [m]");
const auto& wavelength_grid = grids.at("wavelength [m]");

// dimensions
const std::size_t number_of_columns = solar_zenith_angles.size();
const std::size_t number_of_wavelengths = wavelength_grid.NumberOfColumns();
const std::size_t number_of_layers = vertical_grid.NumberOfColumns();
const std::size_t& number_of_columns = solar_zenith_angles.size();
const std::size_t& number_of_wavelengths = wavelength_grid.NumberOfColumns();
const std::size_t& number_of_layers = vertical_grid.NumberOfColumns();

// Delta scaling
T f;
Expand All @@ -63,10 +55,10 @@ namespace tuvx
{
for (std::size_t k = 0; k < number_of_columns; k++)
{
f = omega[i][j][j] * omega[i][j][k];
g[i][j][k] = (g[i][j][k] - f) / (1 - f);
omega[i][j][k] = (1 - f) * omega[i][j][k] / (1 - omega[i][j][k] * f);
tau[i][j][k] = (1 - omega[i][j][k] * f) * tau[i][j][k];
f = omega(i, j, k) * omega(i, j, k);
g(i, j, k) = (g(i, j, k) - f) / (1 - f);
omega(i, j, k) = (1 - f) * omega(i, j, k) / (1 - omega(i, j, k) * f);
tau(i, j, k) = (1 - omega(i, j, k) * f) * tau(i, j, k);
}
}
}
Expand All @@ -78,39 +70,29 @@ namespace tuvx
//}
}

template<typename T, typename GridPolicy, typename ArrayPolicy, typename SourceFunctionPolicy>
template<typename T, typename ArrayPolicy, typename RadiatorStatePolicy, typename SourceFunctionPolicy>
inline void BuildSourceFunctions(
GridPolicy grids,
std::vector<T> solar_zenith_angles,
std::map<std::string, std::vector<T>> solver_variables,
const std::size_t& number_of_columns,
const std::size_t& number_of_wavelengths,
const std::size_t& number_of_layers,
const std::vector<T>& solar_zenith_angles,
const RadiatorStatePolicy& accumulated_radiator_states,
const ApproximationVariables<ArrayPolicy> approximation_variables,
const ArrayPolicy& surface_reflectivity,
std::map<std::string, SourceFunctionPolicy> source_functions)
{

// grid and dimensions
const auto& vertical_grid = grids.at("altitude [m]");
const auto& wavelength_grid = grids.at("wavelength [m]");

// dimensions
const std::size_t& number_of_columns = solar_zenith_angles.size();
const std::size_t& number_of_wavelengths = wavelength_grid.NumberOfColumns();
const std::size_t& number_of_layers = vertical_grid.NumberOfColumns();

// solver parameters (environment related variables)
const auto& tau = solver_variables.at("Optical Depth");
const auto& omega = solver_variables.at("Assymetry Parameter");
const auto& tau = accumulated_radiator_states.optical_depth_;
const auto& omega = accumulated_radiator_states.asymmetry_parameter_;

// parameters used to compute the soulution
const auto& lambda = solver_variables.at("lambda");
const auto& gamma1 = solver_variables.at("gamma1");
const auto& gamma2 = solver_variables.at("gamma2");
const auto& gamma3 = solver_variables.at("gamma3");
const auto& gamma4 = solver_variables.at("gamma4");
const auto& mu = solver_variables.at("mu");

// solution parameters
auto& S_sfc_i = solver_variables.at("Infrared Source Flux");
auto& S_sfc_s = solver_variables.at("Solar Source Flux");
auto& R_sfc = solver_variables.at("source flux");
const auto& mu = approximation_variables.mu_;
const auto& lambda = approximation_variables.lambda_;
const auto& gamma1 = approximation_variables.gamma1_;
const auto& gamma2 = approximation_variables.gamma2_;
const auto& gamma3 = approximation_variables.gamma3_;
const auto& gamma4 = approximation_variables.gamma4_;
const auto& R_sfc = surface_reflectivity;

// temporary variables
T tau_cumulative = 0;
Expand All @@ -120,38 +102,35 @@ namespace tuvx
auto& C_upwelling = source_functions.at("C_upwelling");
auto& C_downwelling = source_functions.at("C_downwelling");

// source terms (C equations from 16, 290; eqns 23, 24)
for (std::size_t i = 0; i < number_of_layers; i++)
for (std::size_t i = 0; i < number_of_columns; i++)
{
mu_0 = std::acos(solar_zenith_angles[i]);
for (std::size_t j = 0; j < number_of_wavelengths; j++)
{
for (std::size_t k = 0; k < number_of_columns; k++)
for (std::size_t k = 0; k < number_of_layers; k++)
{

mu_0 = std::acos(solar_zenith_angles[i][j][k]);
denominator_term = (lambda[i][j][k] * lambda[i][j][k] - 1 / (mu_0 * mu_0));
tau_cumulative += tau[i][j][k];

S_sfc_i[i][j][k] = R_sfc[i][j][k] * mu_0 * std::exp(-tau_cumulative / mu_0);
S_sfc_s[i][j][k] = M_PI * R_sfc[i][j][k];
denominator_term = (lambda(i, j, k) * lambda(i, j, k) - 1 / (mu_0 * mu_0));
tau_cumulative += tau(i, j, k);

// Source function defined for each grid point
C_downwelling[i][j][k] = [&i, &j, &k, &mu_0](T tau) -> T
C_downwelling(i, j, k) =
[&i, &j, &k, &mu_0, &R_sfc, &gamma1, &gamma2, &gamma3, &gamma4, &lambda, &mu, &omega, &tau_cumulative](
T tau) -> T
{
T exponential_term = omega[i][j][k] * M_PI * R_sfc[i][j][k] * std::exp(-(tau_cumulative - tau) / mu_0);
return exponential_term * (((gamma1[i][j][k] + 1) / mu_0) * gamma4[i][j][k] + gamma2[i][j][k] * gamma3[i][j][k]);
T exponential_term = omega(i, j, k) * M_PI * R_sfc(i, j, k) * std::exp(-(tau_cumulative - tau) / mu_0);
return exponential_term * (((gamma1(i, j, k) + 1) / mu_0) * gamma4(i, j, k) + gamma2(i, j, k) * gamma3(i, j, k));
};

C_upwelling[i][j][k] = [&i, &j, &k, &mu_0](T tau) -> T
C_upwelling(i, j, k) =
[&i, &j, &k, &mu_0, &R_sfc, &gamma1, &gamma2, &gamma3, &gamma4, &lambda, &mu, &omega, &tau_cumulative](
T tau) -> T
{
T exponential_term = omega[i][j][k] * M_PI * R_sfc[i][j][k] * std::exp(-(tau_cumulative - tau) / mu_0);
return exponential_term * (((gamma1[i][j][k] - 1) / mu_0) * gamma3[i][j][k] + gamma4[i][j][j] * gamma2[i][j][k]);
T exponential_term = omega(i, j, k) * M_PI * R_sfc(i, j, k) * std::exp(-(tau_cumulative - tau) / mu_0);
return exponential_term * (((gamma1(i, j, k) - 1) / mu_0) * gamma3(i, j, k) + gamma4(i, j, k) * gamma2(i, j, k));
};

}
}
}

}

} // namespace tuvx
Loading

0 comments on commit 11a904d

Please sign in to comment.