diff --git a/openmc/model/model.py b/openmc/model/model.py index 6e4c1c5856d..e02ff4f887e 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -16,6 +16,8 @@ import lxml.etree as ET import numpy as np from scipy.optimize import curve_fit +from scipy.stats import gmean +from uncertainties import ufloat, UFloat import openmc import openmc._xml as xml @@ -2361,6 +2363,195 @@ def _replace_infinity(value): # Take a wild guess as to how many rays are needed self.settings.particles = 2 * int(max_length) + def add_derivative_tallies( + self, + deriv_variable: str, + deriv_material: int, + deriv_nuclide: str | None = None + ) -> None: + """Add base and derivative tallies to model for keff derivative extraction. + + This is a convenience method that adds the four tallies required for + extracting dk/dx during keff_search: base nu-fission, base absorption, + derivative nu-fission, and derivative absorption. + + Note: This method is automatically called by keff_search when + use_derivative_tallies=True. You only need to call it manually if + you want to set up tallies before calling keff_search, or if you're + using the tallies for other purposes. + + Parameters + ---------- + deriv_variable : str + Type of derivative: 'density', 'nuclide_density', or 'temperature'. + + .. warning:: + **Temperature derivatives have severe limitations:** + + - Require Windowed Multipole (WMP) cross section data + - Only valid in resolved resonance range (~1 eV to ~10 keV) + - Not available for most nuclides in standard data libraries + - Not compatible with cross section interpolation + + Temperature derivatives are **not recommended** for practical + k-eff searches. Use derivative-free search instead. + See `src/tallies/derivative.cpp` lines 283-500 for implementation details. + + deriv_material : int + Material ID to perturb + deriv_nuclide : str, optional + Nuclide name for nuclide_density derivatives (e.g., 'B10', 'U235'). + Required if deriv_variable='nuclide_density'. + + Examples + -------- + >>> model = openmc.examples.pwr_pin_cell() + >>> # Add tallies for boron concentration derivative + >>> model.add_derivative_tallies('nuclide_density', 3, 'B10') + >>> # Add tallies for fuel density derivative + >>> model.add_derivative_tallies('density', 1) + + """ + # Base tallies + t_fission = openmc.Tally(name='base_fission') + t_fission.scores = ['nu-fission'] + + t_absorption = openmc.Tally(name='base_absorption') + t_absorption.scores = ['absorption'] + + tallies = [t_fission, t_absorption] + + # Derivative tallies + deriv = openmc.TallyDerivative( + variable=deriv_variable, + material=deriv_material, + nuclide=deriv_nuclide + ) + + t_fission_deriv = openmc.Tally(name=f'fission_deriv_{deriv_variable}') + t_fission_deriv.scores = ['nu-fission'] + t_fission_deriv.derivative = deriv + + t_absorption_deriv = openmc.Tally(name=f'absorption_deriv_{deriv_variable}') + t_absorption_deriv.scores = ['absorption'] + t_absorption_deriv.derivative = deriv + + tallies.extend([t_fission_deriv, t_absorption_deriv]) + self.tallies = openmc.Tallies(tallies) + + def _extract_derivative_constraint( + self, + sp: 'openmc.StatePoint', + deriv_variable: str, + deriv_material: int, + deriv_nuclide: str | None = None, + deriv_to_x_func: Callable[[UFloat], UFloat] | None = None, + ) -> UFloat | None: + r"""Extract dk_eff/dx from StatePoint using derivative tallies. + + This method implements a generic approach to compute the derivative of + k-effective with respect to any perturbation variable (density, + nuclide_density, temperature) by using base and derivative + tallies and the quotient rule: + + .. math:: + \frac{dk}{dx} = \frac{A \frac{dF}{dx} - F \frac{dA}{dx}}{A^2} + + where :math:`F` is the fission production tally (nu-fission score), + :math:`A` is the absorption tally, and :math:`\frac{dF}{dx}`, + :math:`\frac{dA}{dx}` are their derivative counterparts. + + For **nuclide_density derivatives**, OpenMC's C++ backend computes derivatives + with respect to **number density N** (atoms/cm³). If the search parameter x + is a different quantity (e.g., mass parts-per-million for boron), the + caller must provide ``deriv_to_x_func`` to convert: dk/dx = (dk/dN) × (dN/dx). + + Uncertainties are propagated using linear error propagation. + + Parameters + ---------- + sp : openmc.StatePoint + StatePoint after MC run + deriv_variable : str + Type of derivative: 'density', 'nuclide_density', or 'temperature' + deriv_material : int + Material ID being perturbed + deriv_nuclide : str, optional + Nuclide name (required for 'nuclide_density', ignored otherwise) + deriv_to_x_func : callable, optional + For nuclide_density: function to convert dN/dx. Signature: + ``deriv_to_x_func(deriv_value_dN) -> float`` returns dk/dx given dk/dN. + If not provided, returns dk/dN unchanged. + Ignored for other derivative types. + + Returns + ------- + UFloat | None + dk/dx if base and derivative tallies found, else None. + For nuclide_density without deriv_to_x_func, + returned derivative is dk/dN (where N is number density in atoms/cm³). + """ + try: + # Find base tallies (nu-fission and absorption) + base_fission = None + base_absorption = None + deriv_fission = None + deriv_absorption = None + + for tally_id, tally in sp.tallies.items(): + scores = getattr(tally, 'scores', []) or [] + if tally.derivative is None: + # Base tallies (no derivative) + if 'nu-fission' in scores and base_fission is None: + base_fission = tally + if 'absorption' in scores and base_absorption is None: + base_absorption = tally + else: + # Derivative tallies: check if they match the requested variable + deriv = tally.derivative + if (deriv.variable == deriv_variable and + deriv.material == deriv_material and + (deriv_variable != 'nuclide_density' or + deriv.nuclide == deriv_nuclide)): + if 'nu-fission' in scores and deriv_fission is None: + deriv_fission = tally + if 'absorption' in scores and deriv_absorption is None: + deriv_absorption = tally + check_type('base_fission', base_fission, openmc.Tally) + check_type('base_absorption', base_absorption, openmc.Tally) + check_type('deriv_fission', deriv_fission, openmc.Tally) + check_type('deriv_absorption', deriv_absorption, openmc.Tally) + + def tally_to_ufloat(t): + return ufloat(t.mean.squeeze(), + t.std_dev.squeeze()) + + F = tally_to_ufloat(base_fission) + A = tally_to_ufloat(base_absorption) + dF_dx = tally_to_ufloat(deriv_fission) + dA_dx = tally_to_ufloat(deriv_absorption) + + print(f' [DERIV-EXTRACT] Found all 4 tallies for {deriv_variable}') + print(f' [DERIV-EXTRACT] F={F.n:.6e}, A={A.n:.6e}, dF/dx={dF_dx.n:.6e}, dA/dx={dA_dx.n:.6e}') + + # Quotient rule: dk/dx = (A * dF/dx - F * dA/dx) / A^2 + dk_dx = (A * dF_dx - F * dA_dx) / (A * A) + dk_dx_before = dk_dx + print(f' [DERIV-EXTRACT] Computed dk/dx = {dk_dx:.6e} (before any conversion)') + + # For nuclide_density: convert dk/dN to dk/dx if conversion provided + if deriv_variable == 'nuclide_density' and deriv_to_x_func is not None: + # deriv_to_x_func converts one derivative value + # It should return the scaled derivative (dk/dx = (dk/dN) * (dN/dx)) + dk_dx = deriv_to_x_func(dk_dx) + print(f' [DERIV-EXTRACT] Applied deriv_to_x_func: dk/dN={dk_dx_before:.6e} -> dk/dx={dk_dx:.6e}') + return dk_dx + + except TypeError as e: + warnings.warn(f" [DERIV-EXTRACT] ERROR: Could not extract derivative: {e}") + + return None + def keff_search( self, func: ModelModifier, @@ -2379,6 +2570,11 @@ def keff_search( b_max: int | None = None, maxiter: int = 50, output: bool = False, + use_derivative_tallies: bool = False, + deriv_variable: str | None = None, + deriv_material: int | None = None, + deriv_nuclide: str | None = None, + deriv_to_x_func: Callable[[UFloat], UFloat] | None = None, func_kwargs: dict[str, Any] | None = None, run_kwargs: dict[str, Any] | None = None, ) -> SearchResult: @@ -2392,6 +2588,11 @@ def keff_search( point to evaluate. It also adaptively changes the number of batches to meet the target uncertainty value at each iteration. + When derivative tallies are enabled, the gradient-based least-squares + approach is inspired by the methodology developed by Sterling Harper + in `"Calculating Reaction Rate Derivatives in Monte Carlo Neutron + Transport" `_. + The target uncertainty for iteration :math:`n+1` is determined by the following equation (following Eq. (8) in the paper): @@ -2441,6 +2642,50 @@ def keff_search( Maximum number of iterations to perform. output : bool, optional Whether or not to display output showing iteration progress. + use_derivative_tallies : bool, optional + If True, extract derivative tallies from StatePoints and use them + as gradient constraints in the curve fitting process. Requires + deriv_variable and deriv_material to be specified. Default is False. + Derivative tallies are automatically added to the model. + deriv_variable : str, optional + Type of derivative to extract. Supported values: 'density', + 'nuclide_density', 'temperature'. Required if + use_derivative_tallies=True. Example: 'nuclide_density' to + perturb a specific nuclide concentration; 'density' for material + mass density; 'temperature' for Doppler feedback. + + Note: Could be inferred as 'nuclide_density' when deriv_nuclide + is provided, but explicit specification is clearer and allows + for density/temperature derivatives without nuclide ambiguity. + deriv_material : int, optional + Material ID to perturb for derivatives. Required if + use_derivative_tallies=True. Example: Material ID 3 for boron + in coolant, Material ID 1 for fuel. + + Note: Could potentially be inferred by searching all materials + for the specified nuclide, but this would be ambiguous when + multiple materials contain the same element (e.g., oxygen in + both fuel and coolant). Explicit specification avoids ambiguity. + deriv_nuclide : str, optional + Nuclide name (e.g., 'B10', 'U235') for nuclide_density derivatives. + Ignored for other derivative types. Required if + deriv_variable='nuclide_density'. + deriv_to_x_func : callable, optional + For nuclide_density derivatives: a function that computes the conversion + dN/dx where N is the nuclide number density and x is the search parameter. + Signature: ``deriv_to_x_func(deriv_value) -> float`` + + If not provided, returns dk/dN (not dk/dx) for nuclide_density. + Ignored for other derivative types. + + NOTE: When derivative tallies are enabled, derivatives are automatically + used as gradient constraints in the least-squares fitting with their + uncertainties as weights. Derivatives are also normalized by their + magnitude (geometric mean of absolute values) to ensure numerical + stability, handling large magnitudes (e.g., dk/dppm ∼ 10^20) without + requiring manual scaling. Temperature derivatives are not yet supported. + For k-eff searches involving temperature changes, use the standard + derivative-free search (set `use_derivative_tallies=False`). func_kwargs : dict, optional Keyword-based arguments to pass to the `func` function. run_kwargs : dict, optional @@ -2461,6 +2706,20 @@ def keff_search( check_type('target', target, Real) if memory < 2: raise ValueError("memory must be ≥ 2") + + # Validate derivative parameters + if use_derivative_tallies: + # Validate against C++ backend supported types (see src/tallies/derivative.cpp) + check_value('deriv_variable', + deriv_variable, + ('density', 'nuclide_density', 'temperature')) + check_type('deriv_material', deriv_material, int) + if deriv_variable == 'nuclide_density': + check_type('deriv_nuclide', deriv_nuclide, str) + + # Automatically add derivative tallies to the model + self.add_derivative_tallies(deriv_variable, deriv_material, deriv_nuclide) + func_kwargs = {} if func_kwargs is None else dict(func_kwargs) run_kwargs = {} if run_kwargs is None else dict(run_kwargs) run_kwargs.setdefault('output', False) @@ -2470,10 +2729,12 @@ def keff_search( fs: list[float] = [] ss: list[float] = [] gs: list[int] = [] + dks: list[float] = [] # dk/dx derivatives + dks_std: list[float] = [] # uncertainties in derivatives count = 0 # Helper function to evaluate f and store results - def eval_at(x: float, batches: int) -> tuple[float, float]: + def eval_at(x: float, batches: int) -> tuple[float, float, float | None, float | None]: # Modify the model with the current guess func(x, **func_kwargs) @@ -2491,17 +2752,37 @@ def eval_at(x: float, batches: int) -> tuple[float, float]: # Extract keff and its uncertainty with openmc.StatePoint(sp_filepath) as sp: keff = sp.keff + dk_dx = None + + # If requested, extract derivative constraint using generic method + if use_derivative_tallies and deriv_variable and deriv_material: + dk_dx = self._extract_derivative_constraint( + sp, deriv_variable, deriv_material, deriv_nuclide, + deriv_to_x_func + ) + if output and dk_dx is not None: + print(f' [DERIV] Extracted dk/dx={dk_dx:.6e}') if output: nonlocal count count += 1 - print(f'Iteration {count}: {batches=}, {x=:.6g}, {keff=:.5f}') + msg = f'Iteration {count}: {batches=}, {x=:.6g}, {keff=:.5f}' + if dk_dx is not None: + msg = f'{msg}, dk/dx={dk_dx.n:.6g}' + print(msg) xs.append(float(x)) fs.append(float(keff.n - target)) ss.append(float(keff.s)) gs.append(int(batches)) - return fs[-1], ss[-1] + if dk_dx is not None: + dks.append(dk_dx.n) + dks_std.append(dk_dx.s) + return fs[-1], ss[-1], dk_dx.n, dk_dx.s + else: + dks.append(0.0) + dks_std.append(0.0) + return fs[-1], ss[-1], None, None # Default b0 to current model settings if not explicitly provided if b0 is None: @@ -2513,24 +2794,108 @@ def eval_at(x: float, batches: int) -> tuple[float, float]: run_kwargs.setdefault('cwd', tmpdir) # ---- Seed with two evaluations - f0, s0 = eval_at(x0, b0) + f0, s0, dk0, dks0 = eval_at(x0, b0) if abs(f0) <= k_tol and s0 <= sigma_final: return SearchResult(x0, xs, fs, ss, gs, True, "converged") - f1, s1 = eval_at(x1, b0) + f1, s1, dk1, dks1 = eval_at(x1, b0) if abs(f1) <= k_tol and s1 <= sigma_final: return SearchResult(x1, xs, fs, ss, gs, True, "converged") for _ in range(maxiter - 2): - # ------ Step 1: propose next x via GRsecant + # ------ Step 1: propose next x + + # Perform weighted least squares fit on f(x) = a + bx + # If derivative tallies enabled: augment with gradient constraints m = min(memory, len(xs)) - # Perform a curve fit on f(x) = a + bx accounting for - # uncertainties. This is equivalent to minimizing the function - # in Equation (A.14) - (a, b), _ = curve_fit( - lambda x, a, b: a + b*x, - xs[-m:], fs[-m:], sigma=ss[-m:], absolute_sigma=True - ) + xs_fit = np.array(xs[-m:]) + fs_fit = np.array(fs[-m:]) + ss_fit = np.array(ss[-m:]) + + def custom_curve_fit(): + # Perform a curve fit on f(x) = a + bx accounting for uncertainties + # If derivatives are available, augment with gradient constraints + if use_derivative_tallies and any(dks[-m:]): + # Gradient-augmented least squares fit + # Minimize: sum_i (f_i - a - b*x_i)^2 / sigma_i^2 + # + sum_j (b - dk_j/dx_j)^2 / (dk_std_j)^2 + + dks_fit = np.array(dks[-m:]) + dks_std_fit = np.array(dks_std[-m:]) + + # Build augmented system: minimize both point residuals and gradient errors + # Points with valid derivatives contribute dual constraints + valid_derivs = dks_std_fit > 0 + n_pts = len(xs_fit) + n_derivs = np.sum(valid_derivs) + + # Construct augmented system matrix + A = np.vstack([ + np.ones(n_pts) / ss_fit, + xs_fit / ss_fit, + ]).T + b_vec = fs_fit / ss_fit + + # Add gradient constraints (b should match dk/dx at each point) + if n_derivs > 0: + # Gradient constraints: f(x) = a + bx, so df/dx = b + # Constraint: b ≈ dk/dx_j, weighted by 1/dk_std_j + + # AUTO-NORMALIZE DERIVATIVES: When derivatives are very large or very small, + # normalize by their magnitude to avoid ill-conditioned least squares system. + # This is critical for derivatives like dk/dppm which can be O(10^20). + valid_deriv_values = dks_fit[valid_derivs] + valid_deriv_stds = dks_std_fit[valid_derivs] + + if output: + print(f' [DERIV-FIT] Using {n_derivs} derivative constraints in curve fit') + print(f' [DERIV-FIT] Raw derivatives: {valid_deriv_values}') + + # Calculate normalization scale: geometric mean of absolute derivative magnitudes + abs_derivs = np.abs(valid_deriv_values) + abs_derivs = abs_derivs[abs_derivs > 0] # Exclude zeros + deriv_scale = gmean(abs_derivs) if len(abs_derivs) > 0 else 1.0 + + if output: + print(f' [DERIV-FIT] Normalization scale factor: {deriv_scale:.6e}') + + # Apply scaling to derivatives and their uncertainties + scaled_derivs = valid_deriv_values / deriv_scale + scaled_deriv_stds = valid_deriv_stds / deriv_scale + + # Build constraint rows with normalized derivatives + deriv_rows = np.zeros((n_derivs, 2)) + deriv_rows[:, 0] = 0.0 # a coefficient in gradient = 0 + deriv_rows[:, 1] = 1.0 # b coefficient + + # Normalized targets: scale-invariant constraint weighted by uncertainty + deriv_targets = scaled_derivs / scaled_deriv_stds + + A = np.vstack([A, deriv_rows]) + b_vec = np.hstack([b_vec, deriv_targets]) + + # Solve least squares: (A^T A)^{-1} A^T b + try: + coeffs, residuals, rank, s = np.linalg.lstsq(A, b_vec, rcond=None) + a, b = float(coeffs[0]), float(coeffs[1]) + if output: + print(f' [DERIV-FIT] Fitted line (with derivative constraints): f(x) = {a:.6e} + {b:.6e}*x') + return a, b + except np.linalg.LinAlgError: + pass + + # Perform a curve fit on f(x) = a + bx accounting for + # uncertainties. This is equivalent to minimizing the function + # in Equation (A.14) + (a, b), _ = curve_fit( + lambda x, a, b: a + b*x, + xs_fit, fs_fit, sigma=ss_fit, absolute_sigma=True + ) + if output: + print(f' [NO-DERIV-FIT] Standard fit: f(x) = {a:.6e} + {b:.6e}*x (no derivatives)') + return a, b + + a, b = custom_curve_fit() x_new = float(-a / b) # Clamp x_new to the bounds if provided @@ -2570,7 +2935,7 @@ def eval_at(x: float, batches: int) -> tuple[float, float]: b_new = min(b_new, b_max) # Evaluate at proposed x with batches determined above - f_new, s_new = eval_at(x_new, b_new) + f_new, s_new, _, _ = eval_at(x_new, b_new) # Termination based on both criteria (|f| and σ) if abs(f_new) <= k_tol and s_new <= sigma_final: diff --git a/tests/unit_tests/test_criticality_search.py b/tests/unit_tests/test_criticality_search.py new file mode 100644 index 00000000000..db6ce424dcb --- /dev/null +++ b/tests/unit_tests/test_criticality_search.py @@ -0,0 +1,268 @@ +from pathlib import Path +import openmc + + +def test_keff_search_with_derivative_tallies(run_in_tmpdir): + """Test keff_search with derivative tallies enabled for fuel density perturbation.""" + # Build a simple PWR pin-cell model + fuel = openmc.Material(name='Fuel', material_id=1) + fuel.set_density('g/cm3', 10.31341) + fuel.add_element('U', 1., enrichment=1.6) + fuel.add_element('O', 2.) + + clad = openmc.Material(name='Clad', material_id=2) + clad.set_density('g/cm3', 6.55) + clad.add_element('Zr', 1.) + + coolant = openmc.Material(name='Coolant', material_id=3) + coolant.set_density('g/cm3', 0.741) + coolant.add_element('H', 2.) + coolant.add_element('O', 1.) + coolant.add_element('B', 150 * 1e-6) # 150 ppm boron + + materials = openmc.Materials([fuel, clad, coolant]) + + fuel_r = openmc.ZCylinder(r=0.39218) + clad_r = openmc.ZCylinder(r=0.45720) + min_x = openmc.XPlane(x0=-0.63, boundary_type='reflective') + max_x = openmc.XPlane(x0=+0.63, boundary_type='reflective') + min_y = openmc.YPlane(y0=-0.63, boundary_type='reflective') + max_y = openmc.YPlane(y0=+0.63, boundary_type='reflective') + + fuel_cell = openmc.Cell(fill=fuel, region=-fuel_r) + clad_cell = openmc.Cell(fill=clad, region=+fuel_r & -clad_r) + coolant_cell = openmc.Cell(fill=coolant, region=+clad_r & +min_x & -max_x & +min_y & -max_y) + + root = openmc.Universe(cells=[fuel_cell, clad_cell, coolant_cell]) + geometry = openmc.Geometry(root) + + settings = openmc.Settings() + settings.batches = 50 + settings.inactive = 5 + settings.particles = 500 + settings.run_mode = 'eigenvalue' + + bounds = [-0.63, -0.63, -10, 0.63, 0.63, 10.] + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) + settings.source = openmc.Source(space=uniform_dist) + + model = openmc.model.Model(geometry, materials, settings) + + def set_density(x): + # Modify fuel density by removing and re-adding elements + fuel_mat = model.materials[0] + fuel_mat.remove_element('U') + fuel_mat.remove_element('O') + fuel_mat.set_density('g/cm3', x) + fuel_mat.add_element('U', 1., enrichment=1.6) + fuel_mat.add_element('O', 2.) + + # Perform keff search with derivative tallies enabled + # Model class will create the required derivative tallies internally + k_tol = 5e-3 + sigma_final = 5e-3 + result = model.keff_search( + func=set_density, + x0=9.0, + x1=11.0, + target=1.17, + k_tol=k_tol, + sigma_final=sigma_final, + x_min=5.0, + x_max=12.0, + b0=settings.batches - settings.inactive, + maxiter=10, + output=False, + run_kwargs={'cwd': Path('.')}, + use_derivative_tallies=True, + deriv_variable='density', + deriv_material=1, + ) + + # Check type of result + assert isinstance(result, openmc.model.SearchResult) + + # Check that we have function evaluation history + assert len(result.parameters) >= 2 + assert len(result.means) == len(result.parameters) + assert len(result.stdevs) == len(result.parameters) + assert len(result.batches) == len(result.parameters) + + # Check that function_calls property works + assert result.function_calls == len(result.parameters) + + # Check that total_batches property works + assert result.total_batches == sum(result.batches) + assert result.total_batches > 0 + + # If converged, check tolerances (but don't fail if not converged due to limited iterations) + if result.converged: + final_keff = result.means[-1] + 1.17 # Add back target since means are (keff - target) + final_sigma = result.stdevs[-1] + assert abs(final_keff - 1.17) <= k_tol, \ + f"Final keff {final_keff:.5f} not within k_tol {k_tol}" + assert final_sigma <= sigma_final, \ + f"Final uncertainty {final_sigma:.5f} exceeds sigma_final {sigma_final}" + + +def test_keff_search_with_nuclide_density_derivatives(run_in_tmpdir): + """Test keff_search with nuclide density derivatives for boron concentration.""" + # Build a simple PWR pin-cell model + fuel = openmc.Material(name='Fuel', material_id=1) + fuel.set_density('g/cm3', 10.31341) + fuel.add_element('U', 1., enrichment=1.6) + fuel.add_element('O', 2.) + + clad = openmc.Material(name='Clad', material_id=2) + clad.set_density('g/cm3', 6.55) + clad.add_element('Zr', 1.) + + coolant = openmc.Material(name='Coolant', material_id=3) + coolant.set_density('g/cm3', 0.741) + coolant.add_element('H', 2.) + coolant.add_element('O', 1.) + coolant.add_element('B', 1000 * 1e-6) # 1000 ppm boron + + materials = openmc.Materials([fuel, clad, coolant]) + + fuel_r = openmc.ZCylinder(r=0.39218) + clad_r = openmc.ZCylinder(r=0.45720) + min_x = openmc.XPlane(x0=-0.63, boundary_type='reflective') + max_x = openmc.XPlane(x0=+0.63, boundary_type='reflective') + min_y = openmc.YPlane(y0=-0.63, boundary_type='reflective') + max_y = openmc.YPlane(y0=+0.63, boundary_type='reflective') + + fuel_cell = openmc.Cell(fill=fuel, region=-fuel_r) + clad_cell = openmc.Cell(fill=clad, region=+fuel_r & -clad_r) + coolant_cell = openmc.Cell(fill=coolant, region=+clad_r & +min_x & -max_x & +min_y & -max_y) + + root = openmc.Universe(cells=[fuel_cell, clad_cell, coolant_cell]) + geometry = openmc.Geometry(root) + + settings = openmc.Settings() + settings.batches = 50 + settings.inactive = 5 + settings.particles = 500 + settings.run_mode = 'eigenvalue' + + bounds = [-0.63, -0.63, -10, 0.63, 0.63, 10.] + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) + settings.source = openmc.Source(space=uniform_dist) + + model = openmc.model.Model(geometry, materials, settings) + + def set_boron_ppm(x): + # Modify boron concentration by removing and re-adding all elements + x = max(x, 0.0) # Ensure positive + coolant_mat = model.materials[2] + coolant_mat.remove_element('H') + coolant_mat.remove_element('O') + coolant_mat.remove_element('B') + coolant_mat.set_density('g/cm3', 0.741) + coolant_mat.add_element('H', 2.) + coolant_mat.add_element('O', 1.) + coolant_mat.add_element('B', x * 1e-6) + + # Perform keff search with nuclide density derivatives + # Model class will create the required derivative tallies internally + k_tol = 5e-3 + sigma_final = 5e-3 + result = model.keff_search( + func=set_boron_ppm, + x0=500.0, + x1=1500.0, + target=1.20, + k_tol=k_tol, + sigma_final=sigma_final, + x_min=0.1, + b0=settings.batches - settings.inactive, + maxiter=10, + output=False, + run_kwargs={'cwd': Path('.')}, + use_derivative_tallies=True, + deriv_variable='nuclide_density', + deriv_material=3, + deriv_nuclide='B10', + ) + + # Check type of result + assert isinstance(result, openmc.model.SearchResult) + + # Check that we have function evaluation history + assert len(result.parameters) >= 2 + assert len(result.means) == len(result.parameters) + assert len(result.stdevs) == len(result.parameters) + assert len(result.batches) == len(result.parameters) + + # Check that function_calls property works + assert result.function_calls == len(result.parameters) + + # Check that total_batches property works + assert result.total_batches == sum(result.batches) + assert result.total_batches > 0 + + # If converged, check tolerances (but don't fail if not converged due to limited iterations) + if result.converged: + final_keff = result.means[-1] + 1.20 # Add back target since means are (keff - target) + final_sigma = result.stdevs[-1] + assert abs(final_keff - 1.20) <= k_tol, \ + f"Final keff {final_keff:.5f} not within k_tol {k_tol}" + assert final_sigma <= sigma_final, \ + f"Final uncertainty {final_sigma:.5f} exceeds sigma_final {sigma_final}" + + + +def test_keff_search(run_in_tmpdir): + """Test the Model.keff_search method""" + + # Create model of a sphere of U235 + mat = openmc.Material() + mat.set_density('g/cm3', 18.9) + mat.add_nuclide('U235', 1.0) + sphere = openmc.Sphere(r=10.0, boundary_type='vacuum') + cell = openmc.Cell(fill=mat, region=-sphere) + geometry = openmc.Geometry([cell]) + settings = openmc.Settings(particles=1000, inactive=10, batches=30) + model = openmc.Model(geometry=geometry, settings=settings) + + # Define function to modify sphere radius + def modify_radius(radius): + sphere.r = radius + + # Perform keff search + k_tol = 4e-3 + sigma_final = 2e-3 + result = model.keff_search( + func=modify_radius, + x0=6.0, + x1=9.0, + k_tol=k_tol, + sigma_final=sigma_final, + output=True, + ) + + final_keff = result.means[-1] + 1.0 # Add back target since means are (keff - target) + final_sigma = result.stdevs[-1] + + # Check for convergence and that tolerances are met + assert result.converged, "keff_search did not converge" + assert abs(final_keff - 1.0) <= k_tol, \ + f"Final keff {final_keff:.5f} not within k_tol {k_tol}" + assert final_sigma <= sigma_final, \ + f"Final uncertainty {final_sigma:.5f} exceeds sigma_final {sigma_final}" + + # Check type of result + assert isinstance(result, openmc.model.SearchResult) + + # Check that we have function evaluation history + assert len(result.parameters) >= 2 + assert len(result.means) == len(result.parameters) + assert len(result.stdevs) == len(result.parameters) + assert len(result.batches) == len(result.parameters) + + # Check that function_calls property works + assert result.function_calls == len(result.parameters) + + # Check that total_batches property works + assert result.total_batches == sum(result.batches) + assert result.total_batches > 0 diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 3846ba4fb8d..594404da4a6 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -943,62 +943,6 @@ def test_setter_from_list(): assert isinstance(model.plots, openmc.Plots) -def test_keff_search(run_in_tmpdir): - """Test the Model.keff_search method""" - - # Create model of a sphere of U235 - mat = openmc.Material() - mat.set_density('g/cm3', 18.9) - mat.add_nuclide('U235', 1.0) - sphere = openmc.Sphere(r=10.0, boundary_type='vacuum') - cell = openmc.Cell(fill=mat, region=-sphere) - geometry = openmc.Geometry([cell]) - settings = openmc.Settings(particles=1000, inactive=10, batches=30) - model = openmc.Model(geometry=geometry, settings=settings) - - # Define function to modify sphere radius - def modify_radius(radius): - sphere.r = radius - - # Perform keff search - k_tol = 4e-3 - sigma_final = 2e-3 - result = model.keff_search( - func=modify_radius, - x0=6.0, - x1=9.0, - k_tol=k_tol, - sigma_final=sigma_final, - output=True, - ) - - final_keff = result.means[-1] + 1.0 # Add back target since means are (keff - target) - final_sigma = result.stdevs[-1] - - # Check for convergence and that tolerances are met - assert result.converged, "keff_search did not converge" - assert abs(final_keff - 1.0) <= k_tol, \ - f"Final keff {final_keff:.5f} not within k_tol {k_tol}" - assert final_sigma <= sigma_final, \ - f"Final uncertainty {final_sigma:.5f} exceeds sigma_final {sigma_final}" - - # Check type of result - assert isinstance(result, openmc.model.SearchResult) - - # Check that we have function evaluation history - assert len(result.parameters) >= 2 - assert len(result.means) == len(result.parameters) - assert len(result.stdevs) == len(result.parameters) - assert len(result.batches) == len(result.parameters) - - # Check that function_calls property works - assert result.function_calls == len(result.parameters) - - # Check that total_batches property works - assert result.total_batches == sum(result.batches) - assert result.total_batches > 0 - - def test_id_map_to_rgb(): """Test conversion of ID map to RGB image array.""" # Create a simple model