Merge branch '137-support-for-multiple-homogenization-schemes' into 'development'

Resolve "support for multiple homogenization schemes"

Closes #137

See merge request damask/DAMASK!1244

Author: MarDiehl

python/damask/_result.py

@@ -201,7 +201,7 @@ def __init__(self, fname: str | PathLike):
             self.homogenization   = f['cell_to/homogenization']['label'].astype('str')
             self._homogenizations = sorted(np.unique(self.homogenization),key=util.natural_sort)
             self.phase            = f['cell_to/phase']['label'].astype('str')
-            self._phases          = sorted(np.unique(self.phase),key=util.natural_sort)
+            self._phases          = sorted(np.unique(self.phase[self.phase != '']), key=util.natural_sort)
 
             fields = []
             for kind, labels in zip(['phase','homogenization'],[self._phases,self._homogenizations]):

src/material.f90

@@ -195,7 +195,7 @@ subroutine parse()
       material_entry_homogenization(ce) = counterHomogenization(ho)
     end do
 
-    do co = 1, size(ph_of(ma,:)>0)
+    do co = 1, count(ph_of(ma,:)>0)
 
       v  =  v_of(ma,co)
       ph = ph_of(ma,co)

src/result.f90

@@ -103,7 +103,7 @@ subroutine result_init(restart)
   if (.not. file_exists .or. .not. restart) then
     call result_createJobFile()
     call result_addAttribute('DADF5_version_major',1)
-    call result_addAttribute('DADF5_version_minor',3)
+    call result_addAttribute('DADF5_version_minor',4)
     call get_command_argument(0,commandLine)
     call result_addAttribute('creator',trim(commandLine)//' '//DAMASK_VERSION)
     call result_addAttribute('created',now())
@@ -533,10 +533,11 @@ end subroutine result_writeVectorDataset_int
 !--------------------------------------------------------------------------------------------------
 subroutine result_mapping_phase(ID,entry,label)
 
-  integer,          dimension(:,:), intent(in) :: ID                                                !< phase ID at (co,ce)
+  integer,          dimension(:,:), intent(in) :: ID                                                !< phase ID at (co,ce), 0 for padding
   integer,          dimension(:,:), intent(in) :: entry                                             !< phase entry at (co,ce)
   character(len=*), dimension(:),   intent(in) :: label                                             !< label of each phase section
 
+  character(len=len(label)), dimension(size(ID,1),size(ID,2)) :: labelPadded                        !< phase string at (co,ce), padded with '' for empty constituents
   integer(pI64), dimension(size(entry,1),size(entry,2)) :: entryGlobal
   integer(pI64), dimension(size(label),0:worldsize-1) :: entryOffset                                !< offset in entry counting per process
   integer(MPI_INTEGER_KIND), dimension(0:worldsize-1) :: writeSize                                  !< amount of data written per process
@@ -580,15 +581,18 @@ subroutine result_mapping_phase(ID,entry,label)
   entryOffset = 0_pI64
   do co = 1, size(ID,1)
     do ce = 1, size(ID,2)
-      entryOffset(ID(co,ce),worldrank) = entryOffset(ID(co,ce),worldrank) +1_pI64
+      if (ID(co,ce) /= 0) &
+        entryOffset(ID(co,ce),worldrank) = entryOffset(ID(co,ce),worldrank) +1_pI64
     end do
   end do
   call MPI_Allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INTEGER8,MPI_SUM,MPI_COMM_WORLD,err_MPI)! get offset at each process
   call parallelization_chkerr(err_MPI)
   entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2)
+  entryGlobal = -1_pI64
   do co = 1, size(ID,1)
     do ce = 1, size(ID,2)
-      entryGlobal(co,ce) = int(entry(co,ce),pI64) -1_pI64 + entryOffset(ID(co,ce),worldrank)
+      if (ID(co,ce) /= 0) &
+        entryGlobal(co,ce) = int(entry(co,ce),pI64) -1_pI64 + entryOffset(ID(co,ce),worldrank)
     end do
   end do
 #endif
@@ -601,7 +605,7 @@ subroutine result_mapping_phase(ID,entry,label)
 ! compound type: label(ID) + entry
   call H5Tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
   call HDF5_chkerr(hdferr)
-  call H5Tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
+  call H5Tset_size_f(dt_id, int(len(label),SIZE_T), hdferr)
   call HDF5_chkerr(hdferr)
   call H5Tget_size_f(dt_id, type_size_str, hdferr)
   call HDF5_chkerr(hdferr)
@@ -652,10 +656,17 @@ subroutine result_mapping_phase(ID,entry,label)
   call H5Dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr)
   call HDF5_chkerr(hdferr)
 
-  call H5Dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
+  labelPadded = ''
+  do co = 1, size(ID,1)
+    do ce = 1, size(ID,2)
+      if (ID(co,ce) > 0) &
+        labelPadded(co,ce) = label(ID(co,ce))
+    end do
+  end do
+  call H5Dwrite_f(dset_id, label_id, labelPadded, &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   call HDF5_chkerr(hdferr)
-  call H5Dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
+  call H5Dwrite_f(dset_id, entry_id, entryGlobal, &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   call HDF5_chkerr(hdferr)
 
@@ -750,7 +761,7 @@ subroutine result_mapping_homogenization(ID,entry,label)
 ! compound type: label(ID) + entry
   call H5Tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
   call HDF5_chkerr(hdferr)
-  call H5Tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
+  call H5Tset_size_f(dt_id, int(len(label),SIZE_T), hdferr)
   call HDF5_chkerr(hdferr)
   call H5Tget_size_f(dt_id, type_size_str, hdferr)
   call HDF5_chkerr(hdferr)
@@ -801,10 +812,10 @@ subroutine result_mapping_homogenization(ID,entry,label)
   call H5Dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr)
   call HDF5_chkerr(hdferr)
 
-  call H5Dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
+  call H5Dwrite_f(dset_id, label_id, label(ID), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   call HDF5_chkerr(hdferr)
-  call H5Dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
+  call H5Dwrite_f(dset_id, entry_id, entryGlobal, &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   call HDF5_chkerr(hdferr)
 

tests/integration/test_homogenization.py

@@ -1,6 +1,7 @@
 # SPDX-License-Identifier: AGPL-3.0-or-later
 import copy
 
+import h5py
 import pytest
 import numpy as np
 
@@ -60,6 +61,42 @@ def test_homogenization_equivalent(res_path,tmp_path,copy_files,phase,homogeniza
            0.99 < results[phase]['P']/P_11 < 1.01
 
 
+def test_multiple_homogenization_schemes_and_constituent_numbers(res_path,tmp_path,copy_files,np_rng):
+    copy_files(res_path,tmp_path)
+
+    load = 'compressionX'
+    material = 'material'
+    job = f'combined_{load}'
+
+    schemes = ['RGC','none','Taylor2','RGC']
+    damask.GeomGrid(np.arange(len(schemes)).reshape(len(schemes),1,1),
+                    np.array([len(schemes),1,1])*1e-4).save(tmp_path/'combined.vti')
+
+    material_config = damask.ConfigMaterial.load(tmp_path/f'{material}.yaml')
+    N_constituents = [material_config['homogenization'][ho]['N_constituents'] for ho in schemes]
+    assert len(set(N_constituents)) > 1 # make sure we also test different numbers of constituents
+    O = damask.Rotation.from_random((len(schemes),max(N_constituents)),rng_seed=np_rng)
+    material_config['material'] = [
+        dict(homogenization=ho,
+             constituents=[
+                 dict(phase='Isotropic' if ho=='none' or (ho=='Taylor2' and co==0) else 'Phenopowerlaw',
+                      v=1/N,
+                      O=O[ma,co].as_quaternion().tolist()
+                     ) for co in range(N)]
+            ) for ma,(ho,N) in enumerate(zip(schemes,N_constituents))]
+    material_config.save(tmp_path/f'{material}.yaml')
+
+    damask.util.run(f'damask_grid -l {load}.yaml -g combined.vti -m {material}.yaml -j {job}',
+                    wd=tmp_path)
+    result = damask.Result(tmp_path/f'{job}.hdf5')
+    assert result.phases == ['Isotropic', 'Phenopowerlaw']
+    with h5py.File(tmp_path/f'{job}.hdf5','r') as f:
+        label = f['cell_to/phase']['label'][()].astype(str)
+        entry = f['cell_to/phase']['entry'][()]
+
+    assert np.array_equal(np.sum(label != '',axis=1),N_constituents)
+    assert np.all(entry[label == ''] == -1)
+
 def test_homogenization_many(tmp_path,res_path,copy_files,np_rng):
     grid = damask.GeomGrid(np.zeros([2,1,1],int),np.array([2,1,1])*1e-4)
     grid.save(tmp_path/'small')