Extract eigenvalue source from statepoint

User Support
1

Context

Eigenvalue-mode depletion calculation.

Objective

Restart the depletion from last statepoint to extend time steps starting with the converged source from the last step of the previous calculation.

The idea is perhaps I can reduce the number of required inactive batches to speed up transport calculation.

Problem

Apparently, a Source.from_statepoint() method used to exist but I don’t see anything comparable for the replacement IndependentSource in the documentation (I’m on the develop branch).

2

I don’t think this is going to work. It appears based on initial keff iteration results that each timestep in a depletion calculation starts with the original source specified in settings, not the converged source from the previous time step.

I suppose I could dump the converged source from an initial eigenvalue calculation (always run anyway), and use that as a FileSource for the depletion calculation.

Just need to figure out how to dump the converged source or extract it from a statepoint.

3

Here’s what I came up with to extract source from previous calculation and set it up as a FileSource. Basic idea anyway. I’m no python expert, and I had to fight the linter the whole way.

Show code 
def get_source_from_case_dir(case_dir: Path) -> openmc.FileSource:
    import numpy as np
    from openmc.particle_type import ParticleType
    from openmc.source import ParticleList, SourceParticle, write_source_file

    if not os.path.exists(case_dir / "settings.xml"):
        raise ValueError("settings.xml not found")
    settings = openmc.Settings.from_xml(case_dir / "settings.xml")
    num_batches = settings.batches
    if num_batches is None:
        raise ValueError("batches not found")
    if num_batches == 0:
        raise ValueError("batches is 0")

    # Get the source_bank from statepoint file
    st_file = case_dir / f"statepoint.{num_batches}.h5"
    if not os.path.exists(st_file):
        raise ValueError(f"{st_file} not found")

    # Disable autolink else fails
    st = openmc.StatePoint(st_file, autolink=False)
    source_bank: np.ndarray = st.source  # type: ignore
    assert source_bank is not None
    # Convert to ParticleList (see ParticleList.from_hdf5)
    particles = ParticleList(
        [SourceParticle(*params, ParticleType(int(pdg))) for *params, pdg in source_bank]  # type: ignore
    )

    source_file = "./source.h5"
    write_source_file(particles, source_file)

    return openmc.FileSource(source_file)
4

If running MPI, restrict the source file write to a single rank:

    from mpi4py import MPI
    # Only write source file on rank 0
    source_file = "./source.h5"
    comm = MPI.COMM_WORLD
    if comm.rank == 0:
        write_source_file(particles, source_file)
    comm.Barrier()