Support setting density in g/cm3 through C API

docs/source/capi/index.rst

@@ -333,12 +333,14 @@ Functions
    :return: Return status (negative if an error occurred)
    :rtype: int
 
-.. c:function:: int openmc_material_set_density(int32_t index, double density)
+.. c:function:: int openmc_material_set_density(int32_t index, double density, const char* units)
 
    Set the density of a material.
 
    :param int32_t index: Index in the materials array
-   :param double density: Density of the material in atom/b-cm
+   :param double density: Density of the material
+   :param units: Units for density
+   :type units: const char*
    :return: Return status (negative if an error occurs)
    :rtype: int
 

include/openmc/capi.h

@@ -60,7 +60,7 @@ extern "C" {
   int openmc_material_get_id(int32_t index, int32_t* id);
   int openmc_material_get_fissionable(int32_t index, bool* fissionable);
   int openmc_material_get_volume(int32_t index, double* volume);
-  int openmc_material_set_density(int32_t index, double density);
+  int openmc_material_set_density(int32_t index, double density, const char* units);
   int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density);
   int openmc_material_set_id(int32_t index, int32_t id);
   int openmc_material_set_volume(int32_t index, double volume);

openmc/capi/material.py

@@ -35,7 +35,7 @@ _dll.openmc_material_get_densities.errcheck = _error_handler
 _dll.openmc_material_get_volume.argtypes = [c_int32, POINTER(c_double)]
 _dll.openmc_material_get_volume.restype = c_int
 _dll.openmc_material_get_volume.errcheck = _error_handler
-_dll.openmc_material_set_density.argtypes = [c_int32, c_double]
+_dll.openmc_material_set_density.argtypes = [c_int32, c_double, c_char_p]
 _dll.openmc_material_set_density.restype = c_int
 _dll.openmc_material_set_density.errcheck = _error_handler
 _dll.openmc_material_set_densities.argtypes = [
@@ -178,16 +178,18 @@ class Material(_FortranObjectWithID):
         """
         _dll.openmc_material_add_nuclide(self._index, name.encode(), density)
 
-    def set_density(self, density):
+    def set_density(self, density, units='atom/b-cm'):
         """Set density of a material.
 
         Parameters
         ----------
         density : float
-            Density in atom/b-cm
+            Density
+        units : {'atom/b-cm', 'g/cm3'}
+            Units for density
 
         """
-        _dll.openmc_material_set_density(self._index, density)
+        _dll.openmc_material_set_density(self._index, density, units.encode())
 
     def set_densities(self, nuclides, densities):
         """Set the densities of a list of nuclides in a material

src/material_header.F90

@@ -12,7 +12,7 @@ module material_header
   use sab_header
   use simulation_header, only: log_spacing
   use stl_vector, only: VectorReal, VectorInt
-  use string, only: to_str
+  use string, only: to_str, to_f_string
 
   implicit none
 
@@ -166,34 +166,52 @@ contains
     call material_set_fissionable_c(this % ptr, logical(fissionable, C_BOOL))
   end subroutine material_set_fissionable
 
-  function material_set_density(this, density) result(err)
+  function material_set_density(this, density, units) result(err)
     class(Material), intent(inout) :: this
     real(8), intent(in) :: density
+    character(*), intent(in) :: units
     integer :: err
 
     integer :: i
     real(8) :: sum_percent
     real(8) :: awr
+    real(8) :: previous_density_gpcc
+    real(8) :: f
 
     if (allocated(this % atom_density)) then
-      ! Set total density based on value provided
-      this % density = density
-
-      ! Determine normalized atom percents
-      sum_percent = sum(this % atom_density)
-      this % atom_density(:) = this % atom_density / sum_percent
-
-      ! Recalculate nuclide atom densities based on given density
-      this % atom_density(:) = density * this % atom_density
-
-      ! Calculate density in g/cm^3.
-      this % density_gpcc = ZERO
-      do i = 1, this % n_nuclides
-        awr = nuclides(this % nuclide(i)) % awr
-        this % density_gpcc = this % density_gpcc &
-             + this % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO
-      end do
       err = 0
+      select case (units)
+      case ('atom/b-cm')
+        ! Set total density based on value provided
+        this % density = density
+
+        ! Determine normalized atom percents
+        sum_percent = sum(this % atom_density)
+        this % atom_density(:) = this % atom_density / sum_percent
+
+        ! Recalculate nuclide atom densities based on given density
+        this % atom_density(:) = density * this % atom_density
+
+        ! Calculate density in g/cm^3.
+        this % density_gpcc = ZERO
+        do i = 1, this % n_nuclides
+          awr = nuclides(this % nuclide(i)) % awr
+          this % density_gpcc = this % density_gpcc &
+              + this % atom_density(i) * awr * MASS_NEUTRON / N_AVOGADRO
+        end do
+      case ('g/cm3', 'g/cc')
+        ! Determine factor by which to change densities
+        previous_density_gpcc = this % density_gpcc
+        f = density / previous_density_gpcc
+
+        ! Update densities
+        this % density_gpcc = density
+        this % density = f * this % density
+        this % atom_density(:) = f * this % atom_density(:)
+      case default
+        err = E_INVALID_ARGUMENT
+        call set_errmsg("Invalid units '" // trim(units) // "' specified.")
+      end select
     else
       err = E_ALLOCATE
       call set_errmsg("Material atom density array hasn't been allocated.")
@@ -738,16 +756,22 @@ contains
   end function openmc_material_set_id
 
 
-  function openmc_material_set_density(index, density) result(err) bind(C)
-    ! Set the total density of a material in atom/b-cm
+  function openmc_material_set_density(index, density, units) result(err) bind(C)
+    ! Set the total density of a material
     integer(C_INT32_T), value, intent(in) :: index
     real(C_DOUBLE), value, intent(in) :: density
+    character(kind=C_CHAR), intent(in) :: units(*)
     integer(C_INT) :: err
 
+    character(:), allocatable :: units_
+
+    ! Convert C string to Fortran string
+    units_ = to_f_string(units)
+
     err = E_UNASSIGNED
     if (index >= 1 .and. index <= size(materials)) then
       associate (m => materials(index))
-        err = m % set_density(density)
+        err = m % set_density(density, units_)
       end associate
     else
       err = E_OUT_OF_BOUNDS
@@ -794,7 +818,7 @@ contains
         m % n_nuclides = n
 
         ! Set total density to the sum of the vector
-        err = m % set_density(sum(density))
+        err = m % set_density(sum(density), 'atom/b-cm')
 
         ! Assign S(a,b) tables
         call m % assign_sab_tables()

tests/unit_tests/test_capi.py

@@ -105,6 +105,9 @@ def test_material(capi_init):
     m.volume = 10.0
     assert m.volume == 10.0
 
+    with pytest.raises(exc.InvalidArgumentError):
+        m.set_density(1.0, 'goblins')
+
     rho = 2.25e-2
     m.set_density(rho)
     assert sum(m.densities) == pytest.approx(rho)