25654184358

Committed 09 May 2026 03:58AM UTC coverage: 75.687%. Remained the same

Build # 25654184358

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1042 from wdhawkins/ang_quadrature_refactor

Refactoring angular quadrature classes.

Coverage Stats

230 of 278 new or added lines in 27 files covered. (82.73%)

181 existing lines in 12 files now uncovered.

22258 of 29408 relevant lines covered (75.69%)

64739266.89 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.83

/framework/math/quadratures/angular/angular_quadrature.h

// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
// SPDX-License-Identifier: MIT

#pragma once

#include "framework/data_types/vector3.h"
#include "framework/data_types/ndarray.h"
#include <memory>
#include <string>
#include <vector>

namespace opensn
{
struct QuadraturePointPhiTheta;

/// Angular quadrature type identifier.
enum class AngularQuadratureType
{
  PRODUCT_QUADRATURE = 1,
  SLDFE_SQ = 2,
  LEBEDEV_QUADRATURE = 3,
  TRIANGULAR_QUADRATURE = 4,
};

/// Method used to construct the discrete-to-moment and moment-to-discrete operators.
enum class OperatorConstructionMethod
{
  STANDARD = 0,      ///< Standard weighted quadrature projection.
  GALERKIN_ONE = 1,  ///< Galerkin method 1: D2M = (M2D)^{-1}.
  GALERKIN_THREE = 3 ///< Galerkin method 3: orthogonalized approximate harmonics.
};

/// Quadrature point in spherical coordinates.
struct QuadraturePointPhiTheta
{
  /// Azimuthal angle.
  double phi = 0.0;
  /// Polar angle.
  double theta = 0.0;

  QuadraturePointPhiTheta(const double phi, const double theta) : phi(phi), theta(theta) {}
};

/// Base class for angular quadratures used in discrete ordinates transport calculations.
class AngularQuadrature
{
public:
  /// Spherical harmonic indices for moment-to-direction mappings.
  struct HarmonicIndices
  {
    /// Degree of the spherical harmonic.
    unsigned int ell = 0;
    /// Order of the spherical harmonic.
    int m = 0;

    HarmonicIndices() = default;
    HarmonicIndices(unsigned int ell, int m) : ell(ell), m(m) {}

    bool operator==(const HarmonicIndices& other) const
    {
      return (ell == other.ell and m == other.m);
    }
  };

protected:
  explicit AngularQuadrature(
    AngularQuadratureType type,
    unsigned int dimension,
    unsigned int scattering_order,
    OperatorConstructionMethod method = OperatorConstructionMethod::STANDARD)
    : type_(type),
      dimension_(dimension),
      scattering_order_(scattering_order),
      construction_method_(method)
  {
  }

  /// Sets the operator construction method.
  void SetOperatorConstructionMethod(OperatorConstructionMethod method)
  {
    construction_method_ = method;
  }

  /// Populate the map of moment index to spherical harmonic indices.
  void MakeHarmonicIndices();
  /// Quadrature-specific order parameter used by harmonic-selection rules.
  virtual unsigned int GetQuadratureOrder() const { return 0; }
  /// Number of polar angles used by harmonic-selection rules.
  virtual unsigned int GetNumPolarAngles() const { return 0; }
  /// Number of azimuthal angles used by harmonic-selection rules.
  virtual unsigned int GetNumAzimuthalAngles() const { return 0; }

  /// Discrete-to-moment operator matrix.
  /// OpenSn storage is indexed [angle][moment].
  NDArray<double, 2> d2m_op_ = {};
  /// Moment-to-discrete operator matrix.
  /// OpenSn storage is indexed [angle][moment].
  NDArray<double, 2> m2d_op_ = {};
  /// Mapping from moment index to spherical harmonic indices.
  std::vector<HarmonicIndices> m_to_ell_em_map_ = {};
  /// Quadrature type identifier.
  AngularQuadratureType type_;
  /// Spatial dimension of the quadrature.
  unsigned int dimension_;
  /// Maximum scattering order for moment calculations.
  unsigned int scattering_order_;
  OperatorConstructionMethod construction_method_;
  /// Quadrature point abscissae in spherical coordinates.
  std::vector<QuadraturePointPhiTheta> abscissae_ = {};
  /// Quadrature weights.
  std::vector<double> weights_ = {};
  /// Quadrature point direction vectors in Cartesian coordinates.
  std::vector<Vector3> omegas_ = {};

public:
  virtual ~AngularQuadrature() = default;

  /// Compute the discrete-to-moment operator.
  void BuildDiscreteToMomentOperator();

  /// Compute the moment-to-discrete operator.
  void BuildMomentToDiscreteOperator();

  /// Gets the current operator construction method
  OperatorConstructionMethod GetOperatorConstructionMethod() const { return construction_method_; }

  /// Return a reference to the precomputed discrete-to-moment operator.
  /// The operator is accessed as (d,m), where d is the direction index and m is the moment index.
  NDArray<double, 2> const& GetDiscreteToMomentOperator() const;

  /// Return a reference to the precomputed moment-to-discrete operator.
  /// The operator is accessed as (d,m), where d is the direction index and m is the moment index.
  NDArray<double, 2> const& GetMomentToDiscreteOperator() const;

  /// Return a reference to the precomputed harmonic index map.
  const std::vector<HarmonicIndices>& GetMomentToHarmonicsIndexMap() const;

  unsigned int GetDimension() const { return dimension_; }

  unsigned int GetScatteringOrder() const { return scattering_order_; }

  unsigned int GetNumMoments() const { return m_to_ell_em_map_.size(); }

  size_t GetNumAngles() const { return omegas_.size(); }

  bool Empty() const { return GetNumAngles() == 0; }

  AngularQuadratureType GetType() const { return type_; }

  /// Return a user-facing quadrature name.
  virtual std::string GetName() const = 0;

  /// Return the sum of all quadrature weights.
  virtual double GetWeightSum() const;

  const std::vector<QuadraturePointPhiTheta>& GetAbscissae() const { return abscissae_; }

  const QuadraturePointPhiTheta& GetAbscissa(size_t angle_index) const
  {
    return abscissae_.at(angle_index);
  }

  const std::vector<double>& GetWeights() const { return weights_; }

  double GetWeight(size_t angle_index) const { return weights_.at(angle_index); }

  const std::vector<Vector3>& GetOmegas() const { return omegas_; }

  const Vector3& GetOmega(size_t angle_index) const { return omegas_.at(angle_index); }
};

} // namespace opensn

1	// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#pragma once
5
6	#include "framework/data_types/vector3.h"
7	#include "framework/data_types/ndarray.h"
8	#include <memory>
9	#include <string>
10	#include <vector>
11
12	namespace opensn
13	{
14	struct QuadraturePointPhiTheta;
15
16	/// Angular quadrature type identifier.
17	enum class AngularQuadratureType
18	{
19	PRODUCT_QUADRATURE = 1,
20	SLDFE_SQ = 2,
21	LEBEDEV_QUADRATURE = 3,
22	TRIANGULAR_QUADRATURE = 4,
23	};
24
25	/// Method used to construct the discrete-to-moment and moment-to-discrete operators.
26	enum class OperatorConstructionMethod
27	{
28	STANDARD = 0, ///< Standard weighted quadrature projection.
29	GALERKIN_ONE = 1, ///< Galerkin method 1: D2M = (M2D)^{-1}.
30	GALERKIN_THREE = 3 ///< Galerkin method 3: orthogonalized approximate harmonics.
31	};
32
33	/// Quadrature point in spherical coordinates.
34	struct QuadraturePointPhiTheta
35	{
36	/// Azimuthal angle.
37	double phi = 0.0;
38	/// Polar angle.
39	double theta = 0.0;
40
41	QuadraturePointPhiTheta(const double phi, const double theta) : phi(phi), theta(theta) {}	96,023✔
42	};
43
44	/// Base class for angular quadratures used in discrete ordinates transport calculations.
45	class AngularQuadrature
46	{
47	public:
48	/// Spherical harmonic indices for moment-to-direction mappings.
49	struct HarmonicIndices
50	{
51	/// Degree of the spherical harmonic.
52	unsigned int ell = 0;
53	/// Order of the spherical harmonic.
54	int m = 0;
55
56	HarmonicIndices() = default;
57	HarmonicIndices(unsigned int ell, int m) : ell(ell), m(m) {}	2,457✔
58
59	bool operator==(const HarmonicIndices& other) const
60	{
61	return (ell == other.ell and m == other.m);
62	}
63	};
64
65	protected:
66	explicit AngularQuadrature(	797✔
67	AngularQuadratureType type,
68	unsigned int dimension,
69	unsigned int scattering_order,
70	OperatorConstructionMethod method = OperatorConstructionMethod::STANDARD)
71	: type_(type),	797✔
72	dimension_(dimension),	797✔
73	scattering_order_(scattering_order),	797✔
74	construction_method_(method)	797✔
75	{
76	}	797✔
77
78	/// Sets the operator construction method.
79	void SetOperatorConstructionMethod(OperatorConstructionMethod method)
80	{
81	construction_method_ = method;
82	}
83
84	/// Populate the map of moment index to spherical harmonic indices.
85	void MakeHarmonicIndices();
86	/// Quadrature-specific order parameter used by harmonic-selection rules.
87	virtual unsigned int GetQuadratureOrder() const { return 0; }	4✔
88	/// Number of polar angles used by harmonic-selection rules.
89	virtual unsigned int GetNumPolarAngles() const { return 0; }
90	/// Number of azimuthal angles used by harmonic-selection rules.
91	virtual unsigned int GetNumAzimuthalAngles() const { return 0; }
92
93	/// Discrete-to-moment operator matrix.
94	/// OpenSn storage is indexed [angle][moment].
95	NDArray<double, 2> d2m_op_ = {};
96	/// Moment-to-discrete operator matrix.
97	/// OpenSn storage is indexed [angle][moment].
98	NDArray<double, 2> m2d_op_ = {};
99	/// Mapping from moment index to spherical harmonic indices.
100	std::vector<HarmonicIndices> m_to_ell_em_map_ = {};
101	/// Quadrature type identifier.
102	AngularQuadratureType type_;
103	/// Spatial dimension of the quadrature.
104	unsigned int dimension_;
105	/// Maximum scattering order for moment calculations.
106	unsigned int scattering_order_;
107	OperatorConstructionMethod construction_method_;
108	/// Quadrature point abscissae in spherical coordinates.
109	std::vector<QuadraturePointPhiTheta> abscissae_ = {};
110	/// Quadrature weights.
111	std::vector<double> weights_ = {};
112	/// Quadrature point direction vectors in Cartesian coordinates.
113	std::vector<Vector3> omegas_ = {};
114
115	public:
116	virtual ~AngularQuadrature() = default;	1,502✔
117		2✔
118	/// Compute the discrete-to-moment operator.	2✔
119	void BuildDiscreteToMomentOperator();
120
121	/// Compute the moment-to-discrete operator.
122	void BuildMomentToDiscreteOperator();
123
124	/// Gets the current operator construction method
125	OperatorConstructionMethod GetOperatorConstructionMethod() const { return construction_method_; }
126
127	/// Return a reference to the precomputed discrete-to-moment operator.
128	/// The operator is accessed as (d,m), where d is the direction index and m is the moment index.
129	NDArray<double, 2> const& GetDiscreteToMomentOperator() const;
130
131	/// Return a reference to the precomputed moment-to-discrete operator.
132	/// The operator is accessed as (d,m), where d is the direction index and m is the moment index.
133	NDArray<double, 2> const& GetMomentToDiscreteOperator() const;
134
135	/// Return a reference to the precomputed harmonic index map.
136	const std::vector<HarmonicIndices>& GetMomentToHarmonicsIndexMap() const;
137
138	unsigned int GetDimension() const { return dimension_; }
139
140	unsigned int GetScatteringOrder() const { return scattering_order_; }
141
142	unsigned int GetNumMoments() const { return m_to_ell_em_map_.size(); }
143
144	size_t GetNumAngles() const { return omegas_.size(); }	73,360✔
145
146	bool Empty() const { return GetNumAngles() == 0; }
147
148	AngularQuadratureType GetType() const { return type_; }	830✔
149
150	/// Return a user-facing quadrature name.
151	virtual std::string GetName() const = 0;
152
153	/// Return the sum of all quadrature weights.
154	virtual double GetWeightSum() const;
155
156	const std::vector<QuadraturePointPhiTheta>& GetAbscissae() const { return abscissae_; }
157
158	const QuadraturePointPhiTheta& GetAbscissa(size_t angle_index) const
159	{
160	return abscissae_.at(angle_index);
161	}
162
163	const std::vector<double>& GetWeights() const { return weights_; }
164
165	double GetWeight(size_t angle_index) const { return weights_.at(angle_index); }	2,147,483,647✔
166
167	const std::vector<Vector3>& GetOmegas() const { return omegas_; }
168
169	const Vector3& GetOmega(size_t angle_index) const { return omegas_.at(angle_index); }	2,147,483,647✔
170	};	1,062✔
171		2,791✔
172	} // namespace opensn	2,221✔

Open-Sn / opensn / 25654184358

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous