Segfault with "munmap_chunk(): invalid pointer" when using multiple groups per groupset in fixed-source problem

[hityyds]

Hi everyone,

I’m trying to use OpenSN to solve a fixed-source problem involving multiple groupsets, where each groupset contains four energy groups. However, as soon as the computation begins for groupset 0, I encounter the following error:

[0]  ********** Solving groupset 0 with CLASSIC_RICHARDSON.
[0]  
[0]  Quadrature number of angles: 96
[0]  Groups 0 3
[0]  
*** Error in python3': free(): invalid pointer: 0x0000555653f8b110 ***
======= Backtrace: =========
/lib64/libc.so.6(+0x81329)[0x2aefde2cf329]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZNK6opensn17AngularQuadrature27GetMomentToDiscreteOperatorEv+0x98)[0x2af04ef58bf8]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZN6opensn13AAHSweepChunk5SweepERNS_8AngleSetE+0x99)[0x2af04f0f6bc9]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZN6opensn12AAH_AngleSet15AngleSetAdvanceERNS_10SweepChunkENS_14AngleSetStatusE+0x1c2)[0x2af04f0cc642]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZN6opensn14SweepScheduler15ScheduleAlgoDOGERNS_10SweepChunkE+0x73)[0x2af04f0e7693]

It seems that the issue is related to AngularQuadrature. I tried reducing the number of angles, but I still got a similar error.

[0]  ********** Solving groupset 0 with CLASSIC_RICHARDSON.
[0]  
[0]  Quadrature number of angles: 96
[0]  Groups 0 3
[0]  
*** Error in `python3': munmap_chunk(): invalid pointer: 0x0000555d4ec2be40 ***
*** Error in `python3': free(): invalid size: 0x000055d0bfc2a5d0 ***
======= Backtrace: =========
/lib64/libc.so.6(+0x7f474)[0x2af142322474]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZN6opensn13AAHSweepChunk5SweepERNS_8AngleSetE+0x1a0c)[0x2af1b30f853c]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZN6opensn12AAH_AngleSet15AngleSetAdvanceERNS_10SweepChunkENS_14AngleSetStatusE+0x1c2)[0x2af1b30cc642]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZN6opensn14SweepScheduler15ScheduleAlgoDOGERNS_10SweepChunkE+0x73)[0x2af1b30e7693]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZN6opensn14SweepScheduler5SweepEv+0x6c)[0x2af1b30e7ecc]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZN6opensn15SweepWGSContext29ApplyInverseTransportOperatorENS_11SourceFlagsE+0x8f)[0x2af1b30c4acf]
/public3/home/a6s001028/OpenSN_install/opensn/build/lib.linux-x86_64-cpython-310/pyopensn/libopensn.so.0(_ZN6opensn17ClassicRichardson5SolveEv+0x337)[0x2af1b311d117]

Interestingly, I had previously tested the same case with only one group per groupset, and the simulation ran successfully without any issues. The only difference between the two runs is the number of groups included in each groupset.

Below is the solver configuration section from my Python script:

 # Angular Quadrature
    pquad = GLCProductQuadrature3DXYZ(6, 12)

    # LBS block option
    # inner_linear_method = "petsc_gmres"
    inner_linear_method = "classic_richardson"
    phys = DiscreteOrdinatesProblem(
        mesh=grid,
        num_groups=num_groups,
        groupsets=[
            {
                "groups_from_to": (0, 3),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
                #"gmres_restart_interval": 15, 

                # "apply_wgdsa": True,
                # "wgdsa_l_abs_tol": 1.0e-6,
                # "wgdsa_verbose": True,
            },
            {
                "groups_from_to": (4, 7),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
                #"gmres_restart_interval": 15, 

                # "apply_wgdsa": True,
                # "wgdsa_l_abs_tol": 1.0e-6,
                # "wgdsa_verbose": True,
            },

Does anyone have an idea what might be causing this segmentation fault? Could it be a bug or a limitation related to multi-group groupsets?

Any guidance or suggestions would be greatly appreciated!

Thanks in advance.

[wdhawkins]

Is this with the latest version of the code? Can you send your input so that I can reproduce the problem?

[hityyds]

Thank you for your reply!

The code version I’m using should be the Git version from June 6, 2025.
The original error occurred on a supercomputing platform with a relatively complex problem.

However, I just managed to reproduce the same error on my local machine using a very simple mesh file, with the same energy group and angular quadrature settings.
Below is the Python script I used locally to reproduce the issue:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Infinite, 1-group, pure absorber with balance

import os
import sys
import numpy as np

if "opensn_console" not in globals():
    from mpi4py import MPI
    size = MPI.COMM_WORLD.size
    rank = MPI.COMM_WORLD.rank
    from pyopensn.mesh import FromFileMeshGenerator, PETScGraphPartitioner
    from pyopensn.xs import MultiGroupXS
    from pyopensn.source import VolumetricSource
    from pyopensn.aquad import GLCProductQuadrature3DXYZ
    from pyopensn.solver import DiscreteOrdinatesProblem, SteadyStateSolver
    from pyopensn.logvol import RPPLogicalVolume
    from pyopensn.fieldfunc import FieldFunctionGridBased, FieldFunctionInterpolationVolume

if __name__ == "__main__":

    mesh_file = "rect.msh"
    meshgen = FromFileMeshGenerator(
        filename=mesh_file, 
        partitioner=PETScGraphPartitioner(type='ptscotch')
        )
    grid = meshgen.Execute()

    # Add cross sections to materials
    num_groups = 47
    num_mats = 2
    scat_order = 3
    xss = []
    for m in range(num_mats):
        xss.append(MultiGroupXS())
        xss[m].LoadFromOpenSn("mat/mat_%d.xs" % m)
        
    # Create map between materials and mesh blocks
    xs_map = []

    xs_map.append({"block_ids": [1], "xs": xss[0]})
    xs_map.append({"block_ids": [2], "xs": xss[1]})

    mg_srcs = []
    group_strength = [0.0 for _ in range(num_groups)]
    group_strength[0] = 1000.0
    mg_srcs.append(VolumetricSource(block_ids=[1], group_strength=group_strength))

    # Angular Quadrature
    pquad = GLCProductQuadrature3DXYZ(6, 12)

    # LBS block option
    # inner_linear_method = "petsc_gmres"
    inner_linear_method = "classic_richardson"
    phys = DiscreteOrdinatesProblem(
        mesh=grid,
        num_groups=num_groups,
        groupsets=[
            {
                "groups_from_to": (0, 3),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (4, 7),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (8, 11),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (12, 15),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (16, 19),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (20, 23),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (24, 27),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (28,31),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (32, 35),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (36, 39),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (40, 43),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },
            {
                "groups_from_to": (44, 46),
                "angular_quadrature": pquad,
                "angle_aggregation_type": "single",
                "angle_aggregation_num_subsets": 1,
                "inner_linear_method": inner_linear_method,
                "l_abs_tol": 1.0e-3,
                "l_max_its": 5000,
            },

        ],
        xs_map=xs_map,

        options={
            "boundary_conditions": [
                {"name": "xmin", "type": "reflecting"},
                {"name": "ymin", "type": "reflecting"}
            ],
            "volumetric_sources": mg_srcs,
            "scattering_order": scat_order,
        }
        #sweep_type = "AAH"
    )

    # Initialize and execute solver
    ss_solver = SteadyStateSolver(lbs_problem=phys)
    ss_solver.Initialize()
    ss_solver.Execute()

The mesh file was generated using Gmsh, and the corresponding .geo file is as follows:

//+
SetFactory("OpenCASCADE");
//+
Cylinder(1) = {0, 0, -10, 0, 0, 20, 5, 2*Pi};
//+
Box(2) = {-20, -20, -20, 40, 40, 40};
//+
Physical Volume("core", 1) = {1};

//+
BooleanDifference{ Volume{2}; Delete; }{ Volume{1}; }
Physical Volume("Ref", 2) = {2};

MeshSize { : } = 2;

// Using HXT for parallel meshing
Mesh.Algorithm3D = 10;
Mesh 3;
Coherence Mesh;
OptimizeMesh "Netgen";

The material file is quite long, so I don’t think it’s suitable to paste it directly here.
I'm not sure if you already have a 47-group material file on hand — it shouldn't be critical for reproducing the error.

I also noticed that if I comment out the boundary condition settings, the error does not occur.

Hope this helps narrow down the issue!

Let me know if you need anything else.

[wdhawkins]

Thanks! I'll see if I can reproduce the problem. I'll let you know what I find.

[wdhawkins]

Are you running this with the console or with the Python module?

[wdhawkins]

I'm able to run this successfully with the latest version of the code. I only tried the console app, so let me know if you're using the Python module. To run with the latest version, you'll need to change your quadrature object slightly:
pquad = GLCProductQuadrature3DXYZ(n_polar=6, n_azimuthal=12, scattering_order=scat_order)

If you are using multiple groupsets, you probably need to adjust the convergence parameters for the across-groupset solver to get the correct answer. Without configuring the AGS solver, you will not converge the upscattering between groupsets.

[hityyds]

Thanks for the follow-up!

Yes, I'm using the Python module.

I also have an additional question:
Does OpenSN only support boundary conditions on the six axis-aligned faces (X/Y/Z_min/max)? From what I’ve seen in the source code, it looks like OpenSN determines these faces by geometric calculation, and does not read user-defined boundary face sets from the mesh file. Is that correct?

Thanks again for your help!

[wdhawkins]

I'll double check that I can run the problem with the Python module, but I'm confident it will work. There was a PR that was recently merged that fixed an issue with reflecting boundary conditions with multiple groupsets, and I suspect it resolved this issue. I'll wait to close this until I've tested the Python module.

You are correct that we only support boundary conditions on the six axis-aligned faces. We do have plans to support arbitrary face-based boundary conditions, but we don't currently have anyone working on that capability.