From 32cc4c8d434018d9d14fa90bd09b628bf668701f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 10 Feb 2020 16:01:30 -0600 Subject: [PATCH 01/59] Support timestep units in Integrator.__init__ --- openmc/deplete/abc.py | 82 ++++++++++++++++++++------- openmc/deplete/integrators.py | 104 ++++++++++++++++++++++++++-------- 2 files changed, 142 insertions(+), 44 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 10976dd0ed..a9640852cc 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -31,6 +31,9 @@ __all__ = [ "Integrator", "SIIntegrator", "DepSystemSolver"] +_SECONDS_PER_DAY = 24*60*60 +_SECONDS_PER_YEAR = 365.25*24*60*60 + OperatorResult = namedtuple('OperatorResult', ['k', 'rates']) OperatorResult.__doc__ = """\ Result of applying transport operator @@ -597,9 +600,11 @@ class Integrator(ABC): ---------- operator : openmc.deplete.TransportOperator Operator to perform transport simulations - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -612,6 +617,11 @@ class Integrator(ABC): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). Attributes ---------- @@ -625,7 +635,8 @@ class Integrator(ABC): Power of the reactor in [W] for each interval in :attr:`timesteps` """ - def __init__(self, operator, timesteps, power=None, power_density=None): + def __init__(self, operator, timesteps, power=None, power_density=None, + timestep_units='s'): # Check number of stages previously used if operator.prev_res is not None: res = operator.prev_res[-1] @@ -638,27 +649,51 @@ class Integrator(ABC): self._num_stages)) self.operator = operator self.chain = operator.chain - if not isinstance(timesteps, Iterable): - self.timesteps = [timesteps] - else: - self.timesteps = timesteps + + # Determine power and normalize units to W + mass = operator.heavy_metal if power is None: if power_density is None: raise ValueError("Either power or power density must be set") if not isinstance(power_density, Iterable): - power = power_density * operator.heavy_metal + power = power_density * mass else: - power = [p * operator.heavy_metal for p in power_density] - + power = [p*mass for p in power_density] if not isinstance(power, Iterable): # Ensure that power is single value if that is the case - power = [power] * len(self.timesteps) - elif len(power) != len(self.timesteps): + power = [power] * len(timesteps) + + if len(power) != len(timesteps): raise ValueError( "Number of time steps != number of powers. {} vs {}".format( - len(self.timesteps), len(power))) + len(timesteps), len(power))) - self.power = power + # Get list of times / units + if isinstance(timesteps[0], Iterable): + times, units = zip(*timesteps) + else: + times = timesteps + units = [timestep_units] * len(timesteps) + + # Determine number of seconds for each timestep + seconds = [] + for time, unit, watts in zip(times, units, power): + if unit == 's': + seconds.append(time) + elif unit in ('d', 'day'): + seconds.append(time*_SECONDS_PER_DAY) + elif unit in ('a', 'yr', 'year'): + seconds.append(time*_SECONDS_PER_YEAR) + elif unit.lower() == 'mwd/kg': + watt_days_per_kg = 1e6*time + kilograms = 1e-3*mass + days = watt_days_per_kg * kilograms / watts + seconds.append(days*_SECONDS_PER_DAY) + else: + raise ValueError("Invalid timestep unit: {}".format(unit)) + + self.timesteps = asarray(seconds) + self.power = asarray(power) @abstractmethod def __call__(self, conc, rates, dt, power, i): @@ -772,9 +807,11 @@ class SIIntegrator(Integrator): ---------- operator : openmc.deplete.TransportOperator The operator object to simulate on. - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -787,6 +824,11 @@ class SIIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 @@ -805,10 +847,10 @@ class SIIntegrator(Integrator): Number of stochastic iterations per depletion interval """ def __init__(self, operator, timesteps, power=None, power_density=None, - n_steps=10): + timestep_units='s', n_steps=10): check_type("n_steps", n_steps, Integral) check_greater_than("n_steps", n_steps, 0) - super().__init__(operator, timesteps, power, power_density) + super().__init__(operator, timesteps, power, power_density, timestep_units) self.n_steps = n_steps def _get_bos_data_from_operator(self, step_index, step_power, bos_conc): diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index 67106aa3e1..84f838cb7a 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -31,9 +31,11 @@ class PredictorIntegrator(Integrator): ---------- operator : openmc.deplete.TransportOperator Operator to perform transport simulations - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -46,6 +48,11 @@ class PredictorIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). Attributes ---------- @@ -113,9 +120,11 @@ class CECMIntegrator(Integrator): ---------- operator : openmc.deplete.TransportOperator Operator to perform transport simulations - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -128,6 +137,11 @@ class CECMIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). Attributes ---------- @@ -203,9 +217,11 @@ class CF4Integrator(Integrator): ---------- operator : openmc.deplete.TransportOperator Operator to perform transport simulations - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -218,6 +234,11 @@ class CF4Integrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). Attributes ---------- @@ -310,9 +331,11 @@ class CELIIntegrator(Integrator): ---------- operator : openmc.deplete.TransportOperator Operator to perform transport simulations - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -325,6 +348,11 @@ class CELIIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). Attributes ---------- @@ -404,9 +432,11 @@ class EPCRK4Integrator(Integrator): ---------- operator : openmc.deplete.TransportOperator Operator to perform transport simulations - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -419,6 +449,11 @@ class EPCRK4Integrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). Attributes ---------- @@ -518,9 +553,11 @@ class LEQIIntegrator(Integrator): ---------- operator : openmc.deplete.TransportOperator Operator to perform transport simulations - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -533,6 +570,11 @@ class LEQIIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). Attributes ---------- @@ -629,9 +671,11 @@ class SICELIIntegrator(SIIntegrator): ---------- operator : openmc.deplete.TransportOperator The operator object to simulate on. - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -644,6 +688,11 @@ class SICELIIntegrator(SIIntegrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 @@ -730,9 +779,11 @@ class SILEQIIntegrator(SIIntegrator): ---------- operator : openmc.deplete.TransportOperator The operator object to simulate on. - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not - cumulative. + timesteps : iterable of float or iterable of tuple + Array of timesteps. Note that values are not cumulative. The units are + specified by the `timestep_units` argument when `timesteps` is an + iterable of float. Alternatively, units can be specified for each step + by passing an iterable of (value, unit) tuples. power : float or iterable of float, optional Power of the reactor in [W]. A single value indicates that the power is constant over all timesteps. An iterable @@ -745,6 +796,11 @@ class SILEQIIntegrator(SIIntegrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. + timestep_units : {'s', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that + the values are given in burnup (MW-d of energy deposited per kilogram + of heavy metal). n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 From 8771987bf40ad036107419acf9ea948fc3a1919a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 11 Feb 2020 06:44:05 -0600 Subject: [PATCH 02/59] Add unit test checking depletion timestep units for all integrators --- tests/unit_tests/test_deplete_integrator.py | 69 ++++++++++++++++++++- 1 file changed, 68 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 47c9690ae9..9d5766fd5c 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -7,6 +7,7 @@ will be left unimplemented and testing will be done via regression. """ import copy +from random import uniform from unittest.mock import MagicMock import numpy as np @@ -15,11 +16,24 @@ import pytest from openmc.deplete import ( ReactionRates, Results, ResultsList, comm, OperatorResult, - PredictorIntegrator, SICELIIntegrator) + PredictorIntegrator, CECMIntegrator, CF4Integrator, CELIIntegrator, + EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator) from tests import dummy_operator +INTEGRATORS = [ + PredictorIntegrator, + CECMIntegrator, + CF4Integrator, + CELIIntegrator, + EPCRK4Integrator, + LEQIIntegrator, + SICELIIntegrator, + SILEQIIntegrator +] + + def test_results_save(run_in_tmpdir): """Test data save module""" @@ -159,3 +173,56 @@ def test_integrator(run_in_tmpdir, scheme): dep_time = res.get_depletion_time() assert dep_time.shape == (2, ) assert all(dep_time > 0) + + +@pytest.mark.parametrize("integrator", INTEGRATORS) +def test_timesteps(integrator): + # Crate fake operator + op = MagicMock() + op.prev_res = None + op.chain = None + + # Set heavy metal mass and power randomly + op.heavy_metal = uniform(0, 10000) + power = uniform(0, 1e6) + + # Reference timesteps in seconds + day = 86400.0 + ref_timesteps = [1*day, 2*day, 5*day, 10*day] + + # Case 1, timesteps in second + timesteps = ref_timesteps + x = integrator(op, timesteps, power, timestep_units='s') + assert np.allclose(x.timesteps, ref_timesteps) + + # Case 2, timesteps in days + timesteps = [t / day for t in ref_timesteps] + x = integrator(op, timesteps, power, timestep_units='d') + assert np.allclose(x.timesteps, ref_timesteps) + + # Case 3, timesteps in years + year = 365.25*day + timesteps = [t / year for t in ref_timesteps] + x = integrator(op, timesteps, power, timestep_units='a') + assert np.allclose(x.timesteps, ref_timesteps) + + # Case 4, timesteps in MWd/kg + kilograms = op.heavy_metal / 1000.0 + days = [t/day for t in ref_timesteps] + megawatts = power / 1000000.0 + burnup = [t * megawatts / kilograms for t in days] + x = integrator(op, burnup, power, timestep_units='MWd/kg') + assert np.allclose(x.timesteps, ref_timesteps) + + # Case 5, mixed units + burnup_per_day = (1e-6*power) / kilograms + timesteps = [(burnup_per_day, 'MWd/kg'), (2*day, 's'), (5, 'd'), + (10*burnup_per_day, 'MWd/kg')] + x = integrator(op, timesteps, power) + assert np.allclose(x.timesteps, ref_timesteps) + + # Bad units should raise an exception + with pytest.raises(ValueError, match="unit"): + integrator(op, ref_timesteps, power, timestep_units='🐨') + with pytest.raises(ValueError, match="unit"): + integrator(op, [(800.0, 'gorillas')], power) From 85264aa38c60d6846de5ca7005373c7cb3b9a47c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 11 Feb 2020 07:08:53 -0600 Subject: [PATCH 03/59] Support minutes and hours for depletion timesteps rather than years --- openmc/deplete/abc.py | 29 +++++----- openmc/deplete/integrators.py | 64 ++++++++++----------- tests/unit_tests/test_deplete_integrator.py | 26 +++++---- 3 files changed, 64 insertions(+), 55 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index a9640852cc..d75833a6f1 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -31,8 +31,9 @@ __all__ = [ "Integrator", "SIIntegrator", "DepSystemSolver"] +_SECONDS_PER_MINUTE = 60 +_SECONDS_PER_HOUR = 60*60 _SECONDS_PER_DAY = 24*60*60 -_SECONDS_PER_YEAR = 365.25*24*60*60 OperatorResult = namedtuple('OperatorResult', ['k', 'rates']) OperatorResult.__doc__ = """\ @@ -617,11 +618,11 @@ class Integrator(ABC): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). Attributes ---------- @@ -678,19 +679,21 @@ class Integrator(ABC): # Determine number of seconds for each timestep seconds = [] for time, unit, watts in zip(times, units, power): - if unit == 's': + if unit in ('s', 'sec'): seconds.append(time) + elif unit in ('min', 'minute'): + seconds.append(time*_SECONDS_PER_MINUTE) + elif unit in ('h', 'hr', 'hour'): + seconds.append(time*_SECONDS_PER_HOUR) elif unit in ('d', 'day'): seconds.append(time*_SECONDS_PER_DAY) - elif unit in ('a', 'yr', 'year'): - seconds.append(time*_SECONDS_PER_YEAR) elif unit.lower() == 'mwd/kg': watt_days_per_kg = 1e6*time kilograms = 1e-3*mass days = watt_days_per_kg * kilograms / watts seconds.append(days*_SECONDS_PER_DAY) else: - raise ValueError("Invalid timestep unit: {}".format(unit)) + raise ValueError("Invalid timestep unit '{}'".format(unit)) self.timesteps = asarray(seconds) self.power = asarray(power) @@ -824,11 +827,11 @@ class SIIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index 84f838cb7a..da8509d591 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -48,11 +48,11 @@ class PredictorIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). Attributes ---------- @@ -137,11 +137,11 @@ class CECMIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). Attributes ---------- @@ -234,11 +234,11 @@ class CF4Integrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). Attributes ---------- @@ -348,11 +348,11 @@ class CELIIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). Attributes ---------- @@ -449,11 +449,11 @@ class EPCRK4Integrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). Attributes ---------- @@ -570,11 +570,11 @@ class LEQIIntegrator(Integrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). Attributes ---------- @@ -688,11 +688,11 @@ class SICELIIntegrator(SIIntegrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 @@ -796,11 +796,11 @@ class SILEQIIntegrator(SIIntegrator): Power density of the reactor in [W/gHM]. It is multiplied by initial heavy metal inventory to get total power if ``power`` is not speficied. - timestep_units : {'s', 'd', 'a', 'MWd/kg'} + timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'} Units for values specified in the `timesteps` argument. 's' means - seconds, 'd' means days, 'a' means years, and 'MWd/kg' indicates that - the values are given in burnup (MW-d of energy deposited per kilogram - of heavy metal). + seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates + that the values are given in burnup (MW-d of energy deposited per + kilogram of heavy metal). n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 9d5766fd5c..110894cdad 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -190,23 +190,29 @@ def test_timesteps(integrator): day = 86400.0 ref_timesteps = [1*day, 2*day, 5*day, 10*day] - # Case 1, timesteps in second + # Case 1, timesteps in seconds timesteps = ref_timesteps x = integrator(op, timesteps, power, timestep_units='s') assert np.allclose(x.timesteps, ref_timesteps) - # Case 2, timesteps in days + # Case 2, timesteps in minutes + minute = 60 + timesteps = [t / minute for t in ref_timesteps] + x = integrator(op, timesteps, power, timestep_units='min') + assert np.allclose(x.timesteps, ref_timesteps) + + # Case 3, timesteps in hours + hour = 60*60 + timesteps = [t / hour for t in ref_timesteps] + x = integrator(op, timesteps, power, timestep_units='h') + assert np.allclose(x.timesteps, ref_timesteps) + + # Case 4, timesteps in days timesteps = [t / day for t in ref_timesteps] x = integrator(op, timesteps, power, timestep_units='d') assert np.allclose(x.timesteps, ref_timesteps) - # Case 3, timesteps in years - year = 365.25*day - timesteps = [t / year for t in ref_timesteps] - x = integrator(op, timesteps, power, timestep_units='a') - assert np.allclose(x.timesteps, ref_timesteps) - - # Case 4, timesteps in MWd/kg + # Case 5, timesteps in MWd/kg kilograms = op.heavy_metal / 1000.0 days = [t/day for t in ref_timesteps] megawatts = power / 1000000.0 @@ -214,7 +220,7 @@ def test_timesteps(integrator): x = integrator(op, burnup, power, timestep_units='MWd/kg') assert np.allclose(x.timesteps, ref_timesteps) - # Case 5, mixed units + # Case 6, mixed units burnup_per_day = (1e-6*power) / kilograms timesteps = [(burnup_per_day, 'MWd/kg'), (2*day, 's'), (5, 'd'), (10*burnup_per_day, 'MWd/kg')] From 516ac9e1d372ff2030939dc87fd0699d81bb3c9c Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Tue, 11 Feb 2020 13:37:00 +0000 Subject: [PATCH 04/59] Modified the MAX_LOST_PARTICLES const in particle.h to be a variable stored in settings that can be modified by the c/python/xml layers. --- include/openmc/particle.h | 3 --- include/openmc/settings.h | 9 +++++---- openmc/lib/settings.py | 1 + openmc/settings.py | 26 ++++++++++++++++++++++++++ openmc/statepoint.py | 6 ++++++ src/particle.cpp | 6 +++++- src/settings.cpp | 12 ++++++++++-- 7 files changed, 53 insertions(+), 10 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 8143ee5a30..ee59d73014 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -34,9 +34,6 @@ namespace openmc { // use to store the bins for delayed group tallies. constexpr int MAX_DELAYED_GROUPS {8}; -// Maximum number of lost particles -constexpr int MAX_LOST_PARTICLES {10}; - // Maximum number of lost particles, relative to the total number of particles constexpr double REL_MAX_LOST_PARTICLES {1.0e-6}; diff --git a/include/openmc/settings.h b/include/openmc/settings.h index cbc0c488b9..1dcdea3cc7 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -63,10 +63,11 @@ extern std::string path_source; extern std::string path_sourcepoint; //!< path to a source file extern "C" std::string path_statepoint; //!< path to a statepoint file -extern "C" int32_t n_batches; //!< number of (inactive+active) batches -extern "C" int32_t n_inactive; //!< number of inactive batches -extern "C" int32_t gen_per_batch; //!< number of generations per batch -extern "C" int64_t n_particles; //!< number of particles per generation +extern "C" int32_t n_batches; //!< number of (inactive+active) batches +extern "C" int32_t n_inactive; //!< number of inactive batches +extern "C" int32_t n_max_lost_particles; //!< maximum number of lost particles +extern "C" int32_t gen_per_batch; //!< number of generations per batch +extern "C" int64_t n_particles; //!< number of particles per generation extern int64_t max_particles_in_flight; //!< Max num. event-based particles in flight diff --git a/openmc/lib/settings.py b/openmc/lib/settings.py index 275c73a934..a25ccb962f 100644 --- a/openmc/lib/settings.py +++ b/openmc/lib/settings.py @@ -22,6 +22,7 @@ class _Settings(object): entropy_on = _DLLGlobal(c_bool, 'entropy_on') generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch') inactive = _DLLGlobal(c_int32, 'n_inactive') + max_lost_particles = _DLLGlobal(c_int32, 'n_max_lost_particles') particles = _DLLGlobal(c_int64, 'n_particles') restart_run = _DLLGlobal(c_bool, 'restart_run') run_CE = _DLLGlobal(c_bool, 'run_CE') diff --git a/openmc/settings.py b/openmc/settings.py index eb2f926f3a..79e7719a3d 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -58,6 +58,8 @@ class Settings(object): history-based parallelism. generations_per_batch : int Number of generations per batch + max_lost_particles : int + Maximum number of lost particles inactive : int Number of inactive batches keff_trigger : dict @@ -176,6 +178,7 @@ class Settings(object): self._batches = None self._generations_per_batch = None self._inactive = None + self._max_lost_particles = None self._particles = None self._keff_trigger = None @@ -254,6 +257,10 @@ class Settings(object): def inactive(self): return self._inactive + @property + def max_lost_particles(self): + return self._max_lost_particles + @property def particles(self): return self._particles @@ -417,6 +424,12 @@ class Settings(object): cv.check_greater_than('inactive batches', inactive, 0, True) self._inactive = inactive + @max_lost_particles.setter + def max_lost_particles(self, max_lost_particles): + cv.check_type('max_lost_particles', max_lost_particles, Integral) + cv.check_greater_than('max_lost_particles', max_lost_particles, 0) + self._max_lost_particles = max_lost_particles + @particles.setter def particles(self, particles): cv.check_type('particles', particles, Integral) @@ -763,6 +776,11 @@ class Settings(object): element = ET.SubElement(root, "inactive") element.text = str(self._inactive) + def _create_max_lost_particles_subelement(self, root): + if self._max_lost_particles is not None: + element = ET.SubElement(root, "max_lost_particles") + element.text = str(self._max_lost_particles) + def _create_particles_subelement(self, root): if self._particles is not None: element = ET.SubElement(root, "particles") @@ -1009,6 +1027,7 @@ class Settings(object): self._particles_from_xml_element(elem) self._batches_from_xml_element(elem) self._inactive_from_xml_element(elem) + self._max_lost_particles_from_xml_element(elem) self._generations_per_batch_from_xml_element(elem) def _run_mode_from_xml_element(self, root): @@ -1031,6 +1050,11 @@ class Settings(object): if text is not None: self.inactive = int(text) + def _max_lost_particles_from_xml_element(self, root): + text = get_text(root, 'max_lost_particles') + if text is not None: + self.max_lost_particles = int(text) + def _generations_per_batch_from_xml_element(self, root): text = get_text(root, 'generations_per_batch') if text is not None: @@ -1270,6 +1294,7 @@ class Settings(object): self._create_particles_subelement(root_element) self._create_batches_subelement(root_element) self._create_inactive_subelement(root_element) + self._create_max_lost_particles_subelement(root_element) self._create_generations_per_batch_subelement(root_element) self._create_keff_trigger_subelement(root_element) self._create_source_subelement(root_element) @@ -1340,6 +1365,7 @@ class Settings(object): settings._particles_from_xml_element(root) settings._batches_from_xml_element(root) settings._inactive_from_xml_element(root) + settings._max_lost_particles_from_xml_element(root) settings._generations_per_batch_from_xml_element(root) settings._keff_trigger_from_xml_element(root) settings._source_from_xml_element(root) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 5802a74cae..051195ec96 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -77,6 +77,8 @@ class StatePoint(object): Number of batches n_inactive : int Number of inactive batches + n_max_lost_particles : int + Number of max lost particles n_particles : int Number of particles per generation n_realizations : int @@ -312,6 +314,10 @@ class StatePoint(object): else: return None + @property + def n_max_lost_particles(self): + return self._f['n_max_lost_particles'][()] + @property def n_particles(self): return self._f['n_particles'][()] diff --git a/src/particle.cpp b/src/particle.cpp index bd9a611392..8e66222048 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -164,6 +164,8 @@ Particle::event_calculate_xs() if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell; } + std::cout << "settings::n_max_lost_particles " << settings::n_max_lost_particles << std::endl; + // Write particle track. if (write_track_) write_particle_track(*this); @@ -628,9 +630,11 @@ Particle::mark_as_lost(const char* message) auto n = simulation::current_batch * settings::gen_per_batch * simulation::work_per_rank; + std::cout << simulation::n_lost_particles << "HELLLOOOOOOOOOOOO MS DOUBTFIRE" << std::endl; + // Abort the simulation if the maximum number of lost particles has been // reached - if (simulation::n_lost_particles >= MAX_LOST_PARTICLES && + if (simulation::n_lost_particles >= settings::n_max_lost_particles && simulation::n_lost_particles >= REL_MAX_LOST_PARTICLES*n) { fatal_error("Maximum number of lost particles has been reached."); } diff --git a/src/settings.cpp b/src/settings.cpp index 32a7958d8d..41a55501c2 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -79,6 +79,7 @@ std::string path_statepoint; int32_t n_batches; int32_t n_inactive {0}; +int32_t n_max_lost_particles {10}; int32_t gen_per_batch {1}; int64_t n_particles {-1}; @@ -143,6 +144,11 @@ void get_run_parameters(pugi::xml_node node_base) } if (!trigger_on) n_max_batches = n_batches; + // Get max number of lost particles + if (check_for_node(node_base, "max_lost_particles")) { + n_max_lost_particles = std::stoi(get_node_value(node_base, "max_lost_particles")); + } + // Get number of inactive batches if (run_mode == RunMode::EIGENVALUE) { if (check_for_node(node_base, "inactive")) { @@ -343,14 +349,16 @@ void read_settings_xml() // Read run parameters get_run_parameters(node_mode); - // Check number of active batches, inactive batches, and particles + // Check number of active batches, inactive batches, max lost particles and particles if (n_batches <= n_inactive) { fatal_error("Number of active batches must be greater than zero."); } else if (n_inactive < 0) { fatal_error("Number of inactive batches must be non-negative."); } else if (n_particles <= 0) { fatal_error("Number of particles must be greater than zero."); - } + } else if (n_max_lost_particles <= 0) { + fatal_error("Number of max lost particles must be greater than zero."); + } } // Copy random number seed if specified From 9f64dd1a880a150682ae6140da1a9546b1982a35 Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Tue, 11 Feb 2020 13:37:26 +0000 Subject: [PATCH 05/59] Modified the REL_MAX_LOST_PARTICLES const in particle.h to be a variable stored in settings that can be modified by the c/python/xml layers. --- include/openmc/particle.h | 3 --- include/openmc/settings.h | 11 ++++++----- openmc/lib/settings.py | 1 + openmc/settings.py | 31 +++++++++++++++++++++++++++++-- openmc/statepoint.py | 8 +++++++- src/particle.cpp | 6 +----- src/settings.cpp | 10 +++++++++- 7 files changed, 53 insertions(+), 17 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index ee59d73014..2917b45e4a 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -34,9 +34,6 @@ namespace openmc { // use to store the bins for delayed group tallies. constexpr int MAX_DELAYED_GROUPS {8}; -// Maximum number of lost particles, relative to the total number of particles -constexpr double REL_MAX_LOST_PARTICLES {1.0e-6}; - constexpr double CACHE_INVALID {-1.0}; //============================================================================== diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 1dcdea3cc7..3407eadf4b 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -63,11 +63,12 @@ extern std::string path_source; extern std::string path_sourcepoint; //!< path to a source file extern "C" std::string path_statepoint; //!< path to a statepoint file -extern "C" int32_t n_batches; //!< number of (inactive+active) batches -extern "C" int32_t n_inactive; //!< number of inactive batches -extern "C" int32_t n_max_lost_particles; //!< maximum number of lost particles -extern "C" int32_t gen_per_batch; //!< number of generations per batch -extern "C" int64_t n_particles; //!< number of particles per generation +extern "C" int32_t n_batches; //!< number of (inactive+active) batches +extern "C" int32_t n_inactive; //!< number of inactive batches +extern "C" int32_t n_max_lost_particles; //!< maximum number of lost particles +extern double relative_max_lost_particles; //!< maximum number of lost particles, relative to the total number of particles +extern "C" int32_t gen_per_batch; //!< number of generations per batch +extern "C" int64_t n_particles; //!< number of particles per generation extern int64_t max_particles_in_flight; //!< Max num. event-based particles in flight diff --git a/openmc/lib/settings.py b/openmc/lib/settings.py index a25ccb962f..8345e00881 100644 --- a/openmc/lib/settings.py +++ b/openmc/lib/settings.py @@ -23,6 +23,7 @@ class _Settings(object): generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch') inactive = _DLLGlobal(c_int32, 'n_inactive') max_lost_particles = _DLLGlobal(c_int32, 'n_max_lost_particles') + rel_max_lost_particles = _DLLGlobal(c_double, 'relative_max_lost_particles') particles = _DLLGlobal(c_int64, 'n_particles') restart_run = _DLLGlobal(c_bool, 'restart_run') run_CE = _DLLGlobal(c_bool, 'run_CE') diff --git a/openmc/settings.py b/openmc/settings.py index 79e7719a3d..794c97b618 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -60,6 +60,8 @@ class Settings(object): Number of generations per batch max_lost_particles : int Maximum number of lost particles + rel_max_lost_particles : int + Maximum number of lost particles, relative to the total number of particles inactive : int Number of inactive batches keff_trigger : dict @@ -179,6 +181,7 @@ class Settings(object): self._generations_per_batch = None self._inactive = None self._max_lost_particles = None + self._rel_max_lost_particles = None self._particles = None self._keff_trigger = None @@ -261,6 +264,10 @@ class Settings(object): def max_lost_particles(self): return self._max_lost_particles + @property + def rel_max_lost_particles(self): + return self._rel_max_lost_particles + @property def particles(self): return self._particles @@ -430,6 +437,13 @@ class Settings(object): cv.check_greater_than('max_lost_particles', max_lost_particles, 0) self._max_lost_particles = max_lost_particles + @rel_max_lost_particles.setter + def rel_max_lost_particles(self, rel_max_lost_particles): + cv.check_type('rel_max_lost_particles', rel_max_lost_particles, Real) + cv.check_greater_than('rel_max_lost_particles', rel_max_lost_particles, 0) + cv.check_less_than('rel_max_lost_particles', rel_max_lost_particles, 1) + self._rel_max_lost_particles = rel_max_lost_particles + @particles.setter def particles(self, particles): cv.check_type('particles', particles, Integral) @@ -779,7 +793,12 @@ class Settings(object): def _create_max_lost_particles_subelement(self, root): if self._max_lost_particles is not None: element = ET.SubElement(root, "max_lost_particles") - element.text = str(self._max_lost_particles) + element.text = str(self._max_lost_particles) + + def _create_rel_max_lost_particles_subelement(self, root): + if self._rel_max_lost_particles is not None: + element = ET.SubElement(root, "rel_max_lost_particles") + element.text = str(self._rel_max_lost_particles) def _create_particles_subelement(self, root): if self._particles is not None: @@ -1028,6 +1047,7 @@ class Settings(object): self._batches_from_xml_element(elem) self._inactive_from_xml_element(elem) self._max_lost_particles_from_xml_element(elem) + self._rel_max_lost_particles_from_xml_element(elem) self._generations_per_batch_from_xml_element(elem) def _run_mode_from_xml_element(self, root): @@ -1053,7 +1073,12 @@ class Settings(object): def _max_lost_particles_from_xml_element(self, root): text = get_text(root, 'max_lost_particles') if text is not None: - self.max_lost_particles = int(text) + self.max_lost_particles = int(text) + + def _rel_max_lost_particles_from_xml_element(self, root): + text = get_text(root, 'rel_max_lost_particles') + if text is not None: + self.rel_max_lost_particles = float(text) def _generations_per_batch_from_xml_element(self, root): text = get_text(root, 'generations_per_batch') @@ -1295,6 +1320,7 @@ class Settings(object): self._create_batches_subelement(root_element) self._create_inactive_subelement(root_element) self._create_max_lost_particles_subelement(root_element) + self._create_rel_max_lost_particles_subelement(root_element) self._create_generations_per_batch_subelement(root_element) self._create_keff_trigger_subelement(root_element) self._create_source_subelement(root_element) @@ -1366,6 +1392,7 @@ class Settings(object): settings._batches_from_xml_element(root) settings._inactive_from_xml_element(root) settings._max_lost_particles_from_xml_element(root) + settings._rel_max_lost_particles_from_xml_element(root) settings._generations_per_batch_from_xml_element(root) settings._keff_trigger_from_xml_element(root) settings._source_from_xml_element(root) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 051195ec96..0d78bf68dc 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -79,6 +79,8 @@ class StatePoint(object): Number of inactive batches n_max_lost_particles : int Number of max lost particles + relative_max_lost_particles : float + Number of max lost particles, relative to the total number of particles n_particles : int Number of particles per generation n_realizations : int @@ -316,7 +318,11 @@ class StatePoint(object): @property def n_max_lost_particles(self): - return self._f['n_max_lost_particles'][()] + return self._f['n_max_lost_particles'][()] + + @property + def relative_max_lost_particles(self): + return self._f['relative_max_lost_particles'][()] @property def n_particles(self): diff --git a/src/particle.cpp b/src/particle.cpp index 8e66222048..79e6625771 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -164,8 +164,6 @@ Particle::event_calculate_xs() if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell; } - std::cout << "settings::n_max_lost_particles " << settings::n_max_lost_particles << std::endl; - // Write particle track. if (write_track_) write_particle_track(*this); @@ -630,12 +628,10 @@ Particle::mark_as_lost(const char* message) auto n = simulation::current_batch * settings::gen_per_batch * simulation::work_per_rank; - std::cout << simulation::n_lost_particles << "HELLLOOOOOOOOOOOO MS DOUBTFIRE" << std::endl; - // Abort the simulation if the maximum number of lost particles has been // reached if (simulation::n_lost_particles >= settings::n_max_lost_particles && - simulation::n_lost_particles >= REL_MAX_LOST_PARTICLES*n) { + simulation::n_lost_particles >= settings::relative_max_lost_particles*n) { fatal_error("Maximum number of lost particles has been reached."); } } diff --git a/src/settings.cpp b/src/settings.cpp index 41a55501c2..d72e5f17f7 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -80,6 +80,7 @@ std::string path_statepoint; int32_t n_batches; int32_t n_inactive {0}; int32_t n_max_lost_particles {10}; +double relative_max_lost_particles {1.0e-6}; int32_t gen_per_batch {1}; int64_t n_particles {-1}; @@ -149,6 +150,11 @@ void get_run_parameters(pugi::xml_node node_base) n_max_lost_particles = std::stoi(get_node_value(node_base, "max_lost_particles")); } + // Get relative number of lost particles + if (check_for_node(node_base, "rel_max_lost_particles")) { + relative_max_lost_particles = std::stod(get_node_value(node_base, "rel_max_lost_particles")); + } + // Get number of inactive batches if (run_mode == RunMode::EIGENVALUE) { if (check_for_node(node_base, "inactive")) { @@ -358,7 +364,9 @@ void read_settings_xml() fatal_error("Number of particles must be greater than zero."); } else if (n_max_lost_particles <= 0) { fatal_error("Number of max lost particles must be greater than zero."); - } + } else if (relative_max_lost_particles <= 0.0 || relative_max_lost_particles >= 1.0) { + fatal_error("Relative max lost particles must be between zero and one."); + } } // Copy random number seed if specified From 6a40a535e482aa0ae31949b0836c22f3d6445f48 Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Tue, 11 Feb 2020 13:37:32 +0000 Subject: [PATCH 06/59] Modified the test_settings unit test for new max_lost_particles and rel_max_lost_particles to check they can be assigned/read by the python interface --- tests/unit_tests/test_settings.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index 9c33f6a2d4..c3ff728798 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -9,6 +9,8 @@ def test_export_to_xml(run_in_tmpdir): s.generations_per_batch = 10 s.inactive = 100 s.particles = 1000000 + s.max_lost_particles = 5 + s.rel_max_lost_particles = 1e-4 s.keff_trigger = {'type': 'std_dev', 'threshold': 0.001} s.energy_mode = 'continuous-energy' s.max_order = 5 @@ -62,6 +64,8 @@ def test_export_to_xml(run_in_tmpdir): assert s.generations_per_batch == 10 assert s.inactive == 100 assert s.particles == 1000000 + assert s.max_lost_particles == 5 + assert s.rel_max_lost_particles == 1e-4 assert s.keff_trigger == {'type': 'std_dev', 'threshold': 0.001} assert s.energy_mode == 'continuous-energy' assert s.max_order == 5 From 9224ca318675fcbb06c8eb97d9569466874dda21 Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Tue, 11 Feb 2020 13:37:37 +0000 Subject: [PATCH 07/59] Modified an existing regression test to include the new max_lost_particles and relative version. This is mainly a test that the python/xml input/output routines work through the comparison the regression test harness does in this test, neither of these trigger the end of the simulation. --- tests/regression_tests/photon_source/inputs_true.dat | 2 ++ tests/regression_tests/photon_source/test.py | 2 ++ 2 files changed, 4 insertions(+) diff --git a/tests/regression_tests/photon_source/inputs_true.dat b/tests/regression_tests/photon_source/inputs_true.dat index 89f4de0e0c..83c909ac69 100644 --- a/tests/regression_tests/photon_source/inputs_true.dat +++ b/tests/regression_tests/photon_source/inputs_true.dat @@ -18,6 +18,8 @@ fixed source 10000 1 + 5 + 0.1 0 0 0 diff --git a/tests/regression_tests/photon_source/test.py b/tests/regression_tests/photon_source/test.py index 94b0280031..09a0b4a2d0 100644 --- a/tests/regression_tests/photon_source/test.py +++ b/tests/regression_tests/photon_source/test.py @@ -31,6 +31,8 @@ class SourceTestHarness(PyAPITestHarness): settings = openmc.Settings() settings.particles = 10000 settings.batches = 1 + settings.max_lost_particles = 5 + settings.rel_max_lost_particles = 0.1 settings.photon_transport = True settings.electron_treatment = 'ttb' settings.cutoff = {'energy_photon' : 1000.0} From b6a03949c82aee90824b9daf60d3d74a587b5b19 Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Tue, 11 Feb 2020 13:37:41 +0000 Subject: [PATCH 08/59] Modified the StatePoint output to include the new max_lost_particles and relative version. Further modified an existing regression test to check the StatePoint write-out of these variables --- src/state_point.cpp | 4 ++++ tests/regression_tests/photon_source/results_true.dat | 2 ++ tests/regression_tests/photon_source/test.py | 2 ++ 3 files changed, 8 insertions(+) diff --git a/src/state_point.cpp b/src/state_point.cpp index 60fd3826de..a9213bcb4e 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -93,6 +93,8 @@ openmc_statepoint_write(const char* filename, bool* write_source) write_attribute(file_id, "photon_transport", settings::photon_transport); write_dataset(file_id, "n_particles", settings::n_particles); write_dataset(file_id, "n_batches", settings::n_batches); + write_dataset(file_id, "n_max_lost_particles", settings::n_max_lost_particles); + write_dataset(file_id, "relative_max_lost_particles", settings::relative_max_lost_particles); // Write out current batch number write_dataset(file_id, "current_batch", simulation::current_batch); @@ -373,6 +375,8 @@ void load_state_point() read_dataset(file_id, "n_particles", settings::n_particles); int temp; read_dataset(file_id, "n_batches", temp); + read_dataset(file_id, "n_max_lost_particles", settings::n_max_lost_particles); + read_dataset(file_id, "relative_max_lost_particles", settings::relative_max_lost_particles); // Take maximum of statepoint n_batches and input n_batches settings::n_batches = std::max(settings::n_batches, temp); diff --git a/tests/regression_tests/photon_source/results_true.dat b/tests/regression_tests/photon_source/results_true.dat index 1448ad6db8..ab6ac9ff83 100644 --- a/tests/regression_tests/photon_source/results_true.dat +++ b/tests/regression_tests/photon_source/results_true.dat @@ -1,3 +1,5 @@ tally 1: sum = 2.275713E+02 sum_sq = 5.178870E+04 +max_lost = 5.000000E+00 +rel_max_lost = 1.000000E-01 diff --git a/tests/regression_tests/photon_source/test.py b/tests/regression_tests/photon_source/test.py index 09a0b4a2d0..c5c42d680e 100644 --- a/tests/regression_tests/photon_source/test.py +++ b/tests/regression_tests/photon_source/test.py @@ -54,6 +54,8 @@ class SourceTestHarness(PyAPITestHarness): outstr += 'tally {}:\n'.format(t.id) outstr += 'sum = {:12.6E}\n'.format(t.sum[0, 0, 0]) outstr += 'sum_sq = {:12.6E}\n'.format(t.sum_sq[0, 0, 0]) + outstr += 'max_lost = {:12.6E}\n'.format(sp.n_max_lost_particles) + outstr += 'rel_max_lost = {:12.6E}\n'.format(sp.relative_max_lost_particles) return outstr From 70e914d2458f6c35ef46f25b93f221f056dac813 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 11 Feb 2020 10:12:34 -0600 Subject: [PATCH 09/59] Fix failing tests --- openmc/deplete/abc.py | 7 +++---- tests/unit_tests/test_deplete_integrator.py | 4 ++-- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index d75833a6f1..912c0b1fb4 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -652,14 +652,13 @@ class Integrator(ABC): self.chain = operator.chain # Determine power and normalize units to W - mass = operator.heavy_metal if power is None: if power_density is None: raise ValueError("Either power or power density must be set") if not isinstance(power_density, Iterable): - power = power_density * mass + power = power_density * operator.heavy_metal else: - power = [p*mass for p in power_density] + power = [p*operator.heavy_metal for p in power_density] if not isinstance(power, Iterable): # Ensure that power is single value if that is the case power = [power] * len(timesteps) @@ -689,7 +688,7 @@ class Integrator(ABC): seconds.append(time*_SECONDS_PER_DAY) elif unit.lower() == 'mwd/kg': watt_days_per_kg = 1e6*time - kilograms = 1e-3*mass + kilograms = 1e-3*operator.heavy_metal days = watt_days_per_kg * kilograms / watts seconds.append(days*_SECONDS_PER_DAY) else: diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 110894cdad..6c09b3feac 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -119,14 +119,14 @@ def test_results_save(run_in_tmpdir): np.testing.assert_array_equal(res[1].time, t2) -@pytest.mark.parametrize("timesteps", (1, [1])) -def test_bad_integrator_inputs(timesteps): +def test_bad_integrator_inputs(): """Test failure modes for Integrator inputs""" op = MagicMock() op.prev_res = None op.chain = None op.heavy_metal = 1.0 + timesteps = [1] # No power nor power density given with pytest.raises(ValueError, match="Either power or power density"): From 39a60c361c9111575d9ded988c8459dbc4b82342 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 12 Feb 2020 07:19:31 -0600 Subject: [PATCH 10/59] Clarify that heavy metal is initial Co-Authored-By: Andrew Johnson --- openmc/deplete/abc.py | 4 ++-- openmc/deplete/integrators.py | 16 ++++++++-------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 912c0b1fb4..658e93da5a 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -622,7 +622,7 @@ class Integrator(ABC): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). Attributes ---------- @@ -830,7 +830,7 @@ class SIIntegrator(Integrator): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 diff --git a/openmc/deplete/integrators.py b/openmc/deplete/integrators.py index da8509d591..38a7c33b1f 100644 --- a/openmc/deplete/integrators.py +++ b/openmc/deplete/integrators.py @@ -52,7 +52,7 @@ class PredictorIntegrator(Integrator): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). Attributes ---------- @@ -141,7 +141,7 @@ class CECMIntegrator(Integrator): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). Attributes ---------- @@ -238,7 +238,7 @@ class CF4Integrator(Integrator): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). Attributes ---------- @@ -352,7 +352,7 @@ class CELIIntegrator(Integrator): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). Attributes ---------- @@ -453,7 +453,7 @@ class EPCRK4Integrator(Integrator): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). Attributes ---------- @@ -574,7 +574,7 @@ class LEQIIntegrator(Integrator): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). Attributes ---------- @@ -692,7 +692,7 @@ class SICELIIntegrator(SIIntegrator): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 @@ -800,7 +800,7 @@ class SILEQIIntegrator(SIIntegrator): Units for values specified in the `timesteps` argument. 's' means seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates that the values are given in burnup (MW-d of energy deposited per - kilogram of heavy metal). + kilogram of initial heavy metal). n_steps : int, optional Number of stochastic iterations per depletion interval. Must be greater than zero. Default : 10 From 3682f7adcbf6420198a39e332207818685a9524f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 12 Feb 2020 07:36:29 -0600 Subject: [PATCH 11/59] Add type/value checks on timestep, power, and units --- openmc/deplete/abc.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 658e93da5a..95bd3326b8 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -665,7 +665,7 @@ class Integrator(ABC): if len(power) != len(timesteps): raise ValueError( - "Number of time steps != number of powers. {} vs {}".format( + "Number of time steps ({}) != number of powers ({})".format( len(timesteps), len(power))) # Get list of times / units @@ -678,6 +678,13 @@ class Integrator(ABC): # Determine number of seconds for each timestep seconds = [] for time, unit, watts in zip(times, units, power): + # Make sure values passed make sense + check_type('timestep', time, Real) + check_greater_than('timestep', time, 0.0, True) + check_type('timestep units', unit, str) + check_type('power', watts, Real) + check_greater_than('power', watts, 0.0, True) + if unit in ('s', 'sec'): seconds.append(time) elif unit in ('min', 'minute'): From 4966ff7c5dff243853881c75125994dda6d3a887 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 12 Feb 2020 12:21:25 -0600 Subject: [PATCH 12/59] Ensure timesteps are positive Co-Authored-By: Andrew Johnson --- openmc/deplete/abc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 95bd3326b8..31479f9952 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -680,7 +680,7 @@ class Integrator(ABC): for time, unit, watts in zip(times, units, power): # Make sure values passed make sense check_type('timestep', time, Real) - check_greater_than('timestep', time, 0.0, True) + check_greater_than('timestep', time, 0.0, False) check_type('timestep units', unit, str) check_type('power', watts, Real) check_greater_than('power', watts, 0.0, True) From 4200b0dfbf11cf555d3f791285508bf0a18c87fd Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 12 Feb 2020 16:05:58 -0600 Subject: [PATCH 13/59] Update depletion-related documentation --- docs/source/usersguide/depletion.rst | 17 ++- .../pincell_depletion/restart_depletion.py | 50 +++---- .../python/pincell_depletion/run_depletion.py | 122 ++++++------------ 3 files changed, 67 insertions(+), 122 deletions(-) diff --git a/docs/source/usersguide/depletion.rst b/docs/source/usersguide/depletion.rst index 4a4c7b1f88..f6feb6d245 100644 --- a/docs/source/usersguide/depletion.rst +++ b/docs/source/usersguide/depletion.rst @@ -33,12 +33,11 @@ material compositions over time. Each method appears as a different class. For example, :class:`openmc.deplete.CECMIntegrator` runs a depletion calculation using the CE/CM algorithm (deplete over a timestep using the middle-of-step reaction rates). An instance of :class:`openmc.deplete.Operator` is passed to -one of these functions along with the power level and timesteps:: +one of these functions along with the timesteps and power level:: - power = 1200.0e6 - days = 24*60*60 - timesteps = [10.0*days, 10.0*days, 10.0*days] - openmc.deplete.CECMIntegrator(op, power, timesteps).integrate() + power = 1200.0e6 # watts + timesteps = [10.0, 10.0, 10.0] # days + openmc.deplete.CECMIntegrator(op, timesteps, power, timestep_units='d').integrate() The coupled transport-depletion problem is executed, and once it is done a ``depletion_results.h5`` file is written. The results can be analyzed using the @@ -67,7 +66,7 @@ the energy deposited during a transport calculation will be lower than expected. This causes the reaction rates to be over-adjusted to hit the user-specific power, or power density, leading to an over-depletion of burnable materials. -There are some remedies. First, the fission Q values can be directly set in a +There are some remedies. First, the fission Q values can be directly set in a variety of ways. This requires knowing what the total fission energy release should be, including indirect components. Some examples are provided below:: @@ -99,11 +98,11 @@ Local Spectra and Repeated Materials ------------------------------------ It is not uncommon to explicitly create a single burnable material across many locations. -From a pure transport perspective, there is nothing wrong with creating a single +From a pure transport perspective, there is nothing wrong with creating a single 3.5 wt.% enriched fuel ``fuel_3``, and placing that fuel in every fuel pin in an assembly or even full core problem. This certainly expedites the model making process, but can pose -issues with depletion. -Under this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material using +issues with depletion. +Under this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material using a single set of reaction rates, and produce a single new composition for the next time step. This can be problematic if the same ``fuel_3`` is used in very different regions of the problem. diff --git a/examples/python/pincell_depletion/restart_depletion.py b/examples/python/pincell_depletion/restart_depletion.py index 013a2469e0..f0387e066f 100644 --- a/examples/python/pincell_depletion/restart_depletion.py +++ b/examples/python/pincell_depletion/restart_depletion.py @@ -1,25 +1,7 @@ import openmc import openmc.deplete -import numpy as np import matplotlib.pyplot as plt -############################################################################### -# Simulation Input File Parameters -############################################################################### - -# OpenMC simulation parameters -batches = 100 -inactive = 10 -particles = 1000 - -# Depletion simulation parameters -time_step = 1*24*60*60 # s -final_time = 5*24*60*60 # s -time_steps = np.full(final_time // time_step, time_step) - -chain_file = './chain_simple.xml' -power = 174 # W/cm, for 2D simulations only (use W for 3D) - ############################################################################### # Load previous simulation results ############################################################################### @@ -37,31 +19,34 @@ previous_results = openmc.deplete.ResultsList("depletion_results.h5") ############################################################################### # Instantiate a Settings object, set all runtime parameters -settings_file = openmc.Settings() -settings_file.batches = batches -settings_file.inactive = inactive -settings_file.particles = particles +settings = openmc.Settings() +settings.batches = 100 +settings.inactive = 10 +settings.particles = 10000 # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) -settings_file.source = openmc.source.Source(space=uniform_dist) +settings.source = openmc.source.Source(space=uniform_dist) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] entropy_mesh.dimension = [10, 10, 1] -settings_file.entropy_mesh = entropy_mesh +settings.entropy_mesh = entropy_mesh ############################################################################### # Initialize and run depletion calculation ############################################################################### -op = openmc.deplete.Operator(geometry, settings_file, chain_file, - previous_results) +# Create depletion "operator" +chain_file = './chain_simple.xml' +op = openmc.deplete.Operator(geometry, settings, chain_file, previous_results) # Perform simulation using the predictor algorithm -integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power) +time_steps = [1.0, 1.0, 1.0, 1.0, 1.0] # days +power = 174 # W/cm, for 2D simulations only (use W for 3D) +integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d') integrator.integrate() ############################################################################### @@ -77,27 +62,28 @@ time, keff = results.get_eigenvalue() # Obtain U235 concentration as a function of time time, n_U235 = results.get_atoms('1', 'U235') -# Obtain Xe135 absorption as a function of time -time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') +# Obtain Xe135 capture reaction rate as a function of time +time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') ############################################################################### # Generate plots ############################################################################### +days = 24*60*60 plt.figure() -plt.plot(time/(24*60*60), keff, label="K-effective") +plt.plot(time/days, keff, label="K-effective") plt.xlabel("Time (days)") plt.ylabel("Keff") plt.show() plt.figure() -plt.plot(time/(24*60*60), n_U235, label="U 235") +plt.plot(time/days, n_U235, label="U 235") plt.xlabel("Time (days)") plt.ylabel("n U5 (-)") plt.show() plt.figure() -plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption") +plt.plot(time/days, Xe_capture, label="Xe135 capture") plt.xlabel("Time (days)") plt.ylabel("RR (-)") plt.show() diff --git a/examples/python/pincell_depletion/run_depletion.py b/examples/python/pincell_depletion/run_depletion.py index 9933edd48a..ce4b0cf5e6 100644 --- a/examples/python/pincell_depletion/run_depletion.py +++ b/examples/python/pincell_depletion/run_depletion.py @@ -1,47 +1,31 @@ +from math import pi + import openmc import openmc.deplete -import numpy as np import matplotlib.pyplot as plt -############################################################################### -# Simulation Input File Parameters -############################################################################### - -# OpenMC simulation parameters -batches = 100 -inactive = 10 -particles = 1000 - -# Depletion simulation parameters -time_step = 1*24*60*60 # s -final_time = 5*24*60*60 # s -time_steps = np.full(final_time // time_step, time_step) -chain_file = './chain_simple.xml' -power = 174 # W/cm, for 2D simulations only (use W for 3D) - ############################################################################### # Define materials ############################################################################### # Instantiate some Materials and register the appropriate Nuclides -uo2 = openmc.Material(material_id=1, name='UO2 fuel at 2.4% wt enrichment') +uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment') uo2.set_density('g/cm3', 10.29769) uo2.add_element('U', 1., enrichment=2.4) uo2.add_element('O', 2.) -uo2.depletable = True -helium = openmc.Material(material_id=2, name='Helium for gap') +helium = openmc.Material(name='Helium for gap') helium.set_density('g/cm3', 0.001598) helium.add_element('He', 2.4044e-4) -zircaloy = openmc.Material(material_id=3, name='Zircaloy 4') +zircaloy = openmc.Material(name='Zircaloy 4') zircaloy.set_density('g/cm3', 6.55) -zircaloy.add_element('Sn', 0.014 , 'wo') +zircaloy.add_element('Sn', 0.014, 'wo') zircaloy.add_element('Fe', 0.00165, 'wo') -zircaloy.add_element('Cr', 0.001 , 'wo') +zircaloy.add_element('Cr', 0.001, 'wo') zircaloy.add_element('Zr', 0.98335, 'wo') -borated_water = openmc.Material(material_id=4, name='Borated water') +borated_water = openmc.Material(name='Borated water') borated_water.set_density('g/cm3', 0.740582) borated_water.add_element('B', 4.0e-5) borated_water.add_element('H', 5.0e-2) @@ -52,87 +36,62 @@ borated_water.add_s_alpha_beta('c_H_in_H2O') # Create geometry ############################################################################### -# Instantiate ZCylinder surfaces -fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.39218, name='Fuel OR') -clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=0.40005, name='Clad IR') -clad_or = openmc.ZCylinder(surface_id=3, x0=0, y0=0, r=0.45720, name='Clad OR') -left = openmc.XPlane(surface_id=4, x0=-0.62992, name='left') -right = openmc.XPlane(surface_id=5, x0=0.62992, name='right') -bottom = openmc.YPlane(surface_id=6, y0=-0.62992, name='bottom') -top = openmc.YPlane(surface_id=7, y0=0.62992, name='top') +# Define surfaces +pitch = 1.25984 +fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR') +clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR') +clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR') +box = openmc.model.rectangular_prism(pitch, pitch, boundary_type='reflective') -left.boundary_type = 'reflective' -right.boundary_type = 'reflective' -top.boundary_type = 'reflective' -bottom.boundary_type = 'reflective' +# Define cells +fuel = openmc.Cell(fill=uo2, region=-fuel_or) +gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir) +clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or) +water = openmc.Cell(fill=borated_water, region=+clad_or & box) -# Instantiate Cells -fuel = openmc.Cell(cell_id=1, name='cell 1') -gap = openmc.Cell(cell_id=2, name='cell 2') -clad = openmc.Cell(cell_id=3, name='cell 3') -water = openmc.Cell(cell_id=4, name='cell 4') - -# Use surface half-spaces to define regions -fuel.region = -fuel_or -gap.region = +fuel_or & -clad_ir -clad.region = +clad_ir & -clad_or -water.region = +clad_or & +left & -right & +bottom & -top - -# Register Materials with Cells -fuel.fill = uo2 -gap.fill = helium -clad.fill = zircaloy -water.fill = borated_water - -# Instantiate Universe -root = openmc.Universe(universe_id=0, name='root universe') - -# Register Cells with Universe -root.add_cells([fuel, gap, clad, water]) - -# Instantiate a Geometry, register the root Universe -geometry = openmc.Geometry(root) +# Define overall geometry +geometry = openmc.Geometry([fuel, gap, clad, water]) ############################################################################### # Set volumes of depletable materials ############################################################################### -# Compute cell areas -area = {} -area[fuel] = np.pi * fuel_or.coefficients['r'] ** 2 - -# Set materials volume for depletion. Set to an area for 2D simulations -uo2.volume = area[fuel] +# Set material volume for depletion. For 2D simulations, this should be an area. +uo2.volume = pi * fuel_or.r**2 ############################################################################### # Transport calculation settings ############################################################################### # Instantiate a Settings object, set all runtime parameters, and export to XML -settings_file = openmc.Settings() -settings_file.batches = batches -settings_file.inactive = inactive -settings_file.particles = particles +settings = openmc.Settings() +settings.batches = 100 +settings.inactive = 10 +settings.particles = 1000 # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) -settings_file.source = openmc.source.Source(space=uniform_dist) +settings.source = openmc.source.Source(space=uniform_dist) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50] entropy_mesh.dimension = [10, 10, 1] -settings_file.entropy_mesh = entropy_mesh +settings.entropy_mesh = entropy_mesh ############################################################################### # Initialize and run depletion calculation ############################################################################### -op = openmc.deplete.Operator(geometry, settings_file, chain_file) +# Create depletion "operator" +chain_file = './chain_simple.xml' +op = openmc.deplete.Operator(geometry, settings, chain_file) # Perform simulation using the predictor algorithm -integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power) +time_steps = [1.0, 1.0, 1.0, 1.0, 1.0] # days +power = 174 # W/cm, for 2D simulations only (use W for 3D) +integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d') integrator.integrate() ############################################################################### @@ -148,27 +107,28 @@ time, keff = results.get_eigenvalue() # Obtain U235 concentration as a function of time time, n_U235 = results.get_atoms('1', 'U235') -# Obtain Xe135 absorption as a function of time -time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') +# Obtain Xe135 capture reaction rate as a function of time +time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)') ############################################################################### # Generate plots ############################################################################### +days = 24*60*60 plt.figure() -plt.plot(time/(24*60*60), keff, label="K-effective") +plt.plot(time/days, keff, label="K-effective") plt.xlabel("Time (days)") plt.ylabel("Keff") plt.show() plt.figure() -plt.plot(time/(24*60*60), n_U235, label="U 235") +plt.plot(time/days, n_U235, label="U235") plt.xlabel("Time (days)") plt.ylabel("n U5 (-)") plt.show() plt.figure() -plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption") +plt.plot(time/days, Xe_capture, label="Xe135 capture") plt.xlabel("Time (days)") plt.ylabel("RR (-)") plt.show() From 061fc82e4374f0d810acf017236216e77eced683 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 12:30:19 +0000 Subject: [PATCH 14/59] These fixes allow the Cmake find_package infrastructure to work correctly --- CMakeLists.txt | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ec6b500b45..3cbd27b4e3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -207,7 +207,6 @@ endif() add_library(faddeeva STATIC vendor/faddeeva/Faddeeva.cc) target_include_directories(faddeeva PUBLIC - $ $ ) target_compile_options(faddeeva PRIVATE ${cxxflags}) @@ -389,7 +388,7 @@ add_custom_command(TARGET libopenmc POST_BUILD #=============================================================================== set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC) -install(TARGETS openmc libopenmc faddeeva +install(TARGETS openmc libopenmc faddeeva pugixml gsl-lite EXPORT openmc-targets RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} From aa9acca56ad1d8e43d6774afda37a6d4b8d70a32 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 22:10:41 +0000 Subject: [PATCH 15/59] Revert "These fixes allow the Cmake find_package infrastructure to work correctly" This reverts commit cfee500608a823cdbed751ed4cd7675cef519561. --- CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 3cbd27b4e3..ec6b500b45 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -207,6 +207,7 @@ endif() add_library(faddeeva STATIC vendor/faddeeva/Faddeeva.cc) target_include_directories(faddeeva PUBLIC + $ $ ) target_compile_options(faddeeva PRIVATE ${cxxflags}) @@ -388,7 +389,7 @@ add_custom_command(TARGET libopenmc POST_BUILD #=============================================================================== set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC) -install(TARGETS openmc libopenmc faddeeva pugixml gsl-lite +install(TARGETS openmc libopenmc faddeeva EXPORT openmc-targets RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} From b9c3d7af5424d241be522bcf8f791c9e1dfa4ce3 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 22:11:58 +0000 Subject: [PATCH 16/59] Fix CMake issues per Paul R fixes --- cmake/OpenMCConfig.cmake | 3 +++ 1 file changed, 3 insertions(+) diff --git a/cmake/OpenMCConfig.cmake b/cmake/OpenMCConfig.cmake index 0bc86fa71c..29a0e4542f 100644 --- a/cmake/OpenMCConfig.cmake +++ b/cmake/OpenMCConfig.cmake @@ -1,5 +1,8 @@ get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY) +find_package(fmt REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../fmt) +find_package(gsl-lite REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../gsl-lite) +find_package(pugixml REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../pugixml) find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl) find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor) From 8154eeed13db2a663464e9a204749a6d0bb67c10 Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Thu, 13 Feb 2020 14:41:40 +0000 Subject: [PATCH 17/59] Revert "Modified the StatePoint output to include the new max_lost_particles and relative version. Further modified an existing regression test to check the StatePoint write-out of these variables" This reverts commit b6a03949c82aee90824b9daf60d3d74a587b5b19. --- src/state_point.cpp | 4 ---- tests/regression_tests/photon_source/results_true.dat | 2 -- tests/regression_tests/photon_source/test.py | 2 -- 3 files changed, 8 deletions(-) diff --git a/src/state_point.cpp b/src/state_point.cpp index a9213bcb4e..60fd3826de 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -93,8 +93,6 @@ openmc_statepoint_write(const char* filename, bool* write_source) write_attribute(file_id, "photon_transport", settings::photon_transport); write_dataset(file_id, "n_particles", settings::n_particles); write_dataset(file_id, "n_batches", settings::n_batches); - write_dataset(file_id, "n_max_lost_particles", settings::n_max_lost_particles); - write_dataset(file_id, "relative_max_lost_particles", settings::relative_max_lost_particles); // Write out current batch number write_dataset(file_id, "current_batch", simulation::current_batch); @@ -375,8 +373,6 @@ void load_state_point() read_dataset(file_id, "n_particles", settings::n_particles); int temp; read_dataset(file_id, "n_batches", temp); - read_dataset(file_id, "n_max_lost_particles", settings::n_max_lost_particles); - read_dataset(file_id, "relative_max_lost_particles", settings::relative_max_lost_particles); // Take maximum of statepoint n_batches and input n_batches settings::n_batches = std::max(settings::n_batches, temp); diff --git a/tests/regression_tests/photon_source/results_true.dat b/tests/regression_tests/photon_source/results_true.dat index ab6ac9ff83..1448ad6db8 100644 --- a/tests/regression_tests/photon_source/results_true.dat +++ b/tests/regression_tests/photon_source/results_true.dat @@ -1,5 +1,3 @@ tally 1: sum = 2.275713E+02 sum_sq = 5.178870E+04 -max_lost = 5.000000E+00 -rel_max_lost = 1.000000E-01 diff --git a/tests/regression_tests/photon_source/test.py b/tests/regression_tests/photon_source/test.py index c5c42d680e..09a0b4a2d0 100644 --- a/tests/regression_tests/photon_source/test.py +++ b/tests/regression_tests/photon_source/test.py @@ -54,8 +54,6 @@ class SourceTestHarness(PyAPITestHarness): outstr += 'tally {}:\n'.format(t.id) outstr += 'sum = {:12.6E}\n'.format(t.sum[0, 0, 0]) outstr += 'sum_sq = {:12.6E}\n'.format(t.sum_sq[0, 0, 0]) - outstr += 'max_lost = {:12.6E}\n'.format(sp.n_max_lost_particles) - outstr += 'rel_max_lost = {:12.6E}\n'.format(sp.relative_max_lost_particles) return outstr From 68312dd9818cecd8129206b4a1a9d816432169e5 Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Thu, 13 Feb 2020 14:43:31 +0000 Subject: [PATCH 18/59] Revert "Modified an existing regression test to include the new max_lost_particles and relative version. This is mainly a test that the python/xml input/output routines work through the comparison the regression test harness does in this test, neither of these trigger the end of the simulation." This reverts commit 9224ca318675fcbb06c8eb97d9569466874dda21. --- tests/regression_tests/photon_source/inputs_true.dat | 2 -- tests/regression_tests/photon_source/test.py | 2 -- 2 files changed, 4 deletions(-) diff --git a/tests/regression_tests/photon_source/inputs_true.dat b/tests/regression_tests/photon_source/inputs_true.dat index 83c909ac69..89f4de0e0c 100644 --- a/tests/regression_tests/photon_source/inputs_true.dat +++ b/tests/regression_tests/photon_source/inputs_true.dat @@ -18,8 +18,6 @@ fixed source 10000 1 - 5 - 0.1 0 0 0 diff --git a/tests/regression_tests/photon_source/test.py b/tests/regression_tests/photon_source/test.py index 09a0b4a2d0..94b0280031 100644 --- a/tests/regression_tests/photon_source/test.py +++ b/tests/regression_tests/photon_source/test.py @@ -31,8 +31,6 @@ class SourceTestHarness(PyAPITestHarness): settings = openmc.Settings() settings.particles = 10000 settings.batches = 1 - settings.max_lost_particles = 5 - settings.rel_max_lost_particles = 0.1 settings.photon_transport = True settings.electron_treatment = 'ttb' settings.cutoff = {'energy_photon' : 1000.0} From 18b469bfa5cc4b4d1280655e9166c997db267481 Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Thu, 13 Feb 2020 14:44:53 +0000 Subject: [PATCH 19/59] Removed new variables from the statepoint python file as they are no longer read out --- openmc/statepoint.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 0d78bf68dc..a8e13efd25 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -77,10 +77,6 @@ class StatePoint(object): Number of batches n_inactive : int Number of inactive batches - n_max_lost_particles : int - Number of max lost particles - relative_max_lost_particles : float - Number of max lost particles, relative to the total number of particles n_particles : int Number of particles per generation n_realizations : int @@ -316,10 +312,6 @@ class StatePoint(object): else: return None - @property - def n_max_lost_particles(self): - return self._f['n_max_lost_particles'][()] - @property def relative_max_lost_particles(self): return self._f['relative_max_lost_particles'][()] From a9d1e37aba4aeadb683947861bb0a392c364e878 Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Thu, 13 Feb 2020 14:46:34 +0000 Subject: [PATCH 20/59] Changed the name of the max_lost_particles in the C++ layer to match the python side --- include/openmc/settings.h | 2 +- openmc/lib/settings.py | 2 +- src/particle.cpp | 2 +- src/settings.cpp | 6 +++--- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 3407eadf4b..8cee3c72ab 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -65,7 +65,7 @@ extern "C" std::string path_statepoint; //!< path to a statepoint file extern "C" int32_t n_batches; //!< number of (inactive+active) batches extern "C" int32_t n_inactive; //!< number of inactive batches -extern "C" int32_t n_max_lost_particles; //!< maximum number of lost particles +extern "C" int32_t max_lost_particles; //!< maximum number of lost particles extern double relative_max_lost_particles; //!< maximum number of lost particles, relative to the total number of particles extern "C" int32_t gen_per_batch; //!< number of generations per batch extern "C" int64_t n_particles; //!< number of particles per generation diff --git a/openmc/lib/settings.py b/openmc/lib/settings.py index 8345e00881..fa72e13234 100644 --- a/openmc/lib/settings.py +++ b/openmc/lib/settings.py @@ -22,7 +22,7 @@ class _Settings(object): entropy_on = _DLLGlobal(c_bool, 'entropy_on') generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch') inactive = _DLLGlobal(c_int32, 'n_inactive') - max_lost_particles = _DLLGlobal(c_int32, 'n_max_lost_particles') + max_lost_particles = _DLLGlobal(c_int32, 'max_lost_particles') rel_max_lost_particles = _DLLGlobal(c_double, 'relative_max_lost_particles') particles = _DLLGlobal(c_int64, 'n_particles') restart_run = _DLLGlobal(c_bool, 'restart_run') diff --git a/src/particle.cpp b/src/particle.cpp index 79e6625771..7e98b2ce67 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -630,7 +630,7 @@ Particle::mark_as_lost(const char* message) // Abort the simulation if the maximum number of lost particles has been // reached - if (simulation::n_lost_particles >= settings::n_max_lost_particles && + if (simulation::n_lost_particles >= settings::max_lost_particles && simulation::n_lost_particles >= settings::relative_max_lost_particles*n) { fatal_error("Maximum number of lost particles has been reached."); } diff --git a/src/settings.cpp b/src/settings.cpp index d72e5f17f7..933f9a2172 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -79,7 +79,7 @@ std::string path_statepoint; int32_t n_batches; int32_t n_inactive {0}; -int32_t n_max_lost_particles {10}; +int32_t max_lost_particles {10}; double relative_max_lost_particles {1.0e-6}; int32_t gen_per_batch {1}; int64_t n_particles {-1}; @@ -147,7 +147,7 @@ void get_run_parameters(pugi::xml_node node_base) // Get max number of lost particles if (check_for_node(node_base, "max_lost_particles")) { - n_max_lost_particles = std::stoi(get_node_value(node_base, "max_lost_particles")); + max_lost_particles = std::stoi(get_node_value(node_base, "max_lost_particles")); } // Get relative number of lost particles @@ -362,7 +362,7 @@ void read_settings_xml() fatal_error("Number of inactive batches must be non-negative."); } else if (n_particles <= 0) { fatal_error("Number of particles must be greater than zero."); - } else if (n_max_lost_particles <= 0) { + } else if (max_lost_particles <= 0) { fatal_error("Number of max lost particles must be greater than zero."); } else if (relative_max_lost_particles <= 0.0 || relative_max_lost_particles >= 1.0) { fatal_error("Relative max lost particles must be between zero and one."); From a9fc400863952502bad61fd73bb80993448afc74 Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Thu, 13 Feb 2020 14:51:32 +0000 Subject: [PATCH 21/59] Also removed relative_max_lost_particles from statepoint python as it is not being written out --- openmc/statepoint.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index a8e13efd25..84512d70e7 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -310,11 +310,7 @@ class StatePoint(object): if self.run_mode == 'eigenvalue': return self._f['n_inactive'][()] else: - return None - - @property - def relative_max_lost_particles(self): - return self._f['relative_max_lost_particles'][()] + return None @property def n_particles(self): From d763806c43ea05adb25831b8fd11c828fd40d06a Mon Sep 17 00:00:00 2001 From: stevendargaville Date: Thu, 13 Feb 2020 14:52:44 +0000 Subject: [PATCH 22/59] Also made rel_max_lost_particle name consistent between C++ and python layers --- include/openmc/settings.h | 2 +- openmc/lib/settings.py | 2 +- src/particle.cpp | 2 +- src/settings.cpp | 6 +++--- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 8cee3c72ab..2b4203646f 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -66,7 +66,7 @@ extern "C" std::string path_statepoint; //!< path to a statepoint file extern "C" int32_t n_batches; //!< number of (inactive+active) batches extern "C" int32_t n_inactive; //!< number of inactive batches extern "C" int32_t max_lost_particles; //!< maximum number of lost particles -extern double relative_max_lost_particles; //!< maximum number of lost particles, relative to the total number of particles +extern double rel_max_lost_particles; //!< maximum number of lost particles, relative to the total number of particles extern "C" int32_t gen_per_batch; //!< number of generations per batch extern "C" int64_t n_particles; //!< number of particles per generation diff --git a/openmc/lib/settings.py b/openmc/lib/settings.py index fa72e13234..68d68c29e0 100644 --- a/openmc/lib/settings.py +++ b/openmc/lib/settings.py @@ -23,7 +23,7 @@ class _Settings(object): generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch') inactive = _DLLGlobal(c_int32, 'n_inactive') max_lost_particles = _DLLGlobal(c_int32, 'max_lost_particles') - rel_max_lost_particles = _DLLGlobal(c_double, 'relative_max_lost_particles') + rel_max_lost_particles = _DLLGlobal(c_double, 'rel_max_lost_particles') particles = _DLLGlobal(c_int64, 'n_particles') restart_run = _DLLGlobal(c_bool, 'restart_run') run_CE = _DLLGlobal(c_bool, 'run_CE') diff --git a/src/particle.cpp b/src/particle.cpp index 7e98b2ce67..8f6dba803f 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -631,7 +631,7 @@ Particle::mark_as_lost(const char* message) // Abort the simulation if the maximum number of lost particles has been // reached if (simulation::n_lost_particles >= settings::max_lost_particles && - simulation::n_lost_particles >= settings::relative_max_lost_particles*n) { + simulation::n_lost_particles >= settings::rel_max_lost_particles*n) { fatal_error("Maximum number of lost particles has been reached."); } } diff --git a/src/settings.cpp b/src/settings.cpp index 933f9a2172..706a171458 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -80,7 +80,7 @@ std::string path_statepoint; int32_t n_batches; int32_t n_inactive {0}; int32_t max_lost_particles {10}; -double relative_max_lost_particles {1.0e-6}; +double rel_max_lost_particles {1.0e-6}; int32_t gen_per_batch {1}; int64_t n_particles {-1}; @@ -152,7 +152,7 @@ void get_run_parameters(pugi::xml_node node_base) // Get relative number of lost particles if (check_for_node(node_base, "rel_max_lost_particles")) { - relative_max_lost_particles = std::stod(get_node_value(node_base, "rel_max_lost_particles")); + rel_max_lost_particles = std::stod(get_node_value(node_base, "rel_max_lost_particles")); } // Get number of inactive batches @@ -364,7 +364,7 @@ void read_settings_xml() fatal_error("Number of particles must be greater than zero."); } else if (max_lost_particles <= 0) { fatal_error("Number of max lost particles must be greater than zero."); - } else if (relative_max_lost_particles <= 0.0 || relative_max_lost_particles >= 1.0) { + } else if (rel_max_lost_particles <= 0.0 || rel_max_lost_particles >= 1.0) { fatal_error("Relative max lost particles must be between zero and one."); } } From 39672e79b0fcf7e7cbbd63d7a9dad14f7251b459 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Thu, 13 Jun 2019 09:15:21 +0100 Subject: [PATCH 23/59] Added internal infrastructure for dlopen based shared object sources --- include/openmc/settings.h | 1 + include/openmc/source.h | 4 +++ src/settings.cpp | 1 + src/source.cpp | 54 ++++++++++++++++++++++++++++++++++++++- 4 files changed, 59 insertions(+), 1 deletion(-) diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 2b4203646f..494cfd1a20 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -60,6 +60,7 @@ extern std::string path_input; //!< directory where main .xml files r extern std::string path_output; //!< directory where output files are written extern std::string path_particle_restart; //!< path to a particle restart file extern std::string path_source; +extern std::string path_source_library; //!< path to the source shared object extern std::string path_sourcepoint; //!< path to a source file extern "C" std::string path_statepoint; //!< path to a statepoint file diff --git a/include/openmc/source.h b/include/openmc/source.h index a177995ea8..0a73260fd1 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -59,6 +59,10 @@ private: //! Initialize source bank from file/distribution extern "C" void initialize_source(); +// as yet uncreated function to sample a source from a shared object +// +// extern "C" Particle::Bank sample_source(); + //! Sample a site from all external source distributions in proportion to their //! source strength //! \param[inout] seed Pseudorandom seed pointer diff --git a/src/settings.cpp b/src/settings.cpp index 7808a9d0f6..614dd8f5b4 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -74,6 +74,7 @@ std::string path_input; std::string path_output; std::string path_particle_restart; std::string path_source; +std::string path_source_library; std::string path_sourcepoint; std::string path_statepoint; diff --git a/src/source.cpp b/src/source.cpp index b15a15d23a..4666783d4a 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -1,6 +1,8 @@ #include "openmc/source.h" #include // for move +#include // for stringstream +#include // for dlopen #include #include "xtensor/xadapt.hpp" @@ -71,7 +73,14 @@ SourceDistribution::SourceDistribution(pugi::xml_node node) fatal_error(fmt::format("Source file '{}' does not exist.", settings::path_source)); } - + } else if (check_for_node(node, "library")) { + settings::path_source_library = get_node_value(node, "library", false, true); + // check if it exists + if (!file_exists(settings::path_source_library)) { + std::stringstream msg; + msg << "Library file " << settings::path_source_library << "' does not exist."; + fatal_error(msg); + } } else { // Spatial distribution for external source @@ -261,7 +270,50 @@ void initialize_source() // Close file file_close(file_id); + } else if ( settings::path_source_library != "" ) { + // Get the source from a library object + std::stringstream msg; + msg << "Sampling from library source " << settings::path_source << "..."; + write_message(msg, 6); + + // Open the library + void* source_library = dlopen(settings::path_source_library.c_str(),RTLD_LAZY); + if(!source_library) { + std::stringstream msg("Couldnt open source library " + settings::path_source_library); + fatal_error(msg); + } + + // load the symbol + typedef Particle::Bank (*sample_t)(); + + // reset errors + dlerror(); + + // get the function from the library + sample_t sample_source = (sample_t) dlsym(source_library, "sample_source"); + const char *dlsym_error = dlerror(); + + if (dlsym_error) { + dlclose(source_library); + fatal_error(dlsym_error); + } + + // Generation source sites from specified distribution in the + // library source + for (int64_t i = 0; i < simulation::work_per_rank; ++i) { + // initialize random number seed + int64_t id = simulation::total_gen*settings::n_particles + + simulation::work_index[mpi::rank] + i + 1; + set_particle_seed(id); + + // sample external source distribution + simulation::source_bank[i] = sample_source(); + } + + // release the library + dlclose(source_library); + } else { // Generation source sites from specified distribution in user input for (int64_t i = 0; i < simulation::work_per_rank; ++i) { From e1baa804a639527654e4e2bf226014abb7e2a816 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Thu, 13 Jun 2019 14:40:15 +0100 Subject: [PATCH 24/59] Added missing fixed source subroutine --- src/source.cpp | 50 ++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 46 insertions(+), 4 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 4666783d4a..b856a8d2ad 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -294,9 +294,11 @@ void initialize_source() sample_t sample_source = (sample_t) dlsym(source_library, "sample_source"); const char *dlsym_error = dlerror(); + // check for any dlsym errors if (dlsym_error) { + std::cout << dlsym_error << std::endl; dlclose(source_library); - fatal_error(dlsym_error); + fatal_error("Couldnt open the sample_source symbol"); } // Generation source sites from specified distribution in the @@ -313,7 +315,7 @@ void initialize_source() // release the library dlclose(source_library); - + } else { // Generation source sites from specified distribution in user input for (int64_t i = 0; i < simulation::work_per_rank; ++i) { @@ -375,8 +377,7 @@ void free_memory_source() void fill_source_bank_fixedsource() { - if (settings::path_source.empty()) { - #pragma omp parallel for + if (settings::path_source.empty() && settings::path_source_library.empty()) { for (int64_t i = 0; i < simulation::work_per_rank; ++i) { // initialize random number seed int64_t id = (simulation::total_gen + overall_generation()) * @@ -386,6 +387,47 @@ void fill_source_bank_fixedsource() // sample external source distribution simulation::source_bank[i] = sample_external_source(&seed); } + } else if (settings::path_source.empty() && !settings::path_source.empty()) { + std::stringstream msg; + + // Open the library + void* source_library = dlopen(settings::path_source_library.c_str(),RTLD_LAZY); + if(!source_library) { + std::stringstream msg("Couldnt open source library " + settings::path_source_library); + fatal_error(msg); + } + + // load the symbol + typedef Particle::Bank (*sample_t)(); + + // reset errors + dlerror(); + + // get the function from the library + sample_t sample_source = (sample_t) dlsym(source_library, "sample_source"); + const char *dlsym_error = dlerror(); + + // check for any dlsym errors + if (dlsym_error) { + std::cout << dlsym_error << std::endl; + dlclose(source_library); + fatal_error("Couldnt open the sample_source symbol"); + } + + // Generation source sites from specified distribution in the + // library source + for (int64_t i = 0; i < simulation::work_per_rank; ++i) { + // initialize random number seed + int64_t id = simulation::total_gen*settings::n_particles + + simulation::work_index[mpi::rank] + i + 1; + set_particle_seed(id); + + // sample external source distribution + simulation::source_bank[i] = sample_source(); + } + + // release the library + dlclose(source_library); } } From 995a920854a56ae273fcadad3f20b3b25cc86728 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Thu, 13 Jun 2019 17:49:29 +0100 Subject: [PATCH 25/59] Added example on how to use xml based custom source --- examples/xml/custom_source/geometry.xml | 15 +++++++++++++ examples/xml/custom_source/materials.xml | 16 ++++++++++++++ examples/xml/custom_source/settings.xml | 15 +++++++++++++ examples/xml/custom_source/source_ring.cpp | 25 ++++++++++++++++++++++ examples/xml/custom_source/tallies.xml | 17 +++++++++++++++ 5 files changed, 88 insertions(+) create mode 100644 examples/xml/custom_source/geometry.xml create mode 100644 examples/xml/custom_source/materials.xml create mode 100644 examples/xml/custom_source/settings.xml create mode 100644 examples/xml/custom_source/source_ring.cpp create mode 100644 examples/xml/custom_source/tallies.xml diff --git a/examples/xml/custom_source/geometry.xml b/examples/xml/custom_source/geometry.xml new file mode 100644 index 0000000000..b30884f8ca --- /dev/null +++ b/examples/xml/custom_source/geometry.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/examples/xml/custom_source/materials.xml b/examples/xml/custom_source/materials.xml new file mode 100644 index 0000000000..606c676df8 --- /dev/null +++ b/examples/xml/custom_source/materials.xml @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/examples/xml/custom_source/settings.xml b/examples/xml/custom_source/settings.xml new file mode 100644 index 0000000000..f8d497459c --- /dev/null +++ b/examples/xml/custom_source/settings.xml @@ -0,0 +1,15 @@ + + + + fixed source + 10 + 0 + 100000 + + + + + ./source_ring.so + + + diff --git a/examples/xml/custom_source/source_ring.cpp b/examples/xml/custom_source/source_ring.cpp new file mode 100644 index 0000000000..ef2784ea36 --- /dev/null +++ b/examples/xml/custom_source/source_ring.cpp @@ -0,0 +1,25 @@ +#include +#include "openmc/random_lcg.h" +#include "openmc/source.h" +#include "openmc/particle.h" + +// you must have external C linkage here otherwise +// dlopen will not find the file +extern "C" openmc::Particle::Bank sample_source() { + openmc::Particle::Bank particle; + // wgt + particle.particle = openmc::Particle::Type::neutron; + particle.wgt = 1.0; + // position + + double angle = 2.*M_PI*openmc::prn(); + double radius = 3.0; + particle.r.x = radius*std::cos(angle); + particle.r.y = radius*std::sin(angle); + particle.r.z = 0.; + // angle + particle.u = {1.,0,0}; + particle.E = 14.08e6; + particle.delayed_group = 0; + return particle; +} \ No newline at end of file diff --git a/examples/xml/custom_source/tallies.xml b/examples/xml/custom_source/tallies.xml new file mode 100644 index 0000000000..7f6f299261 --- /dev/null +++ b/examples/xml/custom_source/tallies.xml @@ -0,0 +1,17 @@ + + + + + 100 + + + + 0 20.0e6 + + + + 1 2 + flux + + + From c3fb48bc62ff880b654a5fabef7c808f301943ef Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Fri, 14 Jun 2019 08:33:44 +0100 Subject: [PATCH 26/59] Added python bindngs for library based source --- openmc/source.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/openmc/source.py b/openmc/source.py index 88c2f86119..e0fd2b6c40 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -44,11 +44,12 @@ class Source(object): """ def __init__(self, space=None, angle=None, energy=None, filename=None, - strength=1.0, particle='neutron'): + library=None, strength=1.0, particle='neutron'): self._space = None self._angle = None self._energy = None self._file = None + self._source_library = None if space is not None: self.space = space @@ -58,6 +59,8 @@ class Source(object): self.energy = energy if filename is not None: self.file = filename + if library is not None: + self.source_library = library self.strength = strength self.particle = particle @@ -65,6 +68,10 @@ class Source(object): def file(self): return self._file + @property + def library(self): + return self._source_library + @property def space(self): return self._space @@ -90,6 +97,11 @@ class Source(object): cv.check_type('source file', filename, str) self._file = filename + @library.setter + def library(self, library_name): + cv.check_type('library', library_name, str) + self._source_library = library_name + @space.setter def space(self, space): cv.check_type('spatial distribution', space, Spatial) @@ -131,6 +143,8 @@ class Source(object): element.set("particle", self.particle) if self.file is not None: element.set("file", self.file) + if self.source_library is not None: + element.set("library", self.source_library) if self.space is not None: element.append(self.space.to_xml_element()) if self.angle is not None: @@ -168,6 +182,10 @@ class Source(object): if filename is not None: source.file = filename + library = get_text(elem, 'library') + if library is not None: + source.source_library = library + space = elem.find('space') if space is not None: source.space = Spatial.from_xml_element(space) From d140b9e49978cdcbcdf3dfa110a76c48c28859a9 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Fri, 14 Jun 2019 10:52:55 +0100 Subject: [PATCH 27/59] Fix the python part --- openmc/source.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index e0fd2b6c40..f8ea8ea936 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -60,7 +60,7 @@ class Source(object): if filename is not None: self.file = filename if library is not None: - self.source_library = library + self.library = library self.strength = strength self.particle = particle @@ -143,8 +143,8 @@ class Source(object): element.set("particle", self.particle) if self.file is not None: element.set("file", self.file) - if self.source_library is not None: - element.set("library", self.source_library) + if self.library is not None: + element.set("library", self.library) if self.space is not None: element.append(self.space.to_xml_element()) if self.angle is not None: @@ -184,7 +184,7 @@ class Source(object): library = get_text(elem, 'library') if library is not None: - source.source_library = library + source.library = library space = elem.find('space') if space is not None: From 0324200b61d986103f23b6107e59aa0220234b88 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Thu, 13 Jun 2019 09:15:21 +0100 Subject: [PATCH 28/59] Added internal infrastructure for dlopen based shared object sources --- src/source.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index b856a8d2ad..089dbddef7 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -312,10 +312,9 @@ void initialize_source() // sample external source distribution simulation::source_bank[i] = sample_source(); } - // release the library dlclose(source_library); - + } else { // Generation source sites from specified distribution in user input for (int64_t i = 0; i < simulation::work_per_rank; ++i) { From 89f86883fba7c5a02ad8dcb3ac22286204cf305a Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Thu, 13 Jun 2019 14:40:15 +0100 Subject: [PATCH 29/59] Added missing fixed source subroutine --- src/source.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/source.cpp b/src/source.cpp index 089dbddef7..b7d1e24990 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -314,7 +314,7 @@ void initialize_source() } // release the library dlclose(source_library); - + } else { // Generation source sites from specified distribution in user input for (int64_t i = 0; i < simulation::work_per_rank; ++i) { From 7781b9249574acaa770eba61d057844951e2426b Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Fri, 14 Jun 2019 08:33:44 +0100 Subject: [PATCH 30/59] Added python bindngs for library based source --- openmc/source.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index f8ea8ea936..e0fd2b6c40 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -60,7 +60,7 @@ class Source(object): if filename is not None: self.file = filename if library is not None: - self.library = library + self.source_library = library self.strength = strength self.particle = particle @@ -143,8 +143,8 @@ class Source(object): element.set("particle", self.particle) if self.file is not None: element.set("file", self.file) - if self.library is not None: - element.set("library", self.library) + if self.source_library is not None: + element.set("library", self.source_library) if self.space is not None: element.append(self.space.to_xml_element()) if self.angle is not None: @@ -184,7 +184,7 @@ class Source(object): library = get_text(elem, 'library') if library is not None: - source.library = library + source.source_library = library space = elem.find('space') if space is not None: From 8c3d47a4e02fabd4ec2c9cc455f45983e3dca1d0 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Fri, 14 Jun 2019 10:52:55 +0100 Subject: [PATCH 31/59] Fix the python part --- openmc/source.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index e0fd2b6c40..f8ea8ea936 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -60,7 +60,7 @@ class Source(object): if filename is not None: self.file = filename if library is not None: - self.source_library = library + self.library = library self.strength = strength self.particle = particle @@ -143,8 +143,8 @@ class Source(object): element.set("particle", self.particle) if self.file is not None: element.set("file", self.file) - if self.source_library is not None: - element.set("library", self.source_library) + if self.library is not None: + element.set("library", self.library) if self.space is not None: element.append(self.space.to_xml_element()) if self.angle is not None: @@ -184,7 +184,7 @@ class Source(object): library = get_text(elem, 'library') if library is not None: - source.source_library = library + source.library = library space = elem.find('space') if space is not None: From 98fc03bda13314ee539ada19c9917d014ed79151 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Mon, 5 Aug 2019 22:44:21 +0100 Subject: [PATCH 32/59] Now writes a Makefile for users to build a source --- CMakeLists.txt | 58 +++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 55 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ec6b500b45..6cc9fceaf5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -406,6 +406,58 @@ install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1) install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright) install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) -# Copy headers for vendored dependencies (note that all except faddeeva are handled -# separately since they are managed by CMake) -install(DIRECTORY vendor/faddeeva DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) +#============================= +# write the dl_open Makefile +#============================= +file(WRITE share/Makefile" +# Makefile to build dynamic sources for OpenMC +# this assumes that your source sampling filename is +# source_sampling.cpp +# +# you can add fortran, c and cpp dependencies to this source +# adding them in FC_DEPS, C_DEPS, and CXX_DEPS accordingly + +ifeq ($(FC), f77) +FC = gfortran +endif + +default: all + +ALL = source_sampling +# add your fortran depencies here +FC_DEPS = +# add your c dependencies here +C_DEPS = +# add your cpp dependencies here +CPP_DEPS = +DEPS = $(FC_DEPS) $(C_DEPS) $(CPP_DEPS) +OPT_LEVEL = -O3 +FFLAGS = $(OPT_LEVEL) -fPIC +C_FLAGS = -fPIC +CXX_FLAGS = -fPIC + +OPENMC_DIR = ${CMAKE_INSTALL_PREFIX} +OPENMC_INC_DIR = $(OPENMC_DIR)/include +OPENMC_LIB_DIR = $(OPENMC_DIR)/lib +# setting the so name is important +LINK_FLAGS = $(OPT_LEVEL) -Wall -Wl,-soname,source_sampling.so +# setting shared is important +LINK_FLAGS += -L$(OPENMC_LIB_DIR) -lopenmc -lgfortran -fPIC -shared +OPENMC_INCLUDES = -I$(OPENMC_INC_DIR) -I$(OPENMC_DIR)/vendor/pugixml + +all: $(ALL) + +source_sampling: $(DEPS) + $(CXX) source_sampling.cpp $(DEPS) $(OPENMC_INCLUDES) $(LINK_FLAGS) -o $@.so +# make any fortran objects +%.o : %.F90 + $(FC) -c $(FFLAGS) $*.F90 -o $@ +# make any c objects +%.o : %.c + $(CC) -c $(FFLAGS) $*.c -o $@ +#make any cpp objects +%.o : %.cpp + $(CXX) -c $(FFLAGS) $*.cpp -o $@ +clean: + rm -rf *.o *.mod +") From 0bb36da36bb1410a3341c0c3f94f5bef9046bf1a Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Thu, 13 Jun 2019 09:15:21 +0100 Subject: [PATCH 33/59] Added internal infrastructure for dlopen based shared object sources --- src/source.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index b7d1e24990..51d86fcaa5 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -296,9 +296,8 @@ void initialize_source() // check for any dlsym errors if (dlsym_error) { - std::cout << dlsym_error << std::endl; dlclose(source_library); - fatal_error("Couldnt open the sample_source symbol"); + fatal_error(dlsym_error); } // Generation source sites from specified distribution in the @@ -314,7 +313,7 @@ void initialize_source() } // release the library dlclose(source_library); - + } else { // Generation source sites from specified distribution in user input for (int64_t i = 0; i < simulation::work_per_rank; ++i) { From da3c0e83f0f266a2908b3739c440b099ed777de0 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Thu, 13 Jun 2019 14:40:15 +0100 Subject: [PATCH 34/59] Added missing fixed source subroutine --- src/source.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 51d86fcaa5..b7d1e24990 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -296,8 +296,9 @@ void initialize_source() // check for any dlsym errors if (dlsym_error) { + std::cout << dlsym_error << std::endl; dlclose(source_library); - fatal_error(dlsym_error); + fatal_error("Couldnt open the sample_source symbol"); } // Generation source sites from specified distribution in the @@ -313,7 +314,7 @@ void initialize_source() } // release the library dlclose(source_library); - + } else { // Generation source sites from specified distribution in user input for (int64_t i = 0; i < simulation::work_per_rank; ++i) { From 590e659079c76b08902699e782d0cf89a68f93bd Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Fri, 14 Jun 2019 08:33:44 +0100 Subject: [PATCH 35/59] Added python bindngs for library based source --- openmc/source.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index f8ea8ea936..e0fd2b6c40 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -60,7 +60,7 @@ class Source(object): if filename is not None: self.file = filename if library is not None: - self.library = library + self.source_library = library self.strength = strength self.particle = particle @@ -143,8 +143,8 @@ class Source(object): element.set("particle", self.particle) if self.file is not None: element.set("file", self.file) - if self.library is not None: - element.set("library", self.library) + if self.source_library is not None: + element.set("library", self.source_library) if self.space is not None: element.append(self.space.to_xml_element()) if self.angle is not None: @@ -184,7 +184,7 @@ class Source(object): library = get_text(elem, 'library') if library is not None: - source.library = library + source.source_library = library space = elem.find('space') if space is not None: From 50c55041b3309d8496c57832b585db9a3a1289d9 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Fri, 14 Jun 2019 10:52:55 +0100 Subject: [PATCH 36/59] Fix the python part --- openmc/source.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index e0fd2b6c40..f8ea8ea936 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -60,7 +60,7 @@ class Source(object): if filename is not None: self.file = filename if library is not None: - self.source_library = library + self.library = library self.strength = strength self.particle = particle @@ -143,8 +143,8 @@ class Source(object): element.set("particle", self.particle) if self.file is not None: element.set("file", self.file) - if self.source_library is not None: - element.set("library", self.source_library) + if self.library is not None: + element.set("library", self.library) if self.space is not None: element.append(self.space.to_xml_element()) if self.angle is not None: @@ -184,7 +184,7 @@ class Source(object): library = get_text(elem, 'library') if library is not None: - source.source_library = library + source.library = library space = elem.find('space') if space is not None: From ceb88c2d17eb078bdc6587b8202e67c4943e0dc8 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 7 Aug 2019 08:02:17 +0100 Subject: [PATCH 37/59] missing space --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6cc9fceaf5..605ba0b612 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -409,7 +409,7 @@ install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) #============================= # write the dl_open Makefile #============================= -file(WRITE share/Makefile" +file(WRITE share/Makefile " # Makefile to build dynamic sources for OpenMC # this assumes that your source sampling filename is # source_sampling.cpp From 46c79c0889f1b94a27909e376d46e80102efd718 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 7 Aug 2019 08:02:38 +0100 Subject: [PATCH 38/59] Adding unit test for dlopen source --- .../dlopen_source/source_sampling.cpp | 23 ++++++++ .../dlopen_source/test_dlopen_source.py | 57 +++++++++++++++++++ 2 files changed, 80 insertions(+) create mode 100644 tests/unit_tests/dlopen_source/source_sampling.cpp create mode 100644 tests/unit_tests/dlopen_source/test_dlopen_source.py diff --git a/tests/unit_tests/dlopen_source/source_sampling.cpp b/tests/unit_tests/dlopen_source/source_sampling.cpp new file mode 100644 index 0000000000..4b13c7db92 --- /dev/null +++ b/tests/unit_tests/dlopen_source/source_sampling.cpp @@ -0,0 +1,23 @@ +#include +#include "openmc/random_lcg.h" +#include "openmc/source.h" +#include "openmc/particle.h" + +// you must have external C linkage here otherwise +// dlopen will not find the file +extern "C" openmc::Particle::Bank sample_source() { + openmc::Particle::Bank particle; + // wgt + particle.particle = openmc::Particle::Type::neutron; + particle.wgt = 1.0; + // position + + particle.r.x = 0.; + particle.r.y = 0.; + particle.r.z = 0.; + // angle + particle.u = {1.,0,0}; + particle.E = 14.08e6; + particle.delayed_group = 0; + return particle; +} diff --git a/tests/unit_tests/dlopen_source/test_dlopen_source.py b/tests/unit_tests/dlopen_source/test_dlopen_source.py new file mode 100644 index 0000000000..917f3fd4fe --- /dev/null +++ b/tests/unit_tests/dlopen_source/test_dlopen_source.py @@ -0,0 +1,57 @@ +import openmc +import pytest +import os + +# compile the external source +def compile_source(): + # needs a more robust way to know where the + # Makefile is + status = os.system('cp /usr/local/share/Makefile .') + assert os.WEXITSTATUS(status) == 0 + status = os.system('make') + assert os.WEXITSTATUS(status) == 0 + + return + +# build the test geometry +def make_geometry(): + mats = openmc.Materials() + + natural_lead = openmc.Material(1, "natural_lead") + natural_lead.add_element('Pb', 1,'ao') + natural_lead.set_density('g/cm3', 11.34) + mats.append(natural_lead) + + # surfaces + surface_sph1 = openmc.Sphere(r=100) + volume_sph1 = -surface_sph1 + + # cell + cell_1 = openmc.Cell(region=volume_sph1) + cell_1.fill = natural_lead #assigning a material to a cell + universe = openmc.Universe(cells=[cell_1]) #hint, this list will need to include the new cell + geom = openmc.Geometry(universe) + + # settings + sett = openmc.Settings() + batches = 10 + sett.batches = batches + sett.inactive = 0 + sett.particles = 1000 + sett.particle = "neutron" + sett.run_mode = 'fixed source' + + #source + source = openmc.Source() + source.library = PATH+'source_sampling.so' + sett.source = source + + # run + model = openmc.model.Model(geom,mats,sett) + return model + +def test_dlopen_source(): + compile_source() + model = make_geometry() + model.run() + From f387d2cba9437e41affbd8bbe5ea8a9799b34745 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Tue, 15 Oct 2019 20:41:10 +0100 Subject: [PATCH 39/59] make a cmakefile instead of a makefile for cross platform compatibility --- CMakeLists.txt | 59 +++++++------------------------------------------- 1 file changed, 8 insertions(+), 51 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 605ba0b612..18de941218 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -409,55 +409,12 @@ install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) #============================= # write the dl_open Makefile #============================= -file(WRITE share/Makefile " -# Makefile to build dynamic sources for OpenMC -# this assumes that your source sampling filename is -# source_sampling.cpp -# -# you can add fortran, c and cpp dependencies to this source -# adding them in FC_DEPS, C_DEPS, and CXX_DEPS accordingly - -ifeq ($(FC), f77) -FC = gfortran -endif - -default: all - -ALL = source_sampling -# add your fortran depencies here -FC_DEPS = -# add your c dependencies here -C_DEPS = -# add your cpp dependencies here -CPP_DEPS = -DEPS = $(FC_DEPS) $(C_DEPS) $(CPP_DEPS) -OPT_LEVEL = -O3 -FFLAGS = $(OPT_LEVEL) -fPIC -C_FLAGS = -fPIC -CXX_FLAGS = -fPIC - -OPENMC_DIR = ${CMAKE_INSTALL_PREFIX} -OPENMC_INC_DIR = $(OPENMC_DIR)/include -OPENMC_LIB_DIR = $(OPENMC_DIR)/lib -# setting the so name is important -LINK_FLAGS = $(OPT_LEVEL) -Wall -Wl,-soname,source_sampling.so -# setting shared is important -LINK_FLAGS += -L$(OPENMC_LIB_DIR) -lopenmc -lgfortran -fPIC -shared -OPENMC_INCLUDES = -I$(OPENMC_INC_DIR) -I$(OPENMC_DIR)/vendor/pugixml - -all: $(ALL) - -source_sampling: $(DEPS) - $(CXX) source_sampling.cpp $(DEPS) $(OPENMC_INCLUDES) $(LINK_FLAGS) -o $@.so -# make any fortran objects -%.o : %.F90 - $(FC) -c $(FFLAGS) $*.F90 -o $@ -# make any c objects -%.o : %.c - $(CC) -c $(FFLAGS) $*.c -o $@ -#make any cpp objects -%.o : %.cpp - $(CXX) -c $(FFLAGS) $*.cpp -o $@ -clean: - rm -rf *.o *.mod +file(WRITE ${CMAKE_INSTALL_PREFIX}/share/CMakeLists.txt " +cmake_minimum_required(VERSION 3.3 FATAL_ERROR) +project(openmc_sources C CXX) +add_library(source SHARED \$\{SOURCE_FILES\}) +target_include_directories(source PUBLIC ${PROJECT_SOURCE_DIR}/include + ${PROJECT_SOURCE_DIR}/vendor/pugixml ${HDF5_INCLUDE_DIRS}) +target_link_libraries(source ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) ") + From cba7e9bc87e1757d228aba9e29b271f8808d0f6c Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Tue, 15 Oct 2019 22:00:07 +0100 Subject: [PATCH 40/59] Added regression tests for dlopen source --- .../source_dlopen/inputs_true.dat | 23 ++++++ .../source_dlopen/results_true.dat | 0 .../source_dlopen}/source_sampling.cpp | 0 tests/regression_tests/source_dlopen/test.py | 70 +++++++++++++++++++ .../dlopen_source/test_dlopen_source.py | 57 --------------- 5 files changed, 93 insertions(+), 57 deletions(-) create mode 100644 tests/regression_tests/source_dlopen/inputs_true.dat create mode 100644 tests/regression_tests/source_dlopen/results_true.dat rename tests/{unit_tests/dlopen_source => regression_tests/source_dlopen}/source_sampling.cpp (100%) create mode 100644 tests/regression_tests/source_dlopen/test.py delete mode 100644 tests/unit_tests/dlopen_source/test_dlopen_source.py diff --git a/tests/regression_tests/source_dlopen/inputs_true.dat b/tests/regression_tests/source_dlopen/inputs_true.dat new file mode 100644 index 0000000000..23f3d84384 --- /dev/null +++ b/tests/regression_tests/source_dlopen/inputs_true.dat @@ -0,0 +1,23 @@ + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + 0 + + diff --git a/tests/regression_tests/source_dlopen/results_true.dat b/tests/regression_tests/source_dlopen/results_true.dat new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/unit_tests/dlopen_source/source_sampling.cpp b/tests/regression_tests/source_dlopen/source_sampling.cpp similarity index 100% rename from tests/unit_tests/dlopen_source/source_sampling.cpp rename to tests/regression_tests/source_dlopen/source_sampling.cpp diff --git a/tests/regression_tests/source_dlopen/test.py b/tests/regression_tests/source_dlopen/test.py new file mode 100644 index 0000000000..ac63ca8045 --- /dev/null +++ b/tests/regression_tests/source_dlopen/test.py @@ -0,0 +1,70 @@ +import openmc +import pytest +import os +import glob + +from tests.testing_harness import PyAPITestHarness + +# compile the external source +def compile_source(): + # first one should be CMakelists in share/ + files = glob.glob('../../../*/CMakeLists.txt') + assert len(files) != 0 + + # copy the cmakefile + status = os.system('rm -rf build') + assert os.WEXITSTATUS(status) == 0 + status = os.system('mkdir build ; cd build ; cp ../'+files[0]+' .') + assert os.WEXITSTATUS(status) == 0 + status = os.system('cd build ; cmake -DSOURCE_FILES=../source_sampling.cpp ; make') + assert os.WEXITSTATUS(status) == 0 + status = os.system('cp build/libsource.so .') + assert os.WEXITSTATUS(status) == 0 + status = os.system('rm -rf build') + assert os.WEXITSTATUS(status) == 0 + + return 0 + +class SourceTestHarness(PyAPITestHarness): + # build the test geometry + def _build_inputs(self): + mats = openmc.Materials() + + natural_lead = openmc.Material(1, "natural_lead") + natural_lead.add_element('Pb', 1,'ao') + natural_lead.set_density('g/cm3', 11.34) + mats.append(natural_lead) + + # surfaces + surface_sph1 = openmc.Sphere(r=100,boundary_type='vacuum') + volume_sph1 = -surface_sph1 + + # cell + cell_1 = openmc.Cell(region=volume_sph1) + cell_1.fill = natural_lead #assigning a material to a cell + universe = openmc.Universe(cells=[cell_1]) #hint, this list will need to include the new cell + geom = openmc.Geometry(universe) + + # settings + sett = openmc.Settings() + batches = 10 + sett.batches = batches + sett.inactive = 0 + sett.particles = 1000 + sett.particle = "neutron" + sett.run_mode = 'fixed source' + + #source + source = openmc.Source() + source.library = './libsource.so' + sett.source = source + + # run + model = openmc.model.Model(geom,mats,sett) + model.export_to_xml() + return + +def test_dlopen_source(): + compile_source() + harness = SourceTestHarness('statepoint.10.h5') + harness.main() diff --git a/tests/unit_tests/dlopen_source/test_dlopen_source.py b/tests/unit_tests/dlopen_source/test_dlopen_source.py deleted file mode 100644 index 917f3fd4fe..0000000000 --- a/tests/unit_tests/dlopen_source/test_dlopen_source.py +++ /dev/null @@ -1,57 +0,0 @@ -import openmc -import pytest -import os - -# compile the external source -def compile_source(): - # needs a more robust way to know where the - # Makefile is - status = os.system('cp /usr/local/share/Makefile .') - assert os.WEXITSTATUS(status) == 0 - status = os.system('make') - assert os.WEXITSTATUS(status) == 0 - - return - -# build the test geometry -def make_geometry(): - mats = openmc.Materials() - - natural_lead = openmc.Material(1, "natural_lead") - natural_lead.add_element('Pb', 1,'ao') - natural_lead.set_density('g/cm3', 11.34) - mats.append(natural_lead) - - # surfaces - surface_sph1 = openmc.Sphere(r=100) - volume_sph1 = -surface_sph1 - - # cell - cell_1 = openmc.Cell(region=volume_sph1) - cell_1.fill = natural_lead #assigning a material to a cell - universe = openmc.Universe(cells=[cell_1]) #hint, this list will need to include the new cell - geom = openmc.Geometry(universe) - - # settings - sett = openmc.Settings() - batches = 10 - sett.batches = batches - sett.inactive = 0 - sett.particles = 1000 - sett.particle = "neutron" - sett.run_mode = 'fixed source' - - #source - source = openmc.Source() - source.library = PATH+'source_sampling.so' - sett.source = source - - # run - model = openmc.model.Model(geom,mats,sett) - return model - -def test_dlopen_source(): - compile_source() - model = make_geometry() - model.run() - From 8784fc5d86723969ca232cc792bd21e5a66febf5 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Tue, 15 Oct 2019 22:01:43 +0100 Subject: [PATCH 41/59] Stylistic changes --- src/source.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index b7d1e24990..464a39948f 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -77,9 +77,9 @@ SourceDistribution::SourceDistribution(pugi::xml_node node) settings::path_source_library = get_node_value(node, "library", false, true); // check if it exists if (!file_exists(settings::path_source_library)) { - std::stringstream msg; - msg << "Library file " << settings::path_source_library << "' does not exist."; - fatal_error(msg); + std::stringstream msg; + msg << "Library file " << settings::path_source_library << "' does not exist."; + fatal_error(msg); } } else { From 8c432a2e5471cbb4d8d139ffa0b100e3629a2d5b Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Tue, 15 Oct 2019 22:10:33 +0100 Subject: [PATCH 42/59] added unit test for dlopen source, and changed library variable name --- openmc/source.py | 6 +++--- tests/unit_tests/test_source.py | 8 ++++++++ 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/openmc/source.py b/openmc/source.py index f8ea8ea936..339d41723e 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -49,7 +49,7 @@ class Source(object): self._angle = None self._energy = None self._file = None - self._source_library = None + self._library = None if space is not None: self.space = space @@ -70,7 +70,7 @@ class Source(object): @property def library(self): - return self._source_library + return self._library @property def space(self): @@ -100,7 +100,7 @@ class Source(object): @library.setter def library(self, library_name): cv.check_type('library', library_name, str) - self._source_library = library_name + self._library = library_name @space.setter def space(self, space): diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 1c70e159d2..d4d17a3dab 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -34,3 +34,11 @@ def test_source_file(): elem = src.to_xml_element() assert 'strength' in elem.attrib assert 'file' in elem.attrib + +def test_source_dlopen(): + library = './libsource.so' + src = openmc.Source(library=library) + assert src.library == library + + elem = src.to_xml_element() + assert 'library' in elem.attrib From 7d87979880b7143c7435b61a9f0006964828c55b Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Tue, 15 Oct 2019 22:58:36 +0100 Subject: [PATCH 43/59] Added supporting documentation --- docs/source/io_formats/settings.rst | 69 +++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 2e2648237b..1d35c18aaf 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -459,6 +459,15 @@ attributes/sub-elements: *Default*: None + :library: + If this attribute is given, it indicates that the source is to be instanciated + from an externally compiled source function. This source can be as complex as + is required to define the source for your problem. The only requirement that + is made upon this source, is that there is a function called sample_source() + more documentation on how to build sources can be found in :ref:`custom_source` + + *Deafult*: None + :space: An element specifying the spatial distribution of source sites. This element has the following attributes: @@ -591,6 +600,66 @@ attributes/sub-elements: *Default*: false +.. _custom_source: + +Custom Sources +++++++++++++++++++++++++++++++++++++ + +It is often the case that one may wish to simulate a complex source distribution, +which may include physics not present within OpenMC or to be phase space complex. It +is possible to define complex source with an externally defined source function +and loaded at runtime. A simple example source is shown below. + +.. code-block:: c++ + + #include + #include "openmc/random_lcg.h" + #include "openmc/source.h" + #include "openmc/particle.h" + + // you must have external C linkage here + extern "C" openmc::Particle::Bank sample_source() { + openmc::Particle::Bank particle; + // wgt + particle.particle = openmc::Particle::Type::neutron; + particle.wgt = 1.0; + // position + double angle = 2.*M_PI*openmc::prn(); + double radius = 3.0; + particle.r.x = radius*std::cos(angle); + particle.r.y = radius*std::sin(angle); + particle.r.z = 0.; + // angle + particle.u = {1.,0,0}; + particle.E = 14.08e6; + particle.delayed_group = 0; + return particle; + } + +The above source, creates 14.08 MeV neutrons, with an istropic direction +vector but distributed in a ring of radius three cm. This routine is +not particular complex, but should serve as an example upon which to build +more complicated sources. + + .. note:: The function signature must be declared to be extern. + + .. note:: You should only use the openmc::prn() random number generator + +In order to build your external source you need the CMakeLists.txt file that +was created during the OpenMC installation process (found in the share subdirectory) +and your source file(s). + +.. code-block:: bash + + cp $HOME/openmc/share/CMakeLists.txt. + mkdir bld + cd bld + cmake .. -DSOURCE_FILES=source_sampling.cpp + make + +You will now have a libsouce.so file in this directory, now point the library +attribute of source to this file and you will be able to sample particles. + .. _univariate: Univariate Probability Distributions From 94285c39ee29dfdafa49ed5c453f9f95e77126d5 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 16 Oct 2019 07:41:14 +0100 Subject: [PATCH 44/59] Wrapped the first instance trying to use dlopen in a posix safety ifdef --- src/source.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/source.cpp b/src/source.cpp index 464a39948f..066e3a52cc 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -2,7 +2,10 @@ #include // for move #include // for stringstream + +#if defined (__unix__) || (defined (__APPLE__) && defined (__MACH__)) #include // for dlopen +#endif #include #include "xtensor/xadapt.hpp" @@ -277,12 +280,17 @@ void initialize_source() msg << "Sampling from library source " << settings::path_source << "..."; write_message(msg, 6); + #if defined (__unix__) || (defined (__APPLE__) && defined (__MACH__)) // Open the library void* source_library = dlopen(settings::path_source_library.c_str(),RTLD_LAZY); if(!source_library) { std::stringstream msg("Couldnt open source library " + settings::path_source_library); fatal_error(msg); } + #else + std::stringstream msg("This feature has not yet been implemented for non POSIX systems"); + fatal_error(msg); + #endif // load the symbol typedef Particle::Bank (*sample_t)(); From 649a6336000cbe821007d97a7aa5ba155a7b63f3 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 16 Oct 2019 14:21:27 +0100 Subject: [PATCH 45/59] made the building of the cmake more sensible --- CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 18de941218..d61477f190 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -409,7 +409,7 @@ install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) #============================= # write the dl_open Makefile #============================= -file(WRITE ${CMAKE_INSTALL_PREFIX}/share/CMakeLists.txt " +file(WRITE ${CMAKE_BINARY_DIR}/CMakeLists.txt " cmake_minimum_required(VERSION 3.3 FATAL_ERROR) project(openmc_sources C CXX) add_library(source SHARED \$\{SOURCE_FILES\}) @@ -417,4 +417,4 @@ target_include_directories(source PUBLIC ${PROJECT_SOURCE_DIR}/include ${PROJECT_SOURCE_DIR}/vendor/pugixml ${HDF5_INCLUDE_DIRS}) target_link_libraries(source ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) ") - +install(FILES ${CMAKE_BINARY_DIR}/CMakeLists.txt DESTINATION share) From d81bf9a2817084dc96aa93d226db4261731a2c3d Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 16 Oct 2019 20:27:39 +0100 Subject: [PATCH 46/59] try adding an init py file --- tests/regression_tests/source_dlopen/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 tests/regression_tests/source_dlopen/__init__.py diff --git a/tests/regression_tests/source_dlopen/__init__.py b/tests/regression_tests/source_dlopen/__init__.py new file mode 100644 index 0000000000..e69de29bb2 From b39cbe404d610b7626f270bb5de0b5e1634c889a Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 16 Oct 2019 21:07:31 +0100 Subject: [PATCH 47/59] cxx standards for the source compile --- CMakeLists.txt | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d61477f190..ef67da4e50 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -412,7 +412,12 @@ install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) file(WRITE ${CMAKE_BINARY_DIR}/CMakeLists.txt " cmake_minimum_required(VERSION 3.3 FATAL_ERROR) project(openmc_sources C CXX) -add_library(source SHARED \$\{SOURCE_FILES\}) +add_library(source SHARED \$\{SOURCE_FILES\})") +get_target_property(cxx_std openmc CXX_STANDARD) +get_target_property(cxx_ext openmc CXX_EXTENSIONS) +file(APPEND ${CMAKE_BINARY_DIR}/CMakeLists.txt " +set_target_properties(source PROPERTIES CXX_STANDARD ${cxx_std} + CXX_EXTENSIONS ${cxx_ext}) target_include_directories(source PUBLIC ${PROJECT_SOURCE_DIR}/include ${PROJECT_SOURCE_DIR}/vendor/pugixml ${HDF5_INCLUDE_DIRS}) target_link_libraries(source ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) From f1d8565ef4909013e2c1ffda282d3035852ee7e0 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Tue, 11 Feb 2020 22:07:53 +0000 Subject: [PATCH 48/59] updated signatures according to new changes in how seed is set --- src/source.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 066e3a52cc..b1fd8c71e9 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -293,7 +293,7 @@ void initialize_source() #endif // load the symbol - typedef Particle::Bank (*sample_t)(); + typedef Particle::Bank (*sample_t)(uint64_t *seed); // reset errors dlerror(); @@ -315,10 +315,10 @@ void initialize_source() // initialize random number seed int64_t id = simulation::total_gen*settings::n_particles + simulation::work_index[mpi::rank] + i + 1; - set_particle_seed(id); + uint64_t seed = init_seed(id, STREAM_SOURCE); // sample external source distribution - simulation::source_bank[i] = sample_source(); + simulation::source_bank[i] = sample_source(&seed); } // release the library dlclose(source_library); @@ -405,7 +405,7 @@ void fill_source_bank_fixedsource() } // load the symbol - typedef Particle::Bank (*sample_t)(); + typedef Particle::Bank (*sample_t)(uint64_t *seed); // reset errors dlerror(); @@ -427,10 +427,10 @@ void fill_source_bank_fixedsource() // initialize random number seed int64_t id = simulation::total_gen*settings::n_particles + simulation::work_index[mpi::rank] + i + 1; - set_particle_seed(id); - + uint64_t seed = init_seed(id, STREAM_SOURCE); + // sample external source distribution - simulation::source_bank[i] = sample_source(); + simulation::source_bank[i] = sample_source(&seed); } // release the library From e98d7fb28e1aea224c824342922c60cb9d363c17 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 09:29:02 +0000 Subject: [PATCH 49/59] Updated the source routines according to passing seed to prn --- src/source.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index b1fd8c71e9..465c2aac14 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -293,7 +293,7 @@ void initialize_source() #endif // load the symbol - typedef Particle::Bank (*sample_t)(uint64_t *seed); + typedef Particle::Bank (*sample_t)(uint64_t seed); // reset errors dlerror(); @@ -318,7 +318,7 @@ void initialize_source() uint64_t seed = init_seed(id, STREAM_SOURCE); // sample external source distribution - simulation::source_bank[i] = sample_source(&seed); + simulation::source_bank[i] = sample_source(seed); } // release the library dlclose(source_library); @@ -405,7 +405,7 @@ void fill_source_bank_fixedsource() } // load the symbol - typedef Particle::Bank (*sample_t)(uint64_t *seed); + typedef Particle::Bank (*sample_t)(uint64_t seed); // reset errors dlerror(); @@ -430,7 +430,7 @@ void fill_source_bank_fixedsource() uint64_t seed = init_seed(id, STREAM_SOURCE); // sample external source distribution - simulation::source_bank[i] = sample_source(&seed); + simulation::source_bank[i] = sample_source(seed); } // release the library From 2823dd07e9db9397c657da2edf4b9e496b5cc3d3 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 09:43:35 +0000 Subject: [PATCH 50/59] Updated documentation to reflect changes --- docs/source/io_formats/settings.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 1d35c18aaf..fde374b8d1 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -618,7 +618,7 @@ and loaded at runtime. A simple example source is shown below. #include "openmc/particle.h" // you must have external C linkage here - extern "C" openmc::Particle::Bank sample_source() { + extern "C" openmc::Particle::Bank sample_source(const int64_t seed) { openmc::Particle::Bank particle; // wgt particle.particle = openmc::Particle::Type::neutron; From 8f4b8535a6e9e24a686896b76d3ec9c1eab764c5 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 12:24:29 +0000 Subject: [PATCH 51/59] Added statements for additional cmake targets --- CMakeLists.txt | 20 +------------------- 1 file changed, 1 insertion(+), 19 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ef67da4e50..0d4c371f02 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -207,7 +207,6 @@ endif() add_library(faddeeva STATIC vendor/faddeeva/Faddeeva.cc) target_include_directories(faddeeva PUBLIC - $ $ ) target_compile_options(faddeeva PRIVATE ${cxxflags}) @@ -389,7 +388,7 @@ add_custom_command(TARGET libopenmc POST_BUILD #=============================================================================== set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC) -install(TARGETS openmc libopenmc faddeeva +install(TARGETS openmc libopenmc faddeeva pugixml gsl-lite EXPORT openmc-targets RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} @@ -406,20 +405,3 @@ install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1) install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright) install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) -#============================= -# write the dl_open Makefile -#============================= -file(WRITE ${CMAKE_BINARY_DIR}/CMakeLists.txt " -cmake_minimum_required(VERSION 3.3 FATAL_ERROR) -project(openmc_sources C CXX) -add_library(source SHARED \$\{SOURCE_FILES\})") -get_target_property(cxx_std openmc CXX_STANDARD) -get_target_property(cxx_ext openmc CXX_EXTENSIONS) -file(APPEND ${CMAKE_BINARY_DIR}/CMakeLists.txt " -set_target_properties(source PROPERTIES CXX_STANDARD ${cxx_std} - CXX_EXTENSIONS ${cxx_ext}) -target_include_directories(source PUBLIC ${PROJECT_SOURCE_DIR}/include - ${PROJECT_SOURCE_DIR}/vendor/pugixml ${HDF5_INCLUDE_DIRS}) -target_link_libraries(source ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES}) -") -install(FILES ${CMAKE_BINARY_DIR}/CMakeLists.txt DESTINATION share) From 59b00b2171281a39a21e9a1fe534b0376e4e8daf Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 12:44:46 +0000 Subject: [PATCH 52/59] Updated spelling mistakes --- src/source.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/source.cpp b/src/source.cpp index 465c2aac14..f910b384c9 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -284,7 +284,7 @@ void initialize_source() // Open the library void* source_library = dlopen(settings::path_source_library.c_str(),RTLD_LAZY); if(!source_library) { - std::stringstream msg("Couldnt open source library " + settings::path_source_library); + std::stringstream msg("Couldn't open source library " + settings::path_source_library); fatal_error(msg); } #else @@ -306,7 +306,7 @@ void initialize_source() if (dlsym_error) { std::cout << dlsym_error << std::endl; dlclose(source_library); - fatal_error("Couldnt open the sample_source symbol"); + fatal_error("Couldn't open the sample_source symbol"); } // Generation source sites from specified distribution in the @@ -400,7 +400,7 @@ void fill_source_bank_fixedsource() // Open the library void* source_library = dlopen(settings::path_source_library.c_str(),RTLD_LAZY); if(!source_library) { - std::stringstream msg("Couldnt open source library " + settings::path_source_library); + std::stringstream msg("Couldn't open source library " + settings::path_source_library); fatal_error(msg); } @@ -418,7 +418,7 @@ void fill_source_bank_fixedsource() if (dlsym_error) { std::cout << dlsym_error << std::endl; dlclose(source_library); - fatal_error("Couldnt open the sample_source symbol"); + fatal_error("Couldn't open the sample_source symbol"); } // Generation source sites from specified distribution in the From 13133b87e535319c2986427d2969042e77c30efd Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 12:47:39 +0000 Subject: [PATCH 53/59] Corrected style according to review --- examples/xml/custom_source/source_ring.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/xml/custom_source/source_ring.cpp b/examples/xml/custom_source/source_ring.cpp index ef2784ea36..fb361d1aa9 100644 --- a/examples/xml/custom_source/source_ring.cpp +++ b/examples/xml/custom_source/source_ring.cpp @@ -12,14 +12,14 @@ extern "C" openmc::Particle::Bank sample_source() { particle.wgt = 1.0; // position - double angle = 2.*M_PI*openmc::prn(); + double angle = 2. * M_PI * openmc::prn(); double radius = 3.0; - particle.r.x = radius*std::cos(angle); - particle.r.y = radius*std::sin(angle); - particle.r.z = 0.; + particle.r.x = radius * std::cos(angle); + particle.r.y = radius * std::sin(angle); + particle.r.z = 0.0; // angle - particle.u = {1.,0,0}; + particle.u = {1.0, 0.0, 0.0}; particle.E = 14.08e6; particle.delayed_group = 0; return particle; -} \ No newline at end of file +} From a392353a36fb83547647dbb76eeaa0794c6d8f8b Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 12:48:54 +0000 Subject: [PATCH 54/59] line removal from review --- src/source.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/source.cpp b/src/source.cpp index f910b384c9..e7a476746e 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -275,7 +275,6 @@ void initialize_source() file_close(file_id); } else if ( settings::path_source_library != "" ) { // Get the source from a library object - std::stringstream msg; msg << "Sampling from library source " << settings::path_source << "..."; write_message(msg, 6); From 1fd6e003e803c10febdde8dde58aaf84f1ee7665 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 12:50:53 +0000 Subject: [PATCH 55/59] further changes from review --- docs/source/io_formats/settings.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index fde374b8d1..f224df265f 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -658,7 +658,8 @@ and your source file(s). make You will now have a libsouce.so file in this directory, now point the library -attribute of source to this file and you will be able to sample particles. +attribute of the source XML element to this file and you will be able to sample +particles. .. _univariate: From 728969c5f2c43e429be5bd0ba55cefa64ed85e86 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 12:55:21 +0000 Subject: [PATCH 56/59] further review comments regarding style --- docs/source/io_formats/settings.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index f224df265f..5392efef5f 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -603,12 +603,12 @@ attributes/sub-elements: .. _custom_source: Custom Sources -++++++++++++++++++++++++++++++++++++ +++++++++++++++ It is often the case that one may wish to simulate a complex source distribution, which may include physics not present within OpenMC or to be phase space complex. It -is possible to define complex source with an externally defined source function -and loaded at runtime. A simple example source is shown below. +is possible to define a complex source with an externally defined source function +that is loaded at runtime. A simple example source is shown below. .. code-block:: c++ @@ -637,7 +637,7 @@ and loaded at runtime. A simple example source is shown below. } The above source, creates 14.08 MeV neutrons, with an istropic direction -vector but distributed in a ring of radius three cm. This routine is +vector but distributed in a ring with a 3 cm radius. This routine is not particular complex, but should serve as an example upon which to build more complicated sources. From 7c2b695c982ff2f312789dc9557de080dced4e9e Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 13:06:30 +0000 Subject: [PATCH 57/59] Refactor into function as suggested by review --- include/openmc/source.h | 3 ++ src/source.cpp | 110 +++++++++++++++++----------------------- 2 files changed, 49 insertions(+), 64 deletions(-) diff --git a/include/openmc/source.h b/include/openmc/source.h index 0a73260fd1..e980f2a8c8 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -72,6 +72,9 @@ Particle::Bank sample_external_source(uint64_t* seed); //! Fill source bank at end of generation for fixed source simulations void fill_source_bank_fixedsource(); +//! Fill source bank at the end of a generation for dlopen based source simulation +void fill_source_bank_dlopen_source(); + void free_memory_source(); } // namespace openmc diff --git a/src/source.cpp b/src/source.cpp index e7a476746e..2cda76a5c0 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -297,30 +297,7 @@ void initialize_source() // reset errors dlerror(); - // get the function from the library - sample_t sample_source = (sample_t) dlsym(source_library, "sample_source"); - const char *dlsym_error = dlerror(); - - // check for any dlsym errors - if (dlsym_error) { - std::cout << dlsym_error << std::endl; - dlclose(source_library); - fatal_error("Couldn't open the sample_source symbol"); - } - - // Generation source sites from specified distribution in the - // library source - for (int64_t i = 0; i < simulation::work_per_rank; ++i) { - // initialize random number seed - int64_t id = simulation::total_gen*settings::n_particles + - simulation::work_index[mpi::rank] + i + 1; - uint64_t seed = init_seed(id, STREAM_SOURCE); - - // sample external source distribution - simulation::source_bank[i] = sample_source(seed); - } - // release the library - dlclose(source_library); + fill_source_bank_dlopen_source(); } else { // Generation source sites from specified distribution in user input @@ -381,6 +358,50 @@ void free_memory_source() model::external_sources.clear(); } +// fill the source bank from the external source +void fill_source_bank_dlopen_source() +{ + std::stringstream msg; + + // Open the library + void* source_library = dlopen(settings::path_source_library.c_str(),RTLD_LAZY); + if(!source_library) { + std::stringstream msg("Couldn't open source library " + settings::path_source_library); + fatal_error(msg); + } + + // load the symbol + typedef Particle::Bank (*sample_t)(uint64_t seed); + + // reset errors + dlerror(); + + // get the function from the library + sample_t sample_source = (sample_t) dlsym(source_library, "sample_source"); + const char *dlsym_error = dlerror(); + + // check for any dlsym errors + if (dlsym_error) { + std::cout << dlsym_error << std::endl; + dlclose(source_library); + fatal_error("Couldn't open the sample_source symbol"); + } + + // Generation source sites from specified distribution in the + // library source + for (int64_t i = 0; i < simulation::work_per_rank; ++i) { + // initialize random number seed + int64_t id = simulation::total_gen*settings::n_particles + + simulation::work_index[mpi::rank] + i + 1; + uint64_t seed = init_seed(id, STREAM_SOURCE); + // sample external source distribution + simulation::source_bank[i] = sample_source(seed); + } + + // release the library + dlclose(source_library); +} + void fill_source_bank_fixedsource() { if (settings::path_source.empty() && settings::path_source_library.empty()) { @@ -394,46 +415,7 @@ void fill_source_bank_fixedsource() simulation::source_bank[i] = sample_external_source(&seed); } } else if (settings::path_source.empty() && !settings::path_source.empty()) { - std::stringstream msg; - - // Open the library - void* source_library = dlopen(settings::path_source_library.c_str(),RTLD_LAZY); - if(!source_library) { - std::stringstream msg("Couldn't open source library " + settings::path_source_library); - fatal_error(msg); - } - - // load the symbol - typedef Particle::Bank (*sample_t)(uint64_t seed); - - // reset errors - dlerror(); - - // get the function from the library - sample_t sample_source = (sample_t) dlsym(source_library, "sample_source"); - const char *dlsym_error = dlerror(); - - // check for any dlsym errors - if (dlsym_error) { - std::cout << dlsym_error << std::endl; - dlclose(source_library); - fatal_error("Couldn't open the sample_source symbol"); - } - - // Generation source sites from specified distribution in the - // library source - for (int64_t i = 0; i < simulation::work_per_rank; ++i) { - // initialize random number seed - int64_t id = simulation::total_gen*settings::n_particles + - simulation::work_index[mpi::rank] + i + 1; - uint64_t seed = init_seed(id, STREAM_SOURCE); - - // sample external source distribution - simulation::source_bank[i] = sample_source(seed); - } - - // release the library - dlclose(source_library); + fill_source_bank_dlopen_source(); } } From 09581ebf3816b234127dc2b4bee4b8004f8666e9 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 13:11:22 +0000 Subject: [PATCH 58/59] updated for new cmake instructions --- docs/source/io_formats/settings.rst | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 5392efef5f..6f898f911b 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -645,19 +645,17 @@ more complicated sources. .. note:: You should only use the openmc::prn() random number generator -In order to build your external source you need the CMakeLists.txt file that -was created during the OpenMC installation process (found in the share subdirectory) -and your source file(s). +In order to build your external source you need the following CMakeLists.txt file -.. code-block:: bash +.. code-block:: cmake - cp $HOME/openmc/share/CMakeLists.txt. - mkdir bld - cd bld - cmake .. -DSOURCE_FILES=source_sampling.cpp - make + cmake_minimum_required(VERSION 3.3 FATAL_ERROR) + project(openmc_sources CXX) + add_library(source SHARED source_ring.cpp) + find_package(OpenMC REQUIRED HINTS ) + target_link_libraries(source OpenMC::libopenmc) -You will now have a libsouce.so file in this directory, now point the library +You will now have a libsouce.so (or .dylib) file in this directory, now point the library attribute of the source XML element to this file and you will be able to sample particles. From fba8f6a8e0a932d90a1c09b9a02d8f60fcaadbd7 Mon Sep 17 00:00:00 2001 From: Andrew Davis Date: Wed, 12 Feb 2020 13:11:53 +0000 Subject: [PATCH 59/59] updated the example source according to style --- examples/xml/custom_source/source_ring.cpp | 33 +++++++++++----------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/examples/xml/custom_source/source_ring.cpp b/examples/xml/custom_source/source_ring.cpp index fb361d1aa9..babb2f591c 100644 --- a/examples/xml/custom_source/source_ring.cpp +++ b/examples/xml/custom_source/source_ring.cpp @@ -5,21 +5,22 @@ // you must have external C linkage here otherwise // dlopen will not find the file -extern "C" openmc::Particle::Bank sample_source() { - openmc::Particle::Bank particle; - // wgt - particle.particle = openmc::Particle::Type::neutron; - particle.wgt = 1.0; - // position +extern "C" openmc::Particle::Bank sample_source() +{ + openmc::Particle::Bank particle; + // wgt + particle.particle = openmc::Particle::Type::neutron; + particle.wgt = 1.0; + // position - double angle = 2. * M_PI * openmc::prn(); - double radius = 3.0; - particle.r.x = radius * std::cos(angle); - particle.r.y = radius * std::sin(angle); - particle.r.z = 0.0; - // angle - particle.u = {1.0, 0.0, 0.0}; - particle.E = 14.08e6; - particle.delayed_group = 0; - return particle; + double angle = 2. * M_PI * openmc::prn(); + double radius = 3.0; + particle.r.x = radius * std::cos(angle); + particle.r.y = radius * std::sin(angle); + particle.r.z = 0.0; + // angle + particle.u = {1.0, 0.0, 0.0}; + particle.E = 14.08e6; + particle.delayed_group = 0; + return particle; }