Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding Kobayashi Dogleg model #25

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Conversation

ninadol
Copy link

@ninadol ninadol commented Oct 30, 2024

I worked with Dr. Stephen Coleman at Radiasoft LLC to construct the Kobayashi Dogleg Benchmark for problem 3i and 3ii. Stephen formatted my initial model for this package. The relevant plots are included in the post-processing notebook.
This is novel as a benchmark because it utilizes the new MAGIC Method implementation in OpenMC. This is also novel because it is outside the Sinbad database. Stephen will connect below, and to help with any feedback you might have for us.

Removing large files

Fixing Problem_II

New results with modified Problem II
@sjcrs
Copy link

sjcrs commented Oct 30, 2024

Hi Ethan, we spoke a few months ago about this project coming along this summer. Any comments or suggested changes you might have will be handled by one or both of us. I wanted Nina to issue the pull request, though, since she did most of the hard work, and I just formatted it to conform with your other models.
There are some sizable differences in flux represented here relative to the other simulation codes, and I have shown that to the OpenMC team at ANL. If there are any upstream changes to weight window generation that change the output here, I hope that will propagate through to this model and the plots generated in postprocessing.py will update accordingly.

@eepeterson
Copy link
Owner

Thanks so much for this great addition @ninadol and @sjcrs! I will take a look at this soon!

def _parse_args():
"""Parse and return commandline arguments"""
parser = argparse.ArgumentParser()
parser.add_argument("-b", "--batches", type=int, default=20)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
parser.add_argument("-b", "--batches", type=int, default=20)
parser.add_argument("-b", "--batches", type=int, default=100)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We usually leave 100 batches as default just for convention.

parser = argparse.ArgumentParser()
parser.add_argument("-b", "--batches", type=int, default=20)
parser.add_argument("-p", "--particles", type=int, default=int(1e7))
parser.add_argument("-s", "--threads", type=int,default=5)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
parser.add_argument("-s", "--threads", type=int,default=5)
parser.add_argument("-s", "--threads", type=int)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just by convention, we prefer not to specify any default number of threads in the other benchmarks.

Comment on lines +177 to +190
Source = openmc.Material(name='source_mat')
Source.set_density('macro', 1.0)
Source.add_macroscopic('source')

Absorber = openmc.Material(name='absorber_mat')
Absorber.set_density('macro', 1.0)
Absorber.add_macroscopic('absorber')

Void = openmc.Material(name='void_mat')
Void.set_density('macro', 1.0)
Void.add_macroscopic('void')

# Instantiate a Materials collection and export to XML
materials = openmc.Materials([Source,Absorber,Void])
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Source = openmc.Material(name='source_mat')
Source.set_density('macro', 1.0)
Source.add_macroscopic('source')
Absorber = openmc.Material(name='absorber_mat')
Absorber.set_density('macro', 1.0)
Absorber.add_macroscopic('absorber')
Void = openmc.Material(name='void_mat')
Void.set_density('macro', 1.0)
Void.add_macroscopic('void')
# Instantiate a Materials collection and export to XML
materials = openmc.Materials([Source,Absorber,Void])
source = openmc.Material(name='source_mat')
source.set_density('macro', 1.0)
source.add_macroscopic('source')
absorber = openmc.Material(name='absorber_mat')
absorber.set_density('macro', 1.0)
absorber.add_macroscopic('absorber')
void = openmc.Material(name='void_mat')
void.set_density('macro', 1.0)
void.add_macroscopic('void')
# Instantiate a Materials collection and export to XML
materials = openmc.Materials([source, absorber, void])

Comment on lines +234 to +248
source_cell = openmc.Cell(fill = Source,
region = box2,
name = 'source region')

void_cell = openmc.Cell(fill = Void,
region = (segment1 | segment2 | -segment3 | -segment4),
name = 'void region')

absorber_cell = openmc.Cell(fill = Absorber,
region = box0 & ~void_cell.region & ~source_cell.region,
name = 'absorber region')


# Fill the Cell with the Material
Ucell.fill = Void
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
source_cell = openmc.Cell(fill = Source,
region = box2,
name = 'source region')
void_cell = openmc.Cell(fill = Void,
region = (segment1 | segment2 | -segment3 | -segment4),
name = 'void region')
absorber_cell = openmc.Cell(fill = Absorber,
region = box0 & ~void_cell.region & ~source_cell.region,
name = 'absorber region')
# Fill the Cell with the Material
Ucell.fill = Void
source_cell = openmc.Cell(fill = source,
region = box2,
name = 'source_cell')
void_cell = openmc.Cell(fill = void,
region = (segment1 | segment2 | -segment3 | -segment4),
name = 'void_cell')
absorber_cell = openmc.Cell(fill = absorber,
region = box0 & ~void_cell.region & ~source_cell.region,
name = 'absorber_cell')
# Fill the Cell with the Material
ucell.fill = void

Comment on lines +104 to +111
Ucell = openmc.Cell(name='base universe', region = box)

# domain: the domain for spatial homogenization aka
# averaging the cross sections over volume and flux spectrum

total = mgxs.TotalXS(domain=Ucell, energy_groups=groups)
absorption = mgxs.AbsorptionXS(domain=Ucell, energy_groups=groups)
scattering = mgxs.ScatterXS(domain=Ucell, energy_groups=groups)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not convinced by starting with a capital letter for naming the objects.

Suggested change
Ucell = openmc.Cell(name='base universe', region = box)
# domain: the domain for spatial homogenization aka
# averaging the cross sections over volume and flux spectrum
total = mgxs.TotalXS(domain=Ucell, energy_groups=groups)
absorption = mgxs.AbsorptionXS(domain=Ucell, energy_groups=groups)
scattering = mgxs.ScatterXS(domain=Ucell, energy_groups=groups)
ucell = openmc.Cell(name='base_universe', region = box)
# domain: the domain for spatial homogenization aka
# averaging the cross sections over volume and flux spectrum
total = mgxs.TotalXS(domain=ucell, energy_groups=groups)
absorption = mgxs.AbsorptionXS(domain=ucell, energy_groups=groups)
scattering = mgxs.ScatterXS(domain=ucell, energy_groups=groups)

Comment on lines +33 to +66
# generate weight windows tallies
try:
p = subprocess.Popen(["python3", "openmc_model.py", "--Problem_I","--weight_windows"])
except:
p = subprocess.Popen(["python", "openmc_model.py", "--Problem_I","--weight_windows"])

# wait for the simulation to finish
p.wait()

# run simulation without scattering
try:
p = subprocess.Popen(["python3", "openmc_model.py", "--Problem_I"])
except:
p = subprocess.Popen(["python", "openmc_model.py", "--Problem_I"])

# wait for the simulation to finish
p.wait()

try:
p = subprocess.Popen(["python3", "openmc_model.py", "--weight_windows"])
except:
p = subprocess.Popen(["python", "openmc_model.py", "--weight_windows"])

# wait for the simulation to finish
p.wait()

# run simulation with scattering
try:
p = subprocess.Popen(["python3", "openmc_model.py"])
except:
p = subprocess.Popen(["python", "openmc_model.py"])

# wait for the simulation to finish
p.wait()
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would not default to a ww generation run. Especially for the run_and_store file. I would rather provide already some ww files and leave the ww generation as optional. Or does this benchmark specifically focus on ww generation and performances?

Suggested change
# generate weight windows tallies
try:
p = subprocess.Popen(["python3", "openmc_model.py", "--Problem_I","--weight_windows"])
except:
p = subprocess.Popen(["python", "openmc_model.py", "--Problem_I","--weight_windows"])
# wait for the simulation to finish
p.wait()
# run simulation without scattering
try:
p = subprocess.Popen(["python3", "openmc_model.py", "--Problem_I"])
except:
p = subprocess.Popen(["python", "openmc_model.py", "--Problem_I"])
# wait for the simulation to finish
p.wait()
try:
p = subprocess.Popen(["python3", "openmc_model.py", "--weight_windows"])
except:
p = subprocess.Popen(["python", "openmc_model.py", "--weight_windows"])
# wait for the simulation to finish
p.wait()
# run simulation with scattering
try:
p = subprocess.Popen(["python3", "openmc_model.py"])
except:
p = subprocess.Popen(["python", "openmc_model.py"])
# wait for the simulation to finish
p.wait()
# run simulation without scattering
try:
p = subprocess.Popen(["python3", "openmc_model.py", "--Problem_I"])
except:
p = subprocess.Popen(["python", "openmc_model.py", "--Problem_I"])
# wait for the simulation to finish
p.wait()
# run simulation with scattering
try:
p = subprocess.Popen(["python3", "openmc_model.py"])
except:
p = subprocess.Popen(["python", "openmc_model.py"])
# wait for the simulation to finish
p.wait()

Comment on lines +68 to +72
print("Reading Statepoint file")
# read statepoint file
ProblemI_file = ofb.ResultsFromOpenmc('statepoint.20.h5', 'Problem_I')
# read statepoint file
ProblemII_file = ofb.ResultsFromOpenmc('statepoint.20.h5', 'Problem_II')
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
print("Reading Statepoint file")
# read statepoint file
ProblemI_file = ofb.ResultsFromOpenmc('statepoint.20.h5', 'Problem_I')
# read statepoint file
ProblemII_file = ofb.ResultsFromOpenmc('statepoint.20.h5', 'Problem_II')
print("Reading Statepoint file")
# read statepoint file
problemI_file = ofb.ResultsFromOpenmc('statepoint.100.h5', 'Problem_I')
# read statepoint file
problemII_file = ofb.ResultsFromOpenmc('statepoint.100.h5', 'Problem_II')

@SteSeg
Copy link
Collaborator

SteSeg commented Dec 12, 2024

Hi @ninadol and @sjcrs, thank you so much for contributing to this! The Kobayashi benchmark is a great addition to the repository!

I'll leave a couple of comments here that we think should help better homogenize this benchmark with the other ones, if you want to address them:

Style:

  • filenames and foldernames: we are not using capital letters
  • postprocessing.ipynb: some "/" appear in the text
  • Are the StoreBenchmarkData.ipynb and Kobayashi_3ii.ipynb notebooks strictly necessary to run the simulation and analyse resultas?

I would like to open a discussion about using the magic method in the default run. I believe we can either provide the ww file already and foget about it or maybe having the ww generation run as optional with the argparse module. Or is this benchmark specifically conceived for ww generation and performance analysis? @eepeterson what do you think?

Lastly, and most important, did you happen to have the time to explore why openmc is giving results that are so different from mcnp?

Thank you again for the great work!

@sjcrs
Copy link

sjcrs commented Dec 12, 2024

Hi @SteSeg , thank you for your review and for your suggestions. As you said, these are mostly for style conformity and therefore we have no objection to making these changes. StoreBenchmarkData.ipynb and Kobayashi_3ii.ipynb are not strictly necessary to run the simulation, but Kobayashi_3ii.ipynb is there to generate the comparison plots for analysis. Would you prefer that these be consolidated into a postprocessing.ipynb notebook?

The goal of @ninadol 's work was to benchmark the weight windows generation using the MAGIC method, and the goal with this problem was to demonstrate its utility. I had thought about structuring this so that it first produced an analog run with large error bars before moving on to the MAGIC method weight windows generation, but was concerned about producing too many comparisons.

I think it's reasonable to provide a ww file and optinally generate a new one via argparse, but I'll have to think about how (if at all) we provide the user a way to compare the user-generated file to the 'ground truth' one.

Finally, I have spoken to @jtramm about the differences between the results here, because i know he has been looking at the Kobayashi dogleg for a different variance reduction feature. Movement on this PR will encourage us to meet and puzzle this one out.

@sjcrs
Copy link

sjcrs commented Dec 12, 2024

Actually, I'd like to correct myself, Kobayashi_3ii.ipynb is not required and should not have been committed. The analysis is in postprocessing.ipynb.

@SteSeg
Copy link
Collaborator

SteSeg commented Dec 23, 2024

Thank you for the clarification. I think the best way to go is to provide an original ww file, then add the ww generation as optional and have just two notebooks. One for the actual results analysis (postprocessing.ipynb) and the other for the ww performance analysis. In this way one the postprocessing.ipynb doesn't get affected by the choice of running or not the ww generation.

Anyways, I think we should dig more in the results differences, usually openmc and mcnp results are not so far apart. I'll run the benchmark soon just to add another couple of eyes to search for the problem.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants