Merge pull request #2222 from highness-eu/mixed_ncrystal

Adding support for NCrystal materials in OpenMC

openmc/lib/__init__.py

@@ -42,6 +42,9 @@ else:
 def _dagmc_enabled():
     return c_bool.in_dll(_dll, "DAGMC_ENABLED").value
 
+def _ncrystal_enabled():
+    return c_bool.in_dll(_dll, "NCRYSTAL_ENABLED").value
+
 def _coord_levels():
     return c_int.in_dll(_dll, "n_coord_levels").value
 

openmc/material.py

@@ -99,6 +99,10 @@ class Material(IDManagerMixin):
         [decay/sec].
 
         .. versionadded:: 0.13.2
+    ncrystal_cfg : str
+        NCrystal configuration string
+
+        .. versionadded:: 0.13.3
 
     """
 
@@ -118,6 +122,7 @@ class Material(IDManagerMixin):
         self._volume = None
         self._atoms = {}
         self._isotropic = []
+        self._ncrystal_cfg = None
 
         # A list of tuples (nuclide, percent, percent type)
         self._nuclides = []
@@ -140,6 +145,9 @@ class Material(IDManagerMixin):
 
         string += '{: <16}\n'.format('\tS(a,b) Tables')
 
+        if self._ncrystal_cfg:
+            string += '{: <16}=\t{}\n'.format('\tNCrystal conf', self._ncrystal_cfg)
+
         for sab in self._sab:
             string += '{: <16}=\t{}\n'.format('\tS(a,b)', sab)
 
@@ -219,6 +227,10 @@ class Material(IDManagerMixin):
     def volume(self):
         return self._volume
 
+    @property
+    def ncrystal_cfg(self):
+        return self._ncrystal_cfg
+
     @name.setter
     def name(self, name: Optional[str]):
         if name is not None:
@@ -331,6 +343,64 @@ class Material(IDManagerMixin):
 
         return material
 
+    @classmethod
+    def from_ncrystal(cls, cfg, **kwargs):
+        """Create material from NCrystal configuration string.
+
+        Density, temperature, and material composition, and (ultimately) thermal
+        neutron scattering will be automatically be provided by NCrystal based
+        on this string. The name and material_id parameters are simply passed on
+        to the Material constructor.
+
+        Parameters
+        ----------
+        cfg : str
+            NCrystal configuration string
+        **kwargs
+            Keyword arguments passed to :class:`openmc.Material`
+
+        Returns
+        -------
+        openmc.Material
+            Material instance
+
+        """
+
+        import NCrystal
+        nc_mat = NCrystal.createInfo(cfg)
+
+        def openmc_natabund(Z):
+            #nc_mat.getFlattenedComposition might need natural abundancies.
+            #This call-back function is used so NCrystal can flatten composition
+            #using OpenMC's natural abundancies. In practice this function will
+            #only get invoked in the unlikely case where a material is specified
+            #by referring both to natural elements and specific isotopes of the
+            #same element.
+            elem_name = openmc.data.ATOMIC_SYMBOL[Z]
+            return [
+                (int(iso_name[len(elem_name):]), abund)
+                for iso_name, abund in openmc.data.isotopes(elem_name)
+            ]
+
+        flat_compos = nc_mat.getFlattenedComposition(
+            preferNaturalElements=True, naturalAbundProvider=openmc_natabund)
+
+        # Create the Material
+        material = cls(temperature=nc_mat.getTemperature(), **kwargs)
+
+        for Z, A_vals in flat_compos:
+            elemname = openmc.data.ATOMIC_SYMBOL[Z]
+            for A, frac in A_vals:
+                if A:
+                    material.add_nuclide(f'{elemname}{A}', frac)
+                else:
+                    material.add_element(elemname, frac)
+
+        material.set_density('g/cm3', nc_mat.getDensity())
+        material._ncrystal_cfg = NCrystal.normaliseCfg(cfg)
+
+        return material
+
     def add_volume_information(self, volume_calc):
         """Add volume information to a material.
 
@@ -405,6 +475,9 @@ class Material(IDManagerMixin):
                   'macroscopic data-set has already been added'.format(self._id)
             raise ValueError(msg)
 
+        if self._ncrystal_cfg is not None:
+            raise ValueError("Cannot add nuclides to NCrystal material")
+
         # If nuclide name doesn't look valid, give a warning
         try:
             Z, _, _ = openmc.data.zam(nuclide)
@@ -609,6 +682,9 @@ class Material(IDManagerMixin):
             raise ValueError("Element name should be given by the "
                              "element's symbol or name, e.g., 'Zr', 'zirconium'")
 
+        if self._ncrystal_cfg is not None:
+            raise ValueError("Cannot add elements to NCrystal material")
+
         # Allow for element identifier to be given as a symbol or name
         if len(element) > 2:
             el = element.lower()
@@ -990,7 +1066,7 @@ class Material(IDManagerMixin):
             activity[nuclide] = inv_seconds * 1e24 * atoms_per_bcm * multiplier
 
         return activity if by_nuclide else sum(activity.values())
-    
+
     def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False):
         """Returns the decay heat of the material or for each nuclide in the
         material in units of [W], [W/g] or [W/cm3].
@@ -1024,16 +1100,16 @@ class Material(IDManagerMixin):
             multiplier = 1
         elif units == 'W/g':
             multiplier = 1.0 / self.get_mass_density()
-        
-        decayheat = {} 
+
+        decayheat = {}
         for nuclide, atoms_per_bcm in self.get_nuclide_atom_densities().items():
             decay_erg = openmc.data.decay_energy(nuclide)
             inv_seconds = openmc.data.decay_constant(nuclide)
             decay_erg *= openmc.data.JOULE_PER_EV
             decayheat[nuclide] = inv_seconds * decay_erg * 1e24 * atoms_per_bcm * multiplier
 
-        return decayheat if by_nuclide else sum(decayheat.values()) 
-    
+        return decayheat if by_nuclide else sum(decayheat.values())
+
     def get_nuclide_atoms(self):
         """Return number of atoms of each nuclide in the material
 
@@ -1182,6 +1258,14 @@ class Material(IDManagerMixin):
         if self._volume:
             element.set("volume", str(self._volume))
 
+        if self._ncrystal_cfg:
+            if self._sab:
+                raise ValueError("NCrystal materials are not compatible with S(a,b).")
+            if self._macroscopic is not None:
+                raise ValueError("NCrystal materials are not compatible with macroscopic cross sections.")
+
+            element.set("cfg", str(self._ncrystal_cfg))
+
         # Create temperature XML subelement
         if self.temperature is not None:
             element.set("temperature", str(self.temperature))
@@ -1566,4 +1650,4 @@ class Materials(cv.CheckedList):
         tree = ET.parse(path)
         root = tree.getroot()
 
-        return cls.from_xml_element(root)
+        return cls.from_xml_element(root)