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

framework/math/quadratures/angular/triangle_quadrature.h

@@ -0,0 +1,32 @@
+// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include "framework/math/quadratures/angular/angular_quadrature.h"
+
+namespace opensn
+{
+
+class TriangleQuadrature : public AngularQuadrature
+{
+protected:
+  unsigned int method_;
+  unsigned int sn_;
+  unsigned int moments_;
+
+  void TriangleInit();
+
+  void MakeHarmonicIndices(unsigned int scattering_order, int dimension) override;
+
+public:
+  explicit TriangleQuadrature(unsigned int method, unsigned int sn, unsigned int moments = 0);
+
+  void BuildDiscreteToMomentOperator(unsigned int scattering_order, int dimension) override;
+
+  void BuildMomentToDiscreteOperator(unsigned int scattering_order, int dimension) override;
+
+  void FilterMoments(unsigned int scattering_order);
+};
+
+} // namespace opensn

lua/framework/math/quadratures/create_triangle_quadature.cc

@@ -0,0 +1,46 @@
+// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
+// SPDX-License-Identifier: MIT
+
+#include "framework/lua.h"
+#include "quadratures_lua.h"
+#include "framework/math/quadratures/angular/triangle_quadrature.h"
+#include "framework/logging/log.h"
+#include "framework/console/console.h"
+#include "framework/runtime.h"
+#include <memory>
+
+using namespace opensn;
+
+namespace opensnlua
+{
+
+RegisterLuaFunctionNamespace(CreateTriangleQuadrature, aquad, CreateTriangleQuadrature);
+
+int
+CreateTriangleQuadrature(lua_State* L)
+{
+  int num_args = lua_gettop(L);
+  // Parse argument
+  int method = lua_tonumber(L, 1);
+  int order = lua_tonumber(L, 2);
+  int moments = 0;
+  if (num_args == 3)
+    moments = lua_tonumber(L, 3);
+
+ // if (num_args < 2)
+ //   LuaPostArgAmountError("CreateTriangleQuadrature", 2, num_args);
+
+  opensn::log.Log() << "Creating Triangle Quadrature\n";
+
+  auto new_quad = std::make_shared<TriangleQuadrature>(method, order, moments);
+  opensn::angular_quadrature_stack.push_back(new_quad);
+  const size_t index = opensn::angular_quadrature_stack.size() - 1;
+  lua_pushinteger(L, static_cast<lua_Integer>(index));
+
+  opensn::log.Log() << "Created Triangle Quadrature with method = " << method << " and order "
+                    << order << std::endl;
+
+  return 1;
+}
+
+} // namespace opensnlua

lua/framework/math/quadratures/quadratures_lua.h

@@ -97,6 +97,7 @@ int CreateSphericalProductQuadrature(lua_State* L);
  * \author Jan
  */
 int CreateProductQuadrature(lua_State* L);
+int CreateTriangleQuadrature(lua_State* L);
 
 /** Creates a quadrature.
  *

test/modules/linear_boltzmann_solvers/transport_steady/tests.json

@@ -201,6 +201,36 @@
     ],
     "skip": "Working on a solution - Jan"
   },
+{
+    "file": "transport_2d_triangular_P0.lua",
+    "comment": "2D LinearBSolver Test isotropic scatter with triangular quadrature",
+    "num_procs": 1,
+    "weight_class" : "short",
+    "checks": [
+      {
+        "type": "FloatCompare",
+        "key": "max-grp0(latest)",
+        "wordnum" : 4,
+        "gold": 2.938333e+03,
+        "abs_tol": 1.0e-6
+      }
+     ]
+   },
+   {
+    "file": "transport_2d_triangular_P0_GQ3.lua",
+    "comment": "2D LinearBSolver Test isotropic scatter with triangular quadrature method 3",
+    "num_procs": 1,
+    "weight_class" : "short",
+    "checks": [
+      {
+        "type": "FloatCompare",
+        "key": "max-grp0(latest)",
+        "wordnum" : 4,
+        "gold": 2.932038,
+        "abs_tol": 1.0e-6
+      }
+     ]
+   },
   {
     "file": "transport_3d_1a_extruder.lua",
     "comment": "3D LinearBSolver Test - PWLD",
@@ -451,5 +481,21 @@
         "abs_tol": 1.0e-9
       }
     ]
-  }
+  },
+  {
+    "file": "transport_3d_triangular_P0.lua",
+    "comment": "3D LinearBSolver Test isotropic scatter with triangular quadrature",
+    "num_procs": 8,
+    "weight_class" : "intermediate",
+    "checks": [
+      {
+        "type": "FloatCompare",
+        "key": "max-grp0(latest)",
+        "wordnum" : 4,
+        "gold":  2.096991e+00,
+        "abs_tol": 1.0e-6
+      }
+     ]
+   }
+
 ]

test/modules/linear_boltzmann_solvers/transport_steady/transport_2d_triangular_P0.lua

@@ -0,0 +1,86 @@
+-- 2D LinearBSolver vacuum with isotropic source distributed and incident, triangular set 
+num_procs = 1
+
+if (check_num_procs==nil and number_of_processes ~= num_procs) then
+    log.Log(LOG_0ERROR,"Incorrect amount of processors. " ..
+            "Expected "..tostring(num_procs)..
+            ". Pass check_num_procs=false to override if possible.")
+    os.exit(false)
+end
+
+nodes = {}
+N = 65
+ds = 500.0 / N
+for i = 0, N do
+    nodes[i + 1] = i * ds
+end
+meshgen1 = mesh.OrthogonalMeshGenerator.Create({ node_sets = { nodes, nodes } })
+mesh.MeshGenerator.Execute(meshgen1)
+vol0 = logvol.RPPLogicalVolume.Create({ infx = true, infy = true, infz = true })
+mesh.SetMaterialIDFromLogicalVolume(vol0, 0)
+materials = {}
+materials[1] = mat.AddMaterial("Test Material");
+mat.AddProperty(materials[1], TRANSPORT_XSECTIONS)
+mat.AddProperty(materials[1], ISOTROPIC_MG_SOURCE)
+num_groups = 1
+mat.SetProperty(materials[1], TRANSPORT_XSECTIONS, OPENSN_XSFILE, "simple_scatter.xs")
+src = {}
+for g = 1, num_groups do
+    src[g] = 0.0
+end
+src[1] = 1.0
+mat.SetProperty(materials[1], ISOTROPIC_MG_SOURCE, FROM_ARRAY, src)
+
+method = 0
+scattering_order = 0
+sn = 8
+moments = 0
+tquad = aquad.CreateTriangleQuadrature(method, sn)
+aquad.OptimizeForPolarSymmetry(tquad,4.0*math.pi)
+
+lbs_block = {
+    num_groups = num_groups,
+    groupsets = {
+        {
+            groups_from_to = { 0, 0 },
+            angular_quadrature_handle = tquad,
+            angle_aggregation_type = "single",
+            angle_aggregation_num_subsets = 1,
+            groupset_num_subsets = 1,
+            inner_linear_method = "gmres",
+            l_abs_tol = 1.0e-8,
+            l_max_its = 300,
+            gmres_restart_interval = 50,
+        },
+    }
+}
+
+lbs_options = {
+    boundary_conditions = {
+        {
+            name = "xmin",
+            type = "isotropic",
+            group_strength = { 1.0 },
+        }
+    },
+    scattering_order = 0,
+}
+
+phys1 = lbs.DiscreteOrdinatesSolver.Create(lbs_block)
+lbs.SetOptions(phys1, lbs_options)
+ss_solver = lbs.SteadyStateSolver.Create({ lbs_solver_handle = phys1 })
+solver.Initialize(ss_solver)
+solver.Execute(ss_solver)
+fflist, count = lbs.GetScalarFieldFunctionList(phys1)
+pp1 = post.CellVolumeIntegralPostProcessor.Create
+({
+  name="max-grp0",
+  field_function = fflist[1],
+  compute_volume_average = true,
+  print_numeric_format = "scientific"
+})
+post.Execute({ pp1 })
+if master_export == nil then
+    fieldfunc.ExportToVTKMulti(fflist,"ZPhi")
+end
+

test/modules/linear_boltzmann_solvers/transport_steady/transport_2d_triangular_P0_GQ3.lua

@@ -0,0 +1,79 @@
+-- 2D LinearBSolver Vacuum with isotropic incident, triangular set, Henyey-GS Forward peaked cross-sections
+num_procs = 1
+
+if (check_num_procs==nil and number_of_processes ~= num_procs) then
+    log.Log(LOG_0ERROR,"Incorrect amount of processors. " ..
+            "Expected "..tostring(num_procs)..
+            ". Pass check_num_procs=false to override if possible.")
+    os.exit(false)
+end
+
+nodes = {}
+N = 65
+ds = 1.0 / N
+for i = 0, N do
+    nodes[i + 1] = i * ds
+end
+meshgen1 = mesh.OrthogonalMeshGenerator.Create({ node_sets = { nodes, nodes } })
+mesh.MeshGenerator.Execute(meshgen1)
+vol0 = logvol.RPPLogicalVolume.Create({ infx = true, infy = true, infz = true })
+mesh.SetMaterialIDFromLogicalVolume(vol0, 0)
+materials = {}
+materials[1] = mat.AddMaterial("Test Material");
+mat.AddProperty(materials[1], TRANSPORT_XSECTIONS)
+num_groups = 1
+mat.SetProperty(materials[1], TRANSPORT_XSECTIONS, OPENSN_XSFILE, "xs_Henyey_GS_p4.cxs")
+
+method = 3
+scattering_order = 0
+sn = 8
+moments = 0
+tquad = aquad.CreateTriangleQuadrature(method, sn,scattering_order)
+aquad.OptimizeForPolarSymmetry(tquad,4.0*math.pi)
+
+lbs_block = {
+    num_groups = num_groups,
+    groupsets = {
+        {
+            groups_from_to = { 0, 0 },
+            angular_quadrature_handle = tquad,
+            angle_aggregation_type = "single",
+            angle_aggregation_num_subsets = 1,
+            groupset_num_subsets = 1,
+            inner_linear_method = "gmres",
+            l_abs_tol = 1.0e-8,
+            l_max_its = 300,
+            gmres_restart_interval = 50,
+        },
+    }
+}
+
+lbs_options = {
+    boundary_conditions = {
+        {
+            name = "xmin",
+            type = "isotropic",
+            group_strength = { 1.0 },
+        }
+    },
+    scattering_order = 0,
+}
+
+phys1 = lbs.DiscreteOrdinatesSolver.Create(lbs_block)
+lbs.SetOptions(phys1, lbs_options)
+ss_solver = lbs.SteadyStateSolver.Create({ lbs_solver_handle = phys1 })
+solver.Initialize(ss_solver)
+solver.Execute(ss_solver)
+fflist, count = lbs.GetScalarFieldFunctionList(phys1)
+pp1 = post.CellVolumeIntegralPostProcessor.Create
+({
+  name="max-grp0",
+  field_function = fflist[1],
+  compute_volume_average = true,
+  print_numeric_format = "scientific"
+})
+post.Execute({ pp1 })
+if master_export == nil then
+    fieldfunc.ExportToVTKMulti(fflist,"ZPhi")
+end
+

test/modules/linear_boltzmann_solvers/transport_steady/transport_3d_triangular_P0.lua

@@ -0,0 +1,82 @@
+-- 3D LinearBSolver Vacuum with isotropic incident, triangular set
+num_procs = 8
+
+if (check_num_procs==nil and number_of_processes ~= num_procs) then
+    log.Log(LOG_0ERROR,"Incorrect amount of processors. " ..
+            "Expected "..tostring(num_procs)..
+            ". Pass check_num_procs=false to override if possible.")
+    os.exit(false)
+end
+
+nodes = {}
+N = 65
+ds = 10.0/N
+for i = 0, N do
+    nodes[i+1] = i*ds
+end
+meshgen1 = mesh.OrthogonalMeshGenerator.Create({ node_sets = {nodes,nodes,nodes} })
+mesh.MeshGenerator.Execute(meshgen1)
+vol0 = logvol.RPPLogicalVolume.Create({infx=true, infy=true, infz=true})
+mesh.SetMaterialIDFromLogicalVolume(vol0,0)
+materials = {}
+materials[1] = mat.AddMaterial("Test Material");
+mat.AddProperty(materials[1], TRANSPORT_XSECTIONS)
+
+num_groups = 1
+mat.SetProperty(materials[1], TRANSPORT_XSECTIONS, OPENSN_XSFILE, "simple_scatter.xs")
+method = 3
+sn = 4
+moments = 0
+tquad = aquad.CreateTriangleQuadrature(method,sn)
+
+lbs_block =
+{
+    num_groups = num_groups,
+    groupsets =
+    {
+        {
+            groups_from_to = {0,0},
+            angular_quadrature_handle = tquad,
+            angle_aggregation_type = "single",
+            angle_aggregation_num_subsets = 1,
+            groupset_num_subsets = 1,
+            inner_linear_method = "gmres",
+            l_abs_tol = 1.0e-8,
+            l_max_its = 300,
+            gmres_restart_interval = 10,
+        },
+    }
+}
+
+lbs_options =
+{
+    boundary_conditions =
+    {
+        {
+            name = "xmin",
+            type = "isotropic",
+            group_strength={1.0},
+        }
+    },
+    scattering_order = 0,
+}
+
+phys1 = lbs.DiscreteOrdinatesSolver.Create(lbs_block)
+lbs.SetOptions(phys1, lbs_options)
+ss_solver = lbs.SteadyStateSolver.Create({lbs_solver_handle = phys1})
+solver.Initialize(ss_solver)
+solver.Execute(ss_solver)
+fflist,count = lbs.GetScalarFieldFunctionList(phys1)
+pp1 = post.CellVolumeIntegralPostProcessor.Create
+({
+  name="max-grp0",
+  field_function = fflist[1],
+  compute_volume_average = true,
+  print_numeric_format = "scientific"
+})
+post.Execute({ pp1 })
+if (master_export == nil) then
+    fieldfunc.ExportToVTKMulti(fflist,"ZPhi3D")
+end
+
+