Adding material decay gamma dose

openmc/material.py

@@ -369,6 +369,110 @@ def get_decay_photon_energy(
 
         return combined
 
+
+    def get_decay_photon_dose(
+        self,
+        n_samples: int = 10000,
+        distance: float | None = None
+    ) -> float:
+        """Return decay photon dose rate from unstable nuclides.
+
+        .. versionadded:: 0.15.3
+
+        Parameters
+        ----------
+        n_samples : int, optional
+            Number of energy samples to use for Monte Carlo integration.
+            For Discrete distributions, exact calculation is used instead of
+            sampling. Defaults to 10000.
+        distance : float, optional
+            Distance from the source in cm at which to calculate dose rate.
+            If not specified, the dose is calculated at the surface of an
+            equivalent sphere based on the material's volume.
+
+        Returns
+        -------
+        float
+            Dose rate in pSv/s at the specified distance or at the surface
+            of an equivalent sphere
+
+        """
+        cv.check_type('n_samples', n_samples, int)
+        cv.check_greater_than('n_samples', n_samples, 0)
+        if distance is not None:
+            cv.check_type('distance', distance, Real)
+            cv.check_greater_than('distance', distance, 0.0)
+
+        # Get the decay photon energy distribution
+        dist = self.get_decay_photon_energy(units='Bq')
+
+        # If no photon sources, return zero dose
+        if dist is None:
+            return 0.0
+
+        # Get dose coefficients for isotropic geometry
+        energy_bins, dose_coeffs = openmc.data.dose_coefficients(
+            particle="photon",
+            geometry="AP"
+        )
+
+        # For pure Discrete distributions, compute mean dose coefficient exactly
+        # without sampling to improve performance
+        if isinstance(dist, Discrete):
+            # Check if maximum energy exceeds the dose coefficient range
+            max_energy = np.max(dist.x)
+            if max_energy > energy_bins[-1]:
+                raise ValueError(
+                    f"Maximum photon energy {max_energy:.2e} eV exceeds "
+                    f"dose coefficient upper limit of {energy_bins[-1]:.2e} eV"
+                )
+
+            # Interpolate dose coefficients at discrete energy points
+            interpolated_coeffs = np.interp(dist.x, energy_bins, dose_coeffs)
+
+            # Calculate weighted mean dose coefficient
+            mean_dose_coeff = np.sum(interpolated_coeffs * dist.p) / dist.integral()
+
+        else:
+            # For other distributions (Mixture, etc.), use Monte Carlo sampling
+            # Sample energies from the distribution with fixed seed for reproducibility
+            energies = dist.sample(n_samples=n_samples, seed=1)
+
+            # Check if maximum sampled energy exceeds the dose coefficient range
+            max_energy = np.max(energies)
+            if max_energy > energy_bins[-1]:
+                raise ValueError(
+                    f"Maximum sampled photon energy {max_energy:.2e} eV exceeds "
+                    f"dose coefficient upper limit of {energy_bins[-1]:.2e} eV"
+                )
+
+            # Interpolate dose coefficients at sampled energies
+            # dose_coeffs are in pSv·cm²
+            interpolated_coeffs = np.interp(energies, energy_bins, dose_coeffs)
+
+            # Calculate mean dose coefficient
+            mean_dose_coeff = np.mean(interpolated_coeffs)
+
+        # Total activity in Bq (photons/s)
+        total_activity = dist.integral()
+
+        # If distance not specified, use surface of equivalent sphere
+        if distance is None:
+            if self.volume is None:
+                raise ValueError(
+                    "Either distance must be specified or material volume "
+                    "must be set to calculate dose rate."
+                )
+            # Calculate radius of equivalent sphere: V = (4/3)πr³
+            # r = (3V/4π)^(1/3)
+            distance = (3.0 * self.volume / (4.0 * np.pi)) ** (1.0 / 3.0)
+        # Dose rate at specified distance in pSv/s
+        # dose_rate = (coeff * activity) / (4π·r²)
+        dose_rate_pSv_per_s = mean_dose_coeff * total_activity / (4 * np.pi * distance**2)
+
+        return dose_rate_pSv_per_s
+
+
     @classmethod
     def from_hdf5(cls, group: h5py.Group) -> Material:
         """Create material from HDF5 group

tests/unit_tests/test_material.py

@@ -767,3 +767,46 @@ def test_mean_free_path():
     mat2.add_nuclide('Pb208', 1.0)
     mat2.set_density('g/cm3', 11.34)
     assert mat2.mean_free_path(energy=14e6) == pytest.approx(5.65, abs=1e-2)
+
+
+def test_get_decay_photon_dose():
+    """Test decay photon dose rate calculation"""
+    # Set chain file for testing
+    openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml'
+
+    # Test with stable nuclides this should return zero
+    mat_stable = openmc.Material()
+    mat_stable.add_nuclide('Fe56', 1.0)
+    mat_stable.set_density('g/cm3', 7.87)
+    mat_stable.volume = 1.0
+    assert mat_stable.get_decay_photon_dose() == 0.0
+
+    # Test with unstable nuclide this should return positive dose
+    mat_unstable = openmc.Material()
+    mat_unstable.add_nuclide('I135', 1.0)
+    mat_unstable.set_density('g/cm3', 1.0)
+    mat_unstable.volume = 100.0
+    dose = mat_unstable.get_decay_photon_dose(n_samples=10)
+    assert dose > 0.0
+    assert isinstance(dose, float)
+
+    # Test distance dependence (1/r² law)
+    dose_10cm = mat_unstable.get_decay_photon_dose(distance=10.0)
+    dose_20cm = mat_unstable.get_decay_photon_dose(distance=20.0)
+    ratio = dose_10cm / dose_20cm
+    assert 3.5 < ratio < 4.5  # Should be ~4 due to 1/r²
+
+    expected_radius = (3.0 * 100.0 / (4.0 * np.pi)) ** (1.0 / 3.0)
+    dose_default = mat_unstable.get_decay_photon_dose()
+    dose_at_radius = mat_unstable.get_decay_photon_dose(distance=expected_radius)
+    assert abs(dose_default - dose_at_radius) / dose_default < 1e-10
+
+    # Test invalid inputs
+    with pytest.raises(ValueError):
+        mat_unstable.get_decay_photon_dose(n_samples=-100)
+    with pytest.raises(ValueError):
+        mat_unstable.get_decay_photon_dose(n_samples=0)
+    with pytest.raises(ValueError):
+        mat_unstable.get_decay_photon_dose(distance=-10.0)
+    with pytest.raises(ValueError):
+        mat_unstable.get_decay_photon_dose(distance=0.0)