Merge remote-tracking branch 'origin/development' into thermal-partioning

PRIVATE

@@ -1 +1 @@
-Subproject commit 857b994fb7222ab15a2b8c4ded2bba8787d7feb6
+Subproject commit 4c8116ba3b9e9fbb325a580705028e8310139117

examples/config/homogenization/8grains.yaml

@@ -0,0 +1 @@
+N_constituents: 8

examples/config/homogenization/Homogenization_Isostrain_Parallel3.config

@@ -1,4 +0,0 @@
-[Parallel3]
-mech	        isostrain
-nconstituents	3
-mapping	        sum		# or 'parallel'

examples/config/homogenization/bicrystal.yaml

@@ -0,0 +1 @@
+N_constituents: 2

examples/config/homogenization/damage/pass_direct.yaml

@@ -0,0 +1,3 @@
+# For single point calculations, requires N_constituents = 1
+type: pass
+output: ['T']

examples/config/homogenization/direct.yaml

@@ -0,0 +1 @@
+N_constituents: 1

examples/config/homogenization/mechanical/RGC_8grains.yaml

@@ -1,10 +1,8 @@
-8Grains:
+# For Relaxed Grain Cluster homogenization, requires N_constituents = 8
-  N_constituents: 8
+type: RGC
-  mechanical:
+D_alpha: [4.0e-06, 4.0e-06, 2.0e-06]
-    type: RGC
+a_g: [0.0, 0.0, 0.0]
-    D_alpha: [4.0e-06, 4.0e-06, 2.0e-06]
+c_alpha: 2.0
-    a_g: [0.0, 0.0, 0.0]
+cluster_size: [2, 2, 2]
-    c_alpha: 2.0
+output: [M, Delta_V, avg_dot_a, max_dot_a]
-    cluster_size: [2, 2, 2]
+xi_alpha: 10.0
-    output: [M, Delta_V, avg_dot_a, max_dot_a]
-    xi_alpha: 10.0

examples/config/homogenization/mechanical/isostrain_Taylor2.yaml

@@ -1,3 +0,0 @@
-Taylor2:
-  N_constituents: 2
-  mechanical: {type: isostrain}

examples/config/homogenization/mechanical/isostrain_polycrystal.yaml

@@ -0,0 +1,3 @@
+# For Taylor homogenization with N_constituents > 1
+type: isostrain
+output: ['F', 'P']

examples/config/homogenization/mechanical/pass_direct.yaml

@@ -0,0 +1,3 @@
+# For single point calculations, requires N_constituents = 1
+type: pass
+output: ['F', 'P']

examples/config/homogenization/thermal/isotemperature_polycrystal.yaml

@@ -0,0 +1,3 @@
+# For homogenization with N_constituents > 1
+type: isotemperature
+output: ['T']

examples/config/homogenization/thermal/pass_direct.yaml

@@ -0,0 +1,3 @@
+# For single point calculations, requires N_constituents = 1
+type: pass
+output: ['T']

examples/config/phase/thermal/Al.yaml

@@ -5,6 +5,8 @@ references:
     fit to Tab. 3.4.1 (RRR=1000, T_min=200K, T_max=900K)
   - https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
 
+output: [T]
+
 K_11: 2.380e+2
 K_11,T: 2.068e-3
 K_11,T^2: -7.765e-5

examples/config/phase/thermal/Cu.yaml

@@ -5,6 +5,8 @@ references:
     fit to Tab. 2.4.1 (RRR=1000, T_min=200K, T_max=1000K)
   - https://www.mit.edu/~6.777/matprops/copper.htm
 
+output: [T]
+
 K_11: 4.039e+2
 K_11,T: -8.119e-2
 K_11,T^2: 1.454e-5

examples/config/phase/thermal/Fe.yaml

@@ -5,6 +5,8 @@ references:
     fit to Tab. 4.4.1 (RRR=300, T_min=200K, T_max=1000K)
   - https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
 
+output: [T]
+
 K_11: 8.055e+1
 K_11,T: -1.051e-1
 K_11,T^2: 5.464e-5

examples/config/phase/thermal/Ni.yaml

@@ -5,6 +5,8 @@ references:
     fit to Tab. 35R (T_min=150K, T_max=500K)
   - https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
 
+output: [T]
+
 K_11: 9.132e+1
 K_11,T: -1.525e-1
 K_11,T^2: 3.053e-4

examples/config/phase/thermal/Sn-beta.yaml

@@ -5,6 +5,8 @@ references:
     fit to Tab. 61R (T_min=100K, T_max=400K)
   - https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
 
+output: [T]
+
 K_11: 7.414e+1
 K_11,T: -6.465e-2
 K_11,T^2: 2.066e-4

examples/config/phase/thermal/W.yaml

@@ -5,6 +5,8 @@ references:
     fit to Tab. 5.4.1 (RRR=300, T_min=200K, T_max=1000K)
   - https://www.mit.edu/~6.777/matprops/tungsten.htm
 
+output: [T]
+
 K_11: 1.758e+2
 K_11,T: -1.605e-1
 K_11,T^2: 1.160e-4

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha6-4-gca6a3e786
+v3.0.0-alpha6-14-g3657b2316

python/damask/_crystal.py

@@ -29,7 +29,29 @@ lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
 
 
 class Crystal():
-    """Crystal lattice."""
+    """
+    Representation of a crystal as (general) crystal family or (more specific) as a scaled Bravais lattice.
+
+    Examples
+    --------
+    Cubic crystal family:
+
+    >>> import damask
+    >>> cubic = damask.Crystal(family='cubic')
+    >>> cubic
+    Crystal family: cubic
+
+    Body-centered cubic Bravais lattice with parameters of iron:
+
+    >>> import damask
+    >>> Fe = damask.Crystal(lattice='cI', a=287e-12)
+    >>> Fe
+    Crystal family: cubic
+    Bravais lattice: cI
+    a=2.87e-10 m, b=2.87e-10 m, c=2.87e-10 m
+    α=90°, β=90°, γ=90°
+
+    """
 
     def __init__(self, *,
                  family: CrystalFamily = None,
@@ -38,7 +60,7 @@ class Crystal():
                  alpha: float = None, beta: float = None, gamma: float = None,
                  degrees: bool = False):
         """
-        Representation of crystal in terms of crystal family or Bravais lattice.
+        New representation of a crystal.
 
         Parameters
         ----------
@@ -114,10 +136,10 @@ class Crystal():
         """Represent."""
         family = f'Crystal family: {self.family}'
         return family if self.lattice is None else \
-               '\n'.join([family,
+               util.srepr([family,
-                          f'Bravais lattice: {self.lattice}',
+                           f'Bravais lattice: {self.lattice}',
-                          'a={:.5g}m, b={:.5g}m, c={:.5g}m'.format(*self.parameters[:3]),
+                           'a={:.5g} m, b={:.5g} m, c={:.5g} m'.format(*self.parameters[:3]),
-                          'α={:.5g}°, β={:.5g}°, γ={:.5g}°'.format(*np.degrees(self.parameters[3:]))])
+                           'α={:.5g}°, β={:.5g}°, γ={:.5g}°'.format(*np.degrees(self.parameters[3:]))])
 
 
     def __eq__(self,
@@ -323,7 +345,8 @@ class Crystal():
         Parameters
         ----------
         direction|plane : numpy.ndarray, shape (...,3)
-            Vector along direction or plane normal.
+            Real space vector along direction or
+            reciprocal space vector along plane normal.
 
         Returns
         -------
@@ -344,7 +367,7 @@ class Crystal():
                  uvw: FloatSequence = None,
                  hkl: FloatSequence = None) -> np.ndarray:
         """
-        Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl).
+        Calculate crystal frame vector corresponding to lattice direction [uvw] or plane normal (hkl).
 
         Parameters
         ----------
@@ -354,7 +377,24 @@ class Crystal():
         Returns
         -------
         vector : numpy.ndarray, shape (...,3)
-            Crystal frame vector along [uvw] direction or (hkl) plane normal.
+            Crystal frame vector in real space along [uvw] direction or
+            in reciprocal space along (hkl) plane normal.
+
+        Examples
+        --------
+        Crystal frame vector (real space) of Magnesium corresponding to [1,1,0] direction:
+
+        >>> import damask
+        >>> Mg = damask.Crystal(lattice='hP', a=321e-12, c=521e-12)
+        >>> Mg.to_frame(uvw=[1, 1, 0])
+        array([1.60500000e-10, 2.77994155e-10, 0.00000000e+00])
+
+        Crystal frame vector (reciprocal space) of Titanium along (1,0,0) plane normal:
+
+        >>> import damask
+        >>> Ti = damask.Crystal(lattice='hP', a=0.295e-9, c=0.469e-9)
+        >>> Ti.to_frame(hkl=(1, 0, 0))
+        array([ 3.38983051e+09,  1.95711956e+09, -4.15134508e-07])
 
         """
         if (uvw is not None) ^ (hkl is None):

src/homogenization.f90

@@ -6,9 +6,10 @@
 !--------------------------------------------------------------------------------------------------
 module homogenization
   use prec
+  use math
+  use constants
   use IO
   use config
-  use math
   use material
   use phase
   use discretization
@@ -32,7 +33,7 @@ module homogenization
     HOMOGENIZATION_RGC_ID
   end enum
 
-    type(tState),        allocatable, dimension(:), public :: &
+  type(tState),        allocatable, dimension(:), public :: &
     homogState, &
     damageState_h
 
@@ -66,15 +67,13 @@ module homogenization
 !--------------------------------------------------------------------------------------------------
   interface
 
-    module subroutine mechanical_init(num_homog)
+    module subroutine mechanical_init()
-      class(tNode), pointer, intent(in) :: &
-        num_homog                                                                                   !< pointer to mechanical homogenization numerics data
     end subroutine mechanical_init
 
-    module subroutine thermal_init
+    module subroutine thermal_init()
     end subroutine thermal_init
 
-    module subroutine damage_init
+    module subroutine damage_init()
     end subroutine damage_init
 
     module subroutine mechanical_partition(subF,ce)
@@ -204,15 +203,15 @@ subroutine homogenization_init()
 
   allocate(homogState      (size(material_name_homogenization)))
   allocate(damageState_h   (size(material_name_homogenization)))
-  call material_parseHomogenization()
+  call parseHomogenization()
 
   num_homog        => config_numerics%get('homogenization',defaultVal=emptyDict)
   num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
 
-  num%nMPstate  = num_homogGeneric%get_asInt('nMPstate',defaultVal=10)
+  num%nMPstate = num_homogGeneric%get_asInt('nMPstate',defaultVal=10)
   if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate')
 
-  call mechanical_init(num_homog)
+  call mechanical_init()
   call thermal_init()
   call damage_init()
 
@@ -322,13 +321,13 @@ subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolvin
   elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
     IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
       ce = (el-1)*discretization_nIPs + ip
-    ho = material_homogenizationID(ce)
+      ho = material_homogenizationID(ce)
-    do co = 1, homogenization_Nconstituents(ho)
+      do co = 1, homogenization_Nconstituents(ho)
-      call crystallite_orientations(co,ip,el)
+        call crystallite_orientations(co,ip,el)
-      enddo
+      end do
-    call mechanical_homogenize(Delta_t,ce)
+      call mechanical_homogenize(Delta_t,ce)
-    enddo IpLooping3
+    end do IpLooping3
-  enddo elementLooping3
+  end do elementLooping3
   !$OMP END PARALLEL DO
 
 
@@ -446,7 +445,7 @@ end subroutine homogenization_restartRead
 !--------------------------------------------------------------------------------------------------
 !> @brief parses the homogenization part from the material configuration
 !--------------------------------------------------------------------------------------------------
-subroutine material_parseHomogenization
+subroutine parseHomogenization
 
   class(tNode), pointer :: &
     material_homogenization, &
@@ -458,8 +457,8 @@ subroutine material_parseHomogenization
 
   material_homogenization => config_material%get('homogenization')
 
-  allocate(thermal_type(size(material_name_homogenization)),        source=THERMAL_isothermal_ID)
+  allocate(thermal_type(size(material_name_homogenization)),source=THERMAL_isothermal_ID)
-  allocate(damage_type (size(material_name_homogenization)),        source=DAMAGE_none_ID)
+  allocate(damage_type (size(material_name_homogenization)),source=DAMAGE_none_ID)
 
   do h=1, size(material_name_homogenization)
     homog => material_homogenization%get(h)
@@ -467,7 +466,7 @@ subroutine material_parseHomogenization
     if (homog%contains('thermal')) then
       homogThermal => homog%get('thermal')
         select case (homogThermal%get_asString('type'))
-          case('pass')
+          case('pass','isotemperature')
             thermal_type(h) = THERMAL_conduction_ID
           case default
             call IO_error(500,ext_msg=homogThermal%get_asString('type'))
@@ -485,7 +484,7 @@ subroutine material_parseHomogenization
     endif
   enddo
 
-end subroutine material_parseHomogenization
+end subroutine parseHomogenization
 
 
 end module homogenization

src/homogenization_mechanical.f90

@@ -7,15 +7,13 @@ submodule(homogenization) mechanical
 
   interface
 
-    module subroutine pass_init
+    module subroutine pass_init()
     end subroutine pass_init
 
-    module subroutine isostrain_init
+    module subroutine isostrain_init()
     end subroutine isostrain_init
 
-    module subroutine RGC_init(num_homogMech)
+    module subroutine RGC_init()
-      class(tNode), pointer, intent(in) :: &
-        num_homogMech                                                                               !< pointer to mechanical homogenization numerics data
     end subroutine RGC_init
 
 
@@ -52,6 +50,12 @@ submodule(homogenization) mechanical
 
   end interface
 
+  type :: tOutput                                                                                   !< requested output (per phase)
+    character(len=pStringLen), allocatable, dimension(:) :: &
+      label
+  end type tOutput
+  type(tOutput), allocatable, dimension(:) :: output_mechanical
+
   integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:),   allocatable :: &
     homogenization_type                                                                             !< type of each homogenization
 
@@ -60,27 +64,20 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief Allocate variables and set parameters.
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_init(num_homog)
+module subroutine mechanical_init()
-
-  class(tNode), pointer, intent(in) :: &
-    num_homog
-
-  class(tNode), pointer :: &
-    num_homogMech
 
   print'(/,1x,a)', '<<<+-  homogenization:mechanical init  -+>>>'
 
-  call material_parseHomogenization2()
+  call parseMechanical()
 
-  allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
+  allocate(homogenization_dPdF(3,3,3,3,discretization_Ncells), source=0.0_pReal)
-  homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems)                   ! initialize to identity
+  homogenization_F0 = spread(math_I3,3,discretization_Ncells)
-  homogenization_F = homogenization_F0                                                              ! initialize to identity
+  homogenization_F = homogenization_F0
-  allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems),        source=0.0_pReal)
+  allocate(homogenization_P(3,3,discretization_Ncells),source=0.0_pReal)
 
-  num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
+  if (any(homogenization_type == HOMOGENIZATION_NONE_ID))      call pass_init()
-  if (any(homogenization_type == HOMOGENIZATION_NONE_ID))      call pass_init
+  if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call isostrain_init()
-  if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call isostrain_init
+  if (any(homogenization_type == HOMOGENIZATION_RGC_ID))       call RGC_init()
-  if (any(homogenization_type == HOMOGENIZATION_RGC_ID))       call RGC_init(num_homogMech)
 
 end subroutine mechanical_init
 
@@ -185,8 +182,10 @@ module subroutine mechanical_results(group_base,ho)
   character(len=*), intent(in) :: group_base
   integer, intent(in)          :: ho
 
+  integer :: ou
   character(len=:), allocatable :: group
 
+
   group = trim(group_base)//'/mechanical'
   call results_closeGroup(results_addGroup(group))
 
@@ -197,12 +196,17 @@ module subroutine mechanical_results(group_base,ho)
 
   end select
 
-  !temp = reshape(homogenization_F,[3,3,discretization_nIPs*discretization_Nelems])
+  do ou = 1, size(output_mechanical(1)%label)
-  !call results_writeDataset(group,temp,'F',&
+
-  !                          'deformation gradient','1')
+    select case (output_mechanical(ho)%label(ou))
-  !temp = reshape(homogenization_P,[3,3,discretization_nIPs*discretization_Nelems])
+      case('F')
-  !call results_writeDataset(group,temp,'P',&
+        call results_writeDataset(reshape(homogenization_F,[3,3,discretization_nCells]),group,'F', &
-  !                          '1st Piola-Kirchhoff stress','Pa')
+                                  'deformation gradient','1')
+      case('P')
+        call results_writeDataset(reshape(homogenization_P,[3,3,discretization_nCells]),group,'P', &
+                                  'deformation gradient','1')
+    end select
+  end do
 
 end subroutine mechanical_results
 
@@ -210,35 +214,42 @@ end subroutine mechanical_results
 !--------------------------------------------------------------------------------------------------
 !> @brief parses the homogenization part from the material configuration
 !--------------------------------------------------------------------------------------------------
-subroutine material_parseHomogenization2()
+subroutine parseMechanical()
 
   class(tNode), pointer :: &
     material_homogenization, &
     homog, &
-    homogMech
+    mechanical
+
+  integer :: ho
 
-  integer :: h
 
   material_homogenization => config_material%get('homogenization')
 
   allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
+  allocate(output_mechanical(size(material_name_homogenization)))
 
-  do h=1, size(material_name_homogenization)
+  do ho=1, size(material_name_homogenization)
-    homog => material_homogenization%get(h)
+    homog => material_homogenization%get(ho)
-    homogMech => homog%get('mechanical')
+    mechanical => homog%get('mechanical')
-    select case (homogMech%get_asString('type'))
+#if defined(__GFORTRAN__)
+    output_mechanical(ho)%label = output_as1dString(mechanical)
+#else
+    output_mechanical(ho)%label = mechanical%get_as1dString('output',defaultVal=emptyStringArray)
+#endif
+    select case (mechanical%get_asString('type'))
       case('pass')
-        homogenization_type(h) = HOMOGENIZATION_NONE_ID
+        homogenization_type(ho) = HOMOGENIZATION_NONE_ID
       case('isostrain')
-        homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
+        homogenization_type(ho) = HOMOGENIZATION_ISOSTRAIN_ID
       case('RGC')
-        homogenization_type(h) = HOMOGENIZATION_RGC_ID
+        homogenization_type(ho) = HOMOGENIZATION_RGC_ID
       case default
-        call IO_error(500,ext_msg=homogMech%get_asString('type'))
+        call IO_error(500,ext_msg=mechanical%get_asString('type'))
     end select
   end do
 
-end subroutine material_parseHomogenization2
+end subroutine parseMechanical
 
 
 end submodule mechanical

src/homogenization_mechanical_RGC.f90

@@ -71,10 +71,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief allocates all necessary fields, reads information from material configuration file
 !--------------------------------------------------------------------------------------------------
-module subroutine RGC_init(num_homogMech)
+module subroutine RGC_init()
-
-  class(tNode), pointer, intent(in) :: &
-    num_homogMech                                                                                   !< pointer to mechanical homogenization numerics data
 
   integer :: &
     ho, &
@@ -82,6 +79,8 @@ module subroutine RGC_init(num_homogMech)
     sizeState, nIntFaceTot
 
   class (tNode), pointer :: &
+    num_homogenization, &
+    num_mechanical, &
     num_RGC, &                                                                                      ! pointer to RGC numerics data
     material_homogenization, &
     homog, &
@@ -105,7 +104,9 @@ module subroutine RGC_init(num_homogMech)
   allocate(state0(material_homogenization%length))
   allocate(dependentState(material_homogenization%length))
 
-  num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict)
+  num_homogenization => config_numerics%get('homogenization',defaultVal=emptyDict)
+  num_mechanical => num_homogenization%get('mechanical',defaultVal=emptyDict)
+  num_RGC => num_mechanical%get('RGC',defaultVal=emptyDict)
 
   num%atol         =  num_RGC%get_asFloat('atol',              defaultVal=1.0e+4_pReal)
   num%rtol         =  num_RGC%get_asFloat('rtol',              defaultVal=1.0e-3_pReal)
@@ -171,8 +172,8 @@ module subroutine RGC_init(num_homogMech)
     allocate(homogState(ho)%state0   (sizeState,Nmembers), source=0.0_pReal)
     allocate(homogState(ho)%state    (sizeState,Nmembers), source=0.0_pReal)
 
-    stt%relaxationVector   => homogState(ho)%state(1:nIntFaceTot,:)
+    stt%relaxationVector => homogState(ho)%state(1:nIntFaceTot,:)
-    st0%relaxationVector   => homogState(ho)%state0(1:nIntFaceTot,:)
+    st0%relaxationVector => homogState(ho)%state0(1:nIntFaceTot,:)
 
     allocate(dst%volumeDiscrepancy(   Nmembers), source=0.0_pReal)
     allocate(dst%relaxationRate_avg(  Nmembers), source=0.0_pReal)

src/homogenization_thermal.f90

@@ -52,10 +52,11 @@ module subroutine thermal_init()
   allocate(current(configHomogenizations%length))
 
   do ho = 1, configHomogenizations%length
-    allocate(current(ho)%T(count(material_homogenizationID==ho)), source=300.0_pReal)
+    allocate(current(ho)%T(count(material_homogenizationID==ho)), source=T_ROOM)
     allocate(current(ho)%dot_T(count(material_homogenizationID==ho)), source=0.0_pReal)
     configHomogenization => configHomogenizations%get(ho)
     associate(prm => param(ho))
+
       if (configHomogenization%contains('thermal')) then
         configHomogenizationThermal => configHomogenization%get('thermal')
 #if defined (__GFORTRAN__)
@@ -63,13 +64,22 @@ module subroutine thermal_init()
 #else
         prm%output = configHomogenizationThermal%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
+        select case (configHomogenizationThermal%get_asString('type'))
+
+          case ('pass')
+            call pass_init()
+
+          case ('isothermal')
+            call isotemperature_init()
+
+        end select
       else
         prm%output = emptyStringArray
       end if
+
     end associate
   end do
 
-  call pass_init()
 
 end subroutine thermal_init
 

src/homogenization_thermal_isotemperature.f90

@@ -1,13 +1,14 @@
 !--------------------------------------------------------------------------------------------------
 !> @author Martin Diehl, KU Leuven
-!> @brief Dummy homogenization scheme for 1 constituent per material point
+!> @brief Isotemperature homogenization
 !--------------------------------------------------------------------------------------------------
 submodule(homogenization:thermal) isotemperature
 
 contains
 
-module subroutine isotemperature_init
+module subroutine isotemperature_init()
 
+  print'(/,1x,a)', '<<<+-  homogenization:thermal:isotemperature init  -+>>>'
 
 end subroutine isotemperature_init
 

src/homogenization_thermal_pass.f90

@@ -10,6 +10,9 @@ module subroutine pass_init()
 
   print'(/,1x,a)', '<<<+-  homogenization:thermal:pass init  -+>>>'
 
+  if (homogenization_Nconstituents(1) /= 1) &
+    call IO_error(211,ext_msg='N_constituents (pass)')
+
 end subroutine pass_init
 
 end submodule thermal_pass

src/phase.f90

@@ -91,6 +91,11 @@ module phase
       integer,          intent(in) :: ph
     end subroutine damage_results
 
+    module subroutine thermal_results(group,ph)
+      character(len=*), intent(in) :: group
+      integer,          intent(in) :: ph
+    end subroutine thermal_results
+
     module subroutine mechanical_forward()
     end subroutine mechanical_forward
 
@@ -487,6 +492,7 @@ subroutine phase_results()
 
     call mechanical_results(group,ph)
     call damage_results(group,ph)
+    call thermal_results(group,ph)
 
   end do
 

src/phase_damage.f90

@@ -336,7 +336,7 @@ end subroutine damage_results
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief contains the constitutive equation for calculating the rate of change of microstructure
+!> @brief Constitutive equation for calculating the rate of change of microstructure.
 !--------------------------------------------------------------------------------------------------
 function phase_damage_collectDotState(ph,en) result(broken)
 

src/phase_mechanical.f90

@@ -180,11 +180,12 @@ submodule(phase) mechanical
     end function elastic_nu
 
   end interface
-  type :: tOutput                                                                                   !< new requested output (per phase)
+
+  type :: tOutput                                                                                   !< requested output (per phase)
     character(len=pStringLen), allocatable, dimension(:) :: &
       label
   end type tOutput
-  type(tOutput), allocatable, dimension(:) :: output_constituent
+  type(tOutput), allocatable, dimension(:) :: output_mechanical
 
   procedure(integrateStateFPI), pointer :: integrateState
 
@@ -216,7 +217,7 @@ module subroutine mechanical_init(phases)
   print'(/,1x,a)', '<<<+-  phase:mechanical init  -+>>>'
 
 !-------------------------------------------------------------------------------------------------
-  allocate(output_constituent(phases%length))
+  allocate(output_mechanical(phases%length))
 
   allocate(phase_mechanical_Fe(phases%length))
   allocate(phase_mechanical_Fi(phases%length))
@@ -251,12 +252,12 @@ module subroutine mechanical_init(phases)
     allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
     allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
 
-    phase   => phases%get(ph)
+    phase => phases%get(ph)
-    mech    => phase%get('mechanical')
+    mech  => phase%get('mechanical')
 #if defined(__GFORTRAN__)
-    output_constituent(ph)%label = output_as1dString(mech)
+    output_mechanical(ph)%label = output_as1dString(mech)
 #else
-    output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
+    output_mechanical(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
   enddo
 
@@ -330,7 +331,7 @@ module subroutine mechanical_results(group,ph)
   integer,          intent(in) :: ph
 
 
-  call crystallite_results(group,ph)
+  call results(group,ph)
 
   select case(phase_plasticity(ph))
 
@@ -879,9 +880,9 @@ end function integrateStateRK
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief writes crystallite results to HDF5 output file
+!> @brief Write mechanical results to HDF5 output file.
 !--------------------------------------------------------------------------------------------------
-subroutine crystallite_results(group,ph)
+subroutine results(group,ph)
 
   character(len=*), intent(in) :: group
   integer,          intent(in) :: ph
@@ -891,9 +892,9 @@ subroutine crystallite_results(group,ph)
 
   call results_closeGroup(results_addGroup(group//'/mechanical'))
 
-  do ou = 1, size(output_constituent(ph)%label)
+  do ou = 1, size(output_mechanical(ph)%label)
 
-    select case (output_constituent(ph)%label(ou))
+    select case (output_mechanical(ph)%label(ou))
       case('F')
         call results_writeDataset(phase_mechanical_F(ph)%data,group//'/mechanical/','F',&
                                  'deformation gradient','1')
@@ -919,13 +920,13 @@ subroutine crystallite_results(group,ph)
         call results_writeDataset(phase_mechanical_S(ph)%data,group//'/mechanical/','S', &
                                  'second Piola-Kirchhoff stress','Pa')
       case('O')
-        call results_writeDataset(to_quaternion(phase_O(ph)%data),group//'/mechanical',output_constituent(ph)%label(ou),&
+        call results_writeDataset(to_quaternion(phase_O(ph)%data),group//'/mechanical',output_mechanical(ph)%label(ou),&
                                  'crystal orientation as quaternion','q_0 (q_1 q_2 q_3)')
-        call results_addAttribute('lattice',phase_lattice(ph),group//'/mechanical/'//output_constituent(ph)%label(ou))
+        call results_addAttribute('lattice',phase_lattice(ph),group//'/mechanical/'//output_mechanical(ph)%label(ou))
         if (any(phase_lattice(ph) == ['hP', 'tI'])) &
-          call results_addAttribute('c/a',phase_cOverA(ph),group//'/mechanical/'//output_constituent(ph)%label(ou))
+          call results_addAttribute('c/a',phase_cOverA(ph),group//'/mechanical/'//output_mechanical(ph)%label(ou))
     end select
-  enddo
+  end do
 
 
   contains
@@ -947,7 +948,7 @@ subroutine crystallite_results(group,ph)
 
  end function to_quaternion
 
-end subroutine crystallite_results
+end subroutine results
 
 
 !--------------------------------------------------------------------------------------------------
@@ -1335,5 +1336,4 @@ module subroutine phase_set_F(F,co,ce)
 
 end subroutine phase_set_F
 
-
 end submodule mechanical

src/phase_thermal.f90

@@ -6,6 +6,7 @@ submodule(phase) thermal
   type :: tThermalParameters
     real(pReal) ::                 C_p = 0.0_pReal                                                  !< heat capacity
     real(pReal), dimension(3,3) :: K   = 0.0_pReal                                                  !< thermal conductivity
+    character(len=pStringLen), allocatable, dimension(:) :: output
   end type tThermalParameters
 
   integer, dimension(:), allocatable :: &
@@ -96,7 +97,7 @@ module subroutine thermal_init(phases)
 
   do ph = 1, phases%length
     Nmembers = count(material_phaseID == ph)
-    allocate(current(ph)%T(Nmembers),source=300.0_pReal)
+    allocate(current(ph)%T(Nmembers),source=T_ROOM)
     allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
     phase => phases%get(ph)
     thermal => phase%get('thermal',defaultVal=emptyDict)
@@ -108,6 +109,11 @@ module subroutine thermal_init(phases)
       if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asFloat('K_33')
       param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
 
+#if defined(__GFORTRAN__)
+      param(ph)%output = output_as1dString(thermal)
+#else
+      param(ph)%output = thermal%get_as1dString('output',defaultVal=emptyStringArray)
+#endif
       sources => thermal%get('source',defaultVal=emptyList)
       thermal_Nsources(ph) = sources%length
     else
@@ -381,4 +387,35 @@ function thermal_active(source_label,src_length)  result(active_source)
 end function thermal_active
 
 
+!----------------------------------------------------------------------------------------------
+!< @brief writes damage sources results to HDF5 output file
+!----------------------------------------------------------------------------------------------
+module subroutine thermal_results(group,ph)
+
+  character(len=*), intent(in) :: group
+  integer,          intent(in) :: ph
+
+
+  integer :: ou
+
+  if (allocated(param(ph)%output)) then
+    call results_closeGroup(results_addGroup(group//'thermal'))
+  else
+    return
+  endif
+
+  do ou = 1, size(param(ph)%output)
+
+    select case(trim(param(ph)%output(ou)))
+
+      case ('T')
+        call results_writeDataset(current(ph)%T,group//'thermal','T', 'temperature','T')
+
+    end select
+
+  end do
+
+end subroutine thermal_results
+
+
 end submodule thermal