From 191fea95caec004eb8ee51925b4af0bcf6cc23db Mon Sep 17 00:00:00 2001 From: ragusa Date: Wed, 2 Jul 2025 20:07:20 -0500 Subject: [PATCH] problem-2 file was erased by mistake in an earlier PR --- Six_1g_spherical_benchmarks/Problem_2.py | 126 +++++++++-------------- 1 file changed, 47 insertions(+), 79 deletions(-) diff --git a/Six_1g_spherical_benchmarks/Problem_2.py b/Six_1g_spherical_benchmarks/Problem_2.py index 16e1ed4..3bdef60 100644 --- a/Six_1g_spherical_benchmarks/Problem_2.py +++ b/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,7 +106,7 @@ 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) @@ -140,43 +114,36 @@ def compute_average_flux(N_logvols, r_max): 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()