problem-2 file was erased by mistake in an earlier PR

Six_1g_spherical_benchmarks/Problem_2.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 """
-Problem 1
+Problem 2
 Chang, B. (2021a, February 19).
 Six 1D polar transport test problems for the deterministic and monte-carlo method.
 LLNL-TR-819680
@@ -24,62 +24,27 @@
     from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSolver
     from pyopensn.fieldfunc import FieldFunctionInterpolationVolume
     from pyopensn.logvol import SphereLogicalVolume, BooleanLogicalVolume
+    from pyopensn.source import VolumetricSource
 
 if __name__ == "__main__":
 
     meshgen = FromFileMeshGenerator(
-        filename="./meshes/single_sphere.msh",
+        filename="./meshes/two_region_sphere.msh",
         partitioner=PETScGraphPartitioner(type='parmetis'),
     )
     grid = meshgen.Execute()
 
-    # Get measured volumes
-    volumes_per_block = grid.ComputeVolumePerBlockID()
-
-    # Define radii per block-ID
-    radii = {1: 1.0}
-
-    # Sort block IDs by increasing radius
-    block_ids = sorted(volumes_per_block.keys())
-    if rank == 0:
-        print(block_ids)
-
-    # Check for missing radii
-    missing = [blk for blk in block_ids if blk not in radii]
-    if missing:
-        raise KeyError(f"No radius provided for block IDs: {missing}")
-
-    # Compute exact volumes
-    exact_volumes_per_block = {}
-    prev_R = 0.0
-    for blk in block_ids:
-        R = radii[blk]
-        exact_volumes_per_block[blk] = (4.0 / 3.0) * np.pi * (R ** 3 - prev_R ** 3)
-        prev_R = R
-
-    # Build the ratios array
-    ratios = np.array([
-        exact_volumes_per_block[blk] / volumes_per_block[blk]
-        for blk in block_ids
-    ])
-
-    if rank == 0:
-        for idx, blk in enumerate(block_ids):
-            print(f"Block {blk}: measured = {volumes_per_block[blk]:.11e}, "
-                  f"exact = {exact_volumes_per_block[blk]:.11e}, "
-                  f"ratio = {ratios[idx]:.6f}")
-
-    # Create and modify XS
+    # create XS
     num_groups = 1
-    xs_mat = MultiGroupXS()
-    sigt = 1.0
+    xs_void = MultiGroupXS()
+    sigt = 0.0
     c = 0.0
-    xs_mat.CreateSimpleOneGroup(sigma_t=sigt, c=c)
+    xs_void.CreateSimpleOneGroup(sigma_t=sigt, c=c)
 
-    # Boundary source
-    bsrc = [4.0 * np.pi]
+    # volumetric src
+    mg_src = VolumetricSource(block_ids=[1], group_strength=[1.])
 
-    # Angular quadrature
+    # angular quadrature
     pquad = GLCProductQuadrature3DXYZ(n_polar=2, n_azimuthal=4, scattering_order=0)
 
     # Solver
@@ -95,33 +60,42 @@
                 "angle_aggregation_num_subsets": 1,
                 "l_max_its": 500,
                 "l_abs_tol": 1.0e-6,
-                "gmres_restart_interval": 30,
+                "gmres_restart_interval": 30
             },
         ],
         xs_map=[
-            {"block_ids": [1], "xs": xs_mat},
+            {"block_ids": [1, 2], "xs": xs_void},
         ],
         options={
             "scattering_order": 0,
+            "volumetric_sources": [mg_src],
             "boundary_conditions": [
-                {"name": "xmin", "type": "isotropic", "group_strength": bsrc},
-                {"name": "xmax", "type": "isotropic", "group_strength": bsrc},
-                {"name": "ymin", "type": "isotropic", "group_strength": bsrc},
-                {"name": "ymax", "type": "isotropic", "group_strength": bsrc},
-                {"name": "zmin", "type": "isotropic", "group_strength": bsrc},
-                {"name": "zmax", "type": "isotropic", "group_strength": bsrc},
-            ],
+                {"name": "xmin", "type": "vacuum"},
+                {"name": "xmax", "type": "vacuum"},
+                {"name": "ymin", "type": "vacuum"},
+                {"name": "ymax", "type": "vacuum"},
+                {"name": "zmin", "type": "vacuum"},
+                {"name": "zmax", "type": "vacuum"}
+            ]
         },
+
     )
 
     ss_solver = SteadyStateSolver(lbs_problem=phys)
     ss_solver.Initialize()
     ss_solver.Execute()
 
-    # Export
+    # export
     fflist = phys.GetFieldFunctions()
-
-    # Compute the average flux over a logical volume
+    """
+        vtk_basename = "Problem_2_"
+        # export only the scalar flux of group g
+        FieldFunctionGridBased.ExportMultipleToVTK(
+            [fflist[g] for g in range(num_groups)],
+            vtk_basename
+        )
+    """
+    # compute the average flux over a logical volume
     def average_flx_logvol(logvol):
         ffvol = FieldFunctionInterpolationVolume()
         ffvol.SetOperationType("avg")
@@ -132,51 +106,44 @@ def average_flx_logvol(logvol):
         avgval = ffvol.GetValue()
         return avgval
 
-    # Compute the average flux over spherical shells
+    # compute the average flux over spherical shells
     def compute_average_flux(N_logvols, r_max):
         r_vals = np.linspace(0, r_max, N_logvols + 1)
         avg_flx = np.zeros(N_logvols)
         for i in range(N_logvols):
             if i != 0:
                 inner_vol = SphereLogicalVolume(r=r_vals[i])
                 outer_vol = SphereLogicalVolume(r=r_vals[i + 1])
-                logvol = BooleanLogicalVolume(parts=[
-                    {"op": True, "lv": outer_vol},
-                    {"op": False, "lv": inner_vol},
-                ])
+                logvol = BooleanLogicalVolume(parts=[{"op": True, "lv": outer_vol},
+                                                     {"op": False, "lv": inner_vol}])
             else:
                 logvol = SphereLogicalVolume(r=r_vals[i + 1])
             avg_flx[i] = average_flx_logvol(logvol)
         return r_vals, avg_flx
 
-    # Perform the calculation
+    # perform the calculation
     N_logvols = 15
     r_vals, avg_flx = compute_average_flux(N_logvols, 1)
     if rank == 0:
         for r, v in zip(r_vals[1:], avg_flx):
             print("end-radius: {:.2f}, avg-value: {:.6f}".format(r, v))
 
-    # Compute the analytical answer
-    from scipy.special import exp1
-
-    def E2(x):
-        return np.exp(-x) - x * exp1(x)
-
-    def get_phi(r, psi, r_max, sig):
-        term1 = (np.exp(-sig * (r_max - r)) - np.exp(-sig * (r_max + r))) / sig
-        term2 = (r + r_max) * E2(sig * (r_max - r))
-        term3 = (r_max - r) * E2(sig * (r_max + r))
-        phi = psi / (2 * r) * (term1 + term2 - term3)
+    # compute the analytical answer
+    def get_phi(r, q, a):
+        a2 = a**2
+        r2 = r**2
+        phi = q * (a + (a2 - r2) / (2 * r) * np.log((a + r) / np.abs(r - a)))
         return phi
 
     if rank == 0:
-        r_max = radii[1]
+
+        q = 0.5
+        a = 0.5
         eps = 1e-3
-        psi = 2 * np.pi
-        r_fine = np.linspace(0.0 + eps, r_max - eps, 100)
+        r_fine = np.linspace(0. + eps, 1.0 - eps, 100)
         exact_phi = np.zeros_like(r_fine)
         for i, r in enumerate(r_fine):
-            exact_phi[i] = get_phi(r, psi, r_max, sigt)
+            exact_phi[i] = get_phi(r, q, a)
 
         import matplotlib.pyplot as plt
 
@@ -189,5 +156,6 @@ def get_phi(r, psi, r_max, sig):
         plt.title("Flux as a function of radius")
         plt.grid(True)
         plt.legend()
-        plt.savefig("Problem_1.png", dpi=300, bbox_inches='tight')
+        # Save the figure as a PNG file
+        plt.savefig("Problem_2.png", dpi=300, bbox_inches='tight')
         plt.show()