Adding in Triangular Quadrature set and GQ1-3 methods with accompanying Lua #209

emersonshands · 2024-04-09T16:44:09Z

This merge adds in the Gauss-Chebyshev triangular quadrature set with lua calls for the set. Additionally it has methods for forming the D&M matrices seperately from the standard angular quadrature methods currently in the repository. This includes the standard Sn method, termed method 0, uses the standard full rank harmonics mapping and forms both the D&M seperately. The additional first and second Galerkin Quadrature methods are mentioned in Comparison of Two Galerkin Quadrature Methods, Morel 2017, Ref.. Both form the full range D or M matrices then invert them to obtain the accompanying matrix, which in the case of a lower range expansion requires that both matrices be cut down to the required expansion. The third method, GQ3, has not been published yet and is focused on building the D&M matrices up to the desired degree by use of the modified Gram-Schmidt method. All three GQ methods are purpose built for anisotropic problems and are agnostic to the quadrature set with appropriate harmonic choices.

Further verification is needed against an anisotropic benchmark problem from Ref..

wdhawkins · 2024-04-09T17:17:09Z

Could you also add some regression tests to this PR? They don't need to be complicated. Just simple tests to verify that the methods are working correctly.

emersonshands · 2024-04-09T17:26:42Z

Sure I have a 3d and 2d I can add-- I'll add them to the test/modules/LBS/transport_steady.

andrsd · 2024-04-09T17:37:02Z

You will need to do a few things first:

Rebase your work on the current main. This is to remove those 2 merge commits. Do not do this when developing OpenSn.
Then, you need to properly format your code via clang-format - we supply a clang-format settings that should be picked up by clang-format automatically as long as you do that from the repo directory. You can use interactive rebase and format individual patches as you go if you want to preserve the branch history. Do not create a clang-format patch!
You are missing SPDX tags in your new files. This change happened this morning so you were not aware of it, see more info

andrsd · 2024-04-09T19:06:18Z

You still have 2 merge commits in your branch, so you did not rebase anything. You need to rebase on top of upstream main. You need to do something like this:

git fetch upstream main
git rebase --interactive upstream/main

upstream here is a remote that points to git@github.com:Open-Sn/opensn.git. If you don't have it (you can check that via git remote -v), then create it via git remote add upstream git@github.com:Open-Sn/opensn.git. This needs to happen before the fetch command.

When you do the 2nd command, you will see something like this in a text editor:

pick 8a7f719b Adding in triangular quadrature set to branch. New files added and Lua for quadratures modified.
pick 61c2abbc Changed the configure_dep python file to correctly work with folder structure by touching the local files.
pick 70949b9b Added in previous functionality to form a full 3d set that can be reduced to 2d using OptimizeForPolarSymmetry.
pick d120f1b3 Added in 3d formation for set and fixed harmonic mapping. Testing still required to verify against benchmark.
pick 33985015 Rebased and added SPX labeling.

Then, you will remove line 2 ("Changed the configure_dep"). This commit creates a conflict that you would have to resolve. However, you should not be modifying that script in your PR anyway. Your PR is about a quadrature. That's why you will get rid of it. If those changes are needed, then they can come in a separate PR. Save and quit the editor. Git will rebase your branch and then you will have to force push it into your fork.

I advised you to NOT create a clang-format patch, but you did that anyway, so 👎 .

Last advise for future. Do not develop on main branch. Create a feature branch (for example via git checkout -b my-feature-branch) and do all your work there. It will make your life much simpler.

test/modules/linear_boltzmann_solvers/transport_steady/transport_2d_triangular_P0.lua

ragusa · 2024-04-27T00:27:24Z

framework/math/quadratures/angular/triangle_quadrature.h

+  unsigned int sn_;
+  unsigned int moments_;
+
+  void TriangleInit();


check with @wdhawkins but we may need fo explicit names such as TriangleQuadratueInitialize

Just go with Initialize for right now.

ragusa · 2024-04-27T00:27:55Z

framework/math/quadratures/angular/triangle_quadrature.h

+  void BuildDiscreteToMomentOperator(unsigned int scattering_order, int dimension) override;
+
+  void BuildMomentToDiscreteOperator(unsigned int scattering_order, int dimension) override;


check with @wdhawkins if dim will remain int or not

It should be fine for now. It may change later one, but we can leave it for this commit.

ragusa · 2024-04-27T00:28:46Z

framework/math/quadratures/angular/triangle_quadrature.cc

+}
+
+void
+TriangleQuadrature::TriangleInit()


see comment in .h file

ragusa · 2024-04-27T00:30:12Z

framework/math/quadratures/angular/triangle_quadrature.cc

+  const auto GAUSS_LEG_QUAD = GaussLegendreQuadrature(sn_);
+
+  // Formulate the triangular quadrature
+  VecDbl newZi, newWeights;


check with @wdhawkins if we should name this Zi (technically you meant Xi) or Mu (for consistency).

This should probably be mu.

framework/math/quadratures/angular/triangle_quadrature.cc

ragusa · 2024-04-27T00:42:16Z

@emersonshands also, the golden values for the tests are missing (see how to modify the .json file in the proper folder)

…a for quadratures modified.

…uced to 2d using OptimizeForPolarSymmetry.

…ll required to verify against benchmark.

Need verification due to python errors.

Refactoring variable names for standards

Additional test for GQ3 added, along with associated Henyey-G.S. forward peaked cross section file.

Waiting for response and gold file values.

wdhawkins · 2024-05-03T20:50:54Z

framework/math/quadratures/angular/triangle_quadrature.cc

+  const auto GAUSS_LEG_QUAD = GaussLegendreQuadrature(sn_);
+
+  // Formulate the triangular quadrature
+  VecDbl newZi, newWeights;


I'm guessing VecDbl is a typedef for std::vector<double>? If so, please use std::vector<double>.

wdhawkins · 2024-05-03T20:53:16Z

framework/math/quadratures/angular/triangle_quadrature.cc

+  {
+    log.Log0Error() << "TriangleQuarature: " << method_ << " is not a valid method.\n"
+                    << "Method must be 0, 1, 2, or 3." << std::endl;
+    Exit(510);


This should throw an std::invalid_argument instead of Exit(510).

wdhawkins · 2024-05-03T20:54:07Z

framework/math/quadratures/angular/triangle_quadrature.cc

+
+  // Formulate the triangular quadrature
+  VecDbl newZi, newWeights;
+  newZi.reserve(GAUSS_LEG_QUAD.qpoints.size());


No caps in GAUSS_LEG_QUAD or at least not all caps.

wdhawkins · 2024-05-03T20:55:15Z

framework/math/quadratures/angular/triangle_quadrature.cc

+
+  int num_div = 1;
+  int weightPos = 0;
+  for (const auto& QUAD_POINT : GAUSS_LEG_QUAD.qpoints)


No all caps for QUAD_POINT.

wdhawkins · 2024-05-03T20:55:51Z

framework/math/quadratures/angular/triangle_quadrature.cc

+    for (int v = 0; v < num_div; ++v)
+    {
+      double new_z_value = abs(QUAD_POINT.x);
+      double phi = deltaVPhi / 2.0 + (double)v * deltaVPhi;


Please use static_cast<double>.

wdhawkins · 2024-05-03T20:56:17Z

framework/math/quadratures/angular/triangle_quadrature.cc

+      double phi = deltaVPhi / 2.0 + (double)v * deltaVPhi;
+      double theta = acos(new_z_value);
+      double sinTheta = sqrt(1 - new_z_value * new_z_value);
+      //      double weightCurrent = GAUSS_LEG_QUAD.weights[weightPos] / (num_div);


Remove this commented line if it's not necessary.

wdhawkins · 2024-05-03T20:57:58Z

framework/math/quadratures/angular/triangle_quadrature.cc

+      new_omega.x = sinTheta * cos(phi);
+      new_omega.y = sinTheta * sin(phi);
+      new_omega.z = new_z_value;
+      weights_.push_back(GAUSS_LEG_QUAD.weights[weightPos] * (M_PI / (2.0 * num_div)));


Since you know the number of weights, omegas, and abscissae, reserve their size before the loop.

wdhawkins · 2024-05-03T21:01:04Z

framework/math/quadratures/angular/triangle_quadrature.cc

+    new_omega.x = omegas_[point].x; // omegas[l].x*xsign;
+    new_omega.y = omegas_[point].y; // omegas[l].y*ysign;
+    new_omega.z = -omegas_[point].z;
+    double new_z_value = -omegas_[point].z;


Instead of declaring these new variables, you can use the expressions inline in the lines below.

wdhawkins · 2024-05-07T18:58:11Z

framework/math/quadratures/angular/triangle_quadrature.cc

+  // We need to fix the normalization to be 4*pi over the full sphere
+  double normal = 4.0 * M_PI;
+
+  const size_t NUM_DIRS = omegas_.size();


lower case for NUM_DIRS

wdhawkins · 2024-05-07T19:00:03Z

framework/math/quadratures/angular/triangle_quadrature.cc

+  // need to cut after the formation.
+  if (m2d_op_built_ and d2m_op_built_ and moments_ >= 0 and !m_to_ell_em_map_.empty())
+  {
+    int s_order = static_cast<int>(scattering_order);


Make the parameter an int and remove this cast.

wdhawkins · 2024-05-07T19:00:55Z

framework/math/quadratures/angular/triangle_quadrature.cc

+
+    auto m2d_transposed = m2d_op_;
+    m_to_ell_em_map_.resize(moments_to_keep);
+    std::vector<std::vector<double>> m2dworking;


Reserve m2dworking and d2mworking since you know the size.

wdhawkins · 2024-05-07T19:02:42Z

framework/math/quadratures/angular/triangle_quadrature.cc

+void
+TriangleQuadrature::MakeHarmonicIndices(unsigned int scattering_order, int dimension)
+{
+  int L = static_cast<int>(sn_);


Change the unsigned int's in the header to int's and remove these casts.

wdhawkins · 2024-05-07T19:03:08Z

framework/math/quadratures/angular/triangle_quadrature.cc

+}
+
+void
+TriangleQuadrature::MakeHarmonicIndices(unsigned int scattering_order, int dimension)


Change scattering_order to an int and remove the cast below.

wdhawkins · 2024-05-07T19:04:50Z

framework/math/quadratures/angular/triangle_quadrature.cc

+  {
+    d2m_op_.clear();
+
+    const size_t NUM_ANGLES = abscissae_.size();


NUM_ANGLES -> num_angles. Remove caps throughout the file.

wdhawkins assigned emersonshands Apr 9, 2024

wdhawkins requested review from andrsd, ragusa and zhardy-lanl April 9, 2024 17:13

wdhawkins added the enhancement New feature or request label Apr 9, 2024

emersonshands force-pushed the main branch from 3398501 to 0d66eff Compare April 9, 2024 19:12

emersonshands marked this pull request as draft April 9, 2024 19:55

emersonshands force-pushed the main branch 2 times, most recently from 49f2047 to cf7b2e5 Compare April 19, 2024 21:15

emersonshands marked this pull request as ready for review April 19, 2024 21:22

emersonshands force-pushed the main branch from f72fa42 to 1bdbc32 Compare April 25, 2024 20:54

ragusa reviewed Apr 27, 2024

View reviewed changes

Emerson-Sh and others added 11 commits May 1, 2024 10:00

Adding in triangular quadrature set to branch. New files added and Lu…

5f79b5b

…a for quadratures modified.

Added in previous functionality to form a full 3d set that can be red…

0493d5f

…uced to 2d using OptimizeForPolarSymmetry.

Added in 3d formation for set and fixed harmonic mapping. Testing sti…

0e3dd41

…ll required to verify against benchmark.

Rebased and added SPX labeling.

f5b1458

Fix to 3d and filter moments method

0ccefff

Added 2 P0 Tests

a1754fa

Need verification due to python errors.

Fix for P0 3d case, missed in 9b07752

32831cb

Refactoring for code standards

dc32293

Refactoring variable names for standards

Fixes to tests on output and fflist. Running on cluster.

6c97cb7

Adding in GQ3 Test with Henyey truncated to p0.

2f7e2bd

Additional test for GQ3 added, along with associated Henyey-G.S. forward peaked cross section file.

Added in syntax fixes for tests

c18cb6e

Waiting for response and gold file values.

emersonshands force-pushed the main branch from 1bdbc32 to c18cb6e Compare May 1, 2024 15:49

Fixes for running cluster. Test speed mods.

3a66227

Emerson-Sh added 2 commits May 1, 2024 11:58

Gold file values and xs fix.

116c521

Naming issue with cross-section file.

80ad98c

wdhawkins reviewed May 3, 2024

View reviewed changes

wdhawkins reviewed May 7, 2024

View reviewed changes

		void BuildDiscreteToMomentOperator(unsigned int scattering_order, int dimension) override;

		void BuildMomentToDiscreteOperator(unsigned int scattering_order, int dimension) override;

Adding in Triangular Quadrature set and GQ1-3 methods with accompanying Lua #209

Are you sure you want to change the base?

Adding in Triangular Quadrature set and GQ1-3 methods with accompanying Lua #209

Uh oh!

Conversation

emersonshands commented Apr 9, 2024

Uh oh!

wdhawkins commented Apr 9, 2024

Uh oh!

emersonshands commented Apr 9, 2024

Uh oh!

andrsd commented Apr 9, 2024

Uh oh!

andrsd commented Apr 9, 2024

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

ragusa commented Apr 27, 2024

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wdhawkins May 3, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wdhawkins May 7, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wdhawkins May 7, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

wdhawkins May 3, 2024 •

edited

Loading

wdhawkins May 7, 2024 •

edited

Loading

wdhawkins May 7, 2024 •

edited

Loading