Merge branch 'Fortran-cleaning' into 'development'

Fortran cleaning

See merge request damask/DAMASK!365

src/IO.f90

@@ -604,7 +604,7 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
       msg = 'read only the first document'
     case default
       msg = 'unknown warning number'
-    end select
+  end select
 
   !$OMP CRITICAL (write2out)
   write(IO_STDERR,'(/,a)')                ' ┌'//IO_DIVIDER//'┐'
@@ -658,7 +658,7 @@ subroutine selfTest
   if(any([1,1,1]     /= IO_stringPos('a')))         error stop 'IO_stringPos'
   if(any([2,2,3,5,5] /= IO_stringPos(' aa b')))     error stop 'IO_stringPos'
 
-  str=' 1.0 xxx'
+  str = ' 1.0 xxx'
   chunkPos = IO_stringPos(str)
   if(dNeq(1.0_pReal,IO_floatValue(str,chunkPos,1))) error stop 'IO_floatValue'
 

src/YAML_types.f90

@@ -266,7 +266,7 @@ subroutine selfTest
         if(any(dNeq(l2%get_as1dFloat(1),[2.0_pReal,3.0_pReal]))) error stop 'byIndex_as1dFloat'
         call l2%append(l3)
         x = l2%as2dFloat()
-        if(x(2,1)/= 4.0_pReal)                                   error stop 'byKey_as2dFloat'
+        if(dNeq(x(2,1),4.0_pReal))                               error stop 'byKey_as2dFloat'
         if(any(dNeq(pack(l2%as2dFloat(),.true.),&
                [2.0_pReal,4.0_pReal,3.0_pReal,5.0_pReal])))      error stop 'byKey_as2dFloat'
         n => l2

src/grid/grid_damage_spectral.f90

@@ -196,7 +196,7 @@ function grid_damage_spectral_solution(timeinc) result(solution)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce)
+    call homogenization_set_phi(phi_current(i,j,k),ce)
   enddo; enddo; enddo
 
   call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
@@ -233,7 +233,7 @@ subroutine grid_damage_spectral_forward(cutBack)
     call DMDAVecRestoreArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)
     do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
       ce = ce + 1
-      call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),ce)
+      call homogenization_set_phi(phi_current(i,j,k),ce)
     enddo; enddo; enddo
   else
     phi_lastInc = phi_current
@@ -259,7 +259,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   PetscObject :: dummy
   PetscErrorCode :: ierr
   integer :: i, j, k, ce
-  real(pReal)   :: phiDot, dPhiDot_dPhi, mobility
+  real(pReal)   :: phiDot, mobility
 
   phi_current = x_scal
 !--------------------------------------------------------------------------------------------------
@@ -281,7 +281,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   ce = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
-    call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi_current(i,j,k),ce)
+    call damage_nonlocal_getSourceAndItsTangent(phiDot, phi_current(i,j,k),ce)
     mobility = damage_nonlocal_getMobility(ce)
     scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + phiDot) &
                             + mobility*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &

src/grid/grid_thermal_spectral.f90

@@ -282,9 +282,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
     ce = ce + 1
     call thermal_conduction_getSource(Tdot,1,ce)
     scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) &
-                            + thermal_conduction_getMassDensity (ce)* &
+                            + homogenization_thermal_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
-                              thermal_conduction_getSpecificHeat(ce)*(T_lastInc(i,j,k)  - &
-                                                                          T_current(i,j,k))&
                             + mu_ref*T_current(i,j,k)
   enddo; enddo; enddo
 
@@ -314,7 +312,7 @@ subroutine updateReference
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     ce = ce + 1
     K_ref  = K_ref  + thermal_conduction_getConductivity(ce)
-    mu_ref = mu_ref + thermal_conduction_getMassDensity(ce)* thermal_conduction_getSpecificHeat(ce)
+    mu_ref = mu_ref + homogenization_thermal_mu_T(ce)
   enddo; enddo; enddo
   K_ref = K_ref*wgt
   call MPI_Allreduce(MPI_IN_PLACE,K_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)

src/homogenization.f90

@@ -39,8 +39,6 @@ module homogenization
     thermal_type                                                                                    !< thermal transport model
   integer(kind(DAMAGE_none_ID)),              dimension(:),   allocatable :: &
     damage_type                                                                                     !< nonlocal damage model
-  integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:),   allocatable :: &
-    homogenization_type                                                                             !< type of each homogenization
 
   type, private :: tNumerics_damage
     real(pReal) :: &
@@ -117,6 +115,16 @@ module homogenization
       integer, intent(in)          :: ho
     end subroutine mechanical_results
 
+    module subroutine damage_results(ho,group)
+      integer,          intent(in) :: ho
+      character(len=*), intent(in) :: group
+    end subroutine damage_results
+
+    module subroutine thermal_results(ho,group)
+      integer,          intent(in) :: ho
+      character(len=*), intent(in) :: group
+    end subroutine thermal_results
+
     module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy)
       real(pReal), intent(in) :: &
         subdt                                                                                       !< current time step
@@ -133,26 +141,16 @@ module homogenization
       real(pReal), dimension(3,3) :: K
     end function thermal_conduction_getConductivity
 
-    module function thermal_conduction_getSpecificHeat(ce) result(c_P)
+    module function homogenization_thermal_mu_T(ce) result(mu_T)
       integer, intent(in) :: ce
-      real(pReal) :: c_P
+      real(pReal) :: mu_T
-    end function thermal_conduction_getSpecificHeat
+    end function homogenization_thermal_mu_T
-
-    module function thermal_conduction_getMassDensity(ce) result(rho)
-      integer, intent(in) :: ce
-      real(pReal) :: rho
-    end function thermal_conduction_getMassDensity
 
     module subroutine homogenization_thermal_setField(T,dot_T, ce)
       integer, intent(in) :: ce
       real(pReal),   intent(in) :: T, dot_T
     end subroutine homogenization_thermal_setField
 
-    module subroutine thermal_conduction_results(ho,group)
-      integer,          intent(in) :: ho
-      character(len=*), intent(in) :: group
-    end subroutine thermal_conduction_results
-
     module function homogenization_thermal_T(ce) result(T)
       integer, intent(in) :: ce
       real(pReal) :: T
@@ -170,37 +168,31 @@ module homogenization
       real(pReal) :: M
     end function damage_nonlocal_getMobility
 
-    module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce)
+    module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce)
       integer, intent(in) :: ce
-      real(pReal),   intent(in) :: &
+      real(pReal), intent(in) :: &
         phi
-      real(pReal) :: &
+      real(pReal), intent(out) :: &
-        phiDot, dPhiDot_dPhi
+        phiDot
     end subroutine damage_nonlocal_getSourceAndItsTangent
 
-    module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
+    module subroutine homogenization_set_phi(phi,ce)
       integer, intent(in) :: ce
       real(pReal),   intent(in) :: &
         phi
-    end subroutine damage_nonlocal_putNonLocalDamage
+    end subroutine homogenization_set_phi
-
-    module subroutine damage_nonlocal_results(ho,group)
-      integer,          intent(in) :: ho
-      character(len=*), intent(in) :: group
-    end subroutine damage_nonlocal_results
 
   end interface
 
   public ::  &
     homogenization_init, &
     materialpoint_stressAndItsTangent, &
-    thermal_conduction_getSpecificHeat, &
+    homogenization_thermal_mu_T, &
     thermal_conduction_getConductivity, &
-    thermal_conduction_getMassDensity, &
     thermal_conduction_getSource, &
     damage_nonlocal_getMobility, &
     damage_nonlocal_getSourceAndItsTangent, &
-    damage_nonlocal_putNonLocalDamage, &
+    homogenization_set_phi, &
     homogenization_thermal_setfield, &
     homogenization_thermal_T, &
     homogenization_forward, &
@@ -211,7 +203,6 @@ module homogenization
     DAMAGE_NONLOCAL_ID
 
   public :: &
-    damage_nonlocal_init, &
     damage_nonlocal_getDiffusion
 
 contains
@@ -242,8 +233,6 @@ subroutine homogenization_init()
   call mechanical_init(num_homog)
   call thermal_init()
   call damage_init()
-  call damage_nonlocal_init()
-
 
 end subroutine homogenization_init
 
@@ -371,14 +360,14 @@ subroutine homogenization_results
       case(DAMAGE_NONLOCAL_ID)
         group = trim(group_base)//'/damage'
         call results_closeGroup(results_addGroup(group))
-        call damage_nonlocal_results(ho,group)
+        call damage_results(ho,group)
     end select
 
     select case(thermal_type(ho))
       case(THERMAL_CONDUCTION_ID)
         group = trim(group_base)//'/thermal'
         call results_closeGroup(results_addGroup(group))
-        call thermal_conduction_results(ho,group)
+        call thermal_results(ho,group)
     end select
 
  enddo
@@ -458,41 +447,6 @@ subroutine homogenization_restartRead(fileHandle)
 end subroutine homogenization_restartRead
 
 
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-subroutine damage_nonlocal_init
-
-  integer :: Ninstances,Nmaterialpoints,h
-  class(tNode), pointer :: &
-    num_generic, &
-    material_homogenization
-
-  print'(/,a)', ' <<<+-  damage_nonlocal init  -+>>>'; flush(6)
-
-!------------------------------------------------------------------------------------
-! read numerics parameter
-  num_generic => config_numerics%get('generic',defaultVal= emptyDict)
-  num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
-
-  Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
-
-  material_homogenization => config_material%get('homogenization')
-  do h = 1, material_homogenization%length
-    if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
-
-    Nmaterialpoints = count(material_homogenizationAt == h)
-    damageState_h(h)%sizeState = 1
-    allocate(damageState_h(h)%state0   (1,Nmaterialpoints), source=1.0_pReal)
-    allocate(damageState_h(h)%state    (1,Nmaterialpoints), source=1.0_pReal)
-
-  enddo
-
-end subroutine damage_nonlocal_init
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief returns homogenized non local damage diffusion tensor in reference configuration
 !--------------------------------------------------------------------------------------------------
@@ -521,14 +475,12 @@ end function damage_nonlocal_getDiffusion
 
 !--------------------------------------------------------------------------------------------------
 !> @brief parses the homogenization part from the material configuration
-! ToDo: This should be done in homogenization
 !--------------------------------------------------------------------------------------------------
 subroutine material_parseHomogenization
 
   class(tNode), pointer :: &
     material_homogenization, &
     homog, &
-    homogMech, &
     homogThermal, &
     homogDamage
 
@@ -536,23 +488,11 @@ subroutine material_parseHomogenization
 
   material_homogenization => config_material%get('homogenization')
 
-  allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
   allocate(thermal_type(size(material_name_homogenization)),        source=THERMAL_isothermal_ID)
   allocate(damage_type (size(material_name_homogenization)),        source=DAMAGE_none_ID)
 
   do h=1, size(material_name_homogenization)
     homog => material_homogenization%get(h)
-    homogMech => homog%get('mechanical')
-    select case (homogMech%get_asString('type'))
-      case('pass')
-        homogenization_type(h) = HOMOGENIZATION_NONE_ID
-      case('isostrain')
-        homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
-      case('RGC')
-        homogenization_type(h) = HOMOGENIZATION_RGC_ID
-      case default
-        call IO_error(500,ext_msg=homogMech%get_asString('type'))
-    end select
 
     if (homog%contains('thermal')) then
       homogThermal => homog%get('thermal')

src/homogenization_damage.f90

@@ -37,12 +37,15 @@ module subroutine damage_init()
   class(tNode), pointer :: &
     configHomogenizations, &
     configHomogenization, &
-    configHomogenizationDamage
+    configHomogenizationDamage, &
+    num_generic, &
+    material_homogenization
   integer :: ho
+  integer :: Ninstances,Nmaterialpoints,h
 
 
   print'(/,a)', ' <<<+-  homogenization:damage init  -+>>>'
-  print'(/,a)', ' <<<+-  homogenization:damage:isodamage init  -+>>>'
+  print'(/,a)', ' <<<+-  homogenization:damage:pass init  -+>>>'
 
   configHomogenizations => config_material%get('homogenization')
   allocate(param(configHomogenizations%length))
@@ -65,6 +68,24 @@ module subroutine damage_init()
     end associate
   enddo
 
+!------------------------------------------------------------------------------------
+! read numerics parameter
+  num_generic => config_numerics%get('generic',defaultVal= emptyDict)
+  num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
+
+  Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
+
+  material_homogenization => config_material%get('homogenization')
+  do h = 1, material_homogenization%length
+    if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle
+
+    Nmaterialpoints = count(material_homogenizationAt == h)
+    damageState_h(h)%sizeState = 1
+    allocate(damageState_h(h)%state0   (1,Nmaterialpoints), source=1.0_pReal)
+    allocate(damageState_h(h)%state    (1,Nmaterialpoints), source=1.0_pReal)
+
+  enddo
+
 end subroutine damage_init
 
 
@@ -76,14 +97,10 @@ module subroutine damage_partition(ce)
   real(pReal) :: phi
   integer,     intent(in) :: ce
 
-  integer :: co
-
 
   if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return
   phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce))
-  do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
+  call phase_damage_set_phi(phi,1,ce)
-    call phase_damage_set_phi(phi,co,ce)
-  enddo
 
 end subroutine damage_partition
 
@@ -95,17 +112,9 @@ end subroutine damage_partition
 module function damage_nonlocal_getMobility(ce) result(M)
 
   integer, intent(in) :: ce
-  integer :: &
-    co
   real(pReal) :: M
 
-  M = 0.0_pReal
+  M = lattice_M(material_phaseID(1,ce))
-
-  do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
-    M = M + lattice_M(material_phaseID(co,ce))
-  enddo
-
-  M = M/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
 end function damage_nonlocal_getMobility
 
@@ -113,20 +122,15 @@ end function damage_nonlocal_getMobility
 !--------------------------------------------------------------------------------------------------
 !> @brief  calculates homogenized damage driving forces
 !--------------------------------------------------------------------------------------------------
-module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ce)
+module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, phi, ce)
 
   integer, intent(in) :: ce
-  real(pReal),   intent(in) :: &
+  real(pReal), intent(in) :: &
     phi
-  real(pReal) :: &
+  real(pReal), intent(out) :: &
-    phiDot, dPhiDot_dPhi
+    phiDot
 
-  phiDot = 0.0_pReal
+  phiDot = phase_damage_phi_dot(phi, 1, ce)
-  dPhiDot_dPhi = 0.0_pReal
-
-  call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
-  phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
-  dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
 end subroutine damage_nonlocal_getSourceAndItsTangent
 
@@ -134,7 +138,7 @@ end subroutine damage_nonlocal_getSourceAndItsTangent
 !--------------------------------------------------------------------------------------------------
 !> @brief updated nonlocal damage field with solution from damage phase field PDE
 !--------------------------------------------------------------------------------------------------
-module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
+module subroutine homogenization_set_phi(phi,ce)
 
   integer, intent(in) :: ce
   real(pReal),   intent(in) :: &
@@ -146,14 +150,15 @@ module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
   ho = material_homogenizationID(ce)
   en = material_homogenizationEntry(ce)
   damagestate_h(ho)%state(1,en) = phi
+  current(ho)%phi(en) = phi
 
-end subroutine damage_nonlocal_putNonLocalDamage
+end subroutine homogenization_set_phi
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-module subroutine damage_nonlocal_results(ho,group)
+module subroutine damage_results(ho,group)
 
   integer,          intent(in) :: ho
   character(len=*), intent(in) :: group
@@ -170,6 +175,6 @@ module subroutine damage_nonlocal_results(ho,group)
   enddo outputsLoop
   end associate
 
-end subroutine damage_nonlocal_results
+end subroutine damage_results
 
 end submodule damage

src/homogenization_mechanical.f90

@@ -71,6 +71,9 @@ submodule(homogenization) mechanical
 
   end interface
 
+  integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:),   allocatable :: &
+    homogenization_type                                                                             !< type of each homogenization
+
 contains
 
 !--------------------------------------------------------------------------------------------------
@@ -86,6 +89,8 @@ module subroutine mechanical_init(num_homog)
 
   print'(/,a)', ' <<<+-  homogenization:mechanical init  -+>>>'
 
+  call material_parseHomogenization2()
+
   allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
   homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems)                   ! initialize to identity
   homogenization_F = homogenization_F0                                                              ! initialize to identity
@@ -127,7 +132,7 @@ module subroutine mechanical_partition(subF,ce)
   end select chosenHomogenization
 
   do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
-    call phase_mechanical_setF(Fs(1:3,1:3,co),co,ce)
+    call phase_set_F(Fs(1:3,1:3,co),co,ce)
   enddo
 
 
@@ -150,13 +155,13 @@ module subroutine mechanical_homogenize(dt,ce)
   chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce)))
 
     case (HOMOGENIZATION_NONE_ID) chosenHomogenization
-        homogenization_P(1:3,1:3,ce)            = phase_mechanical_getP(1,ce)
+        homogenization_P(1:3,1:3,ce)            = phase_P(1,ce)
         homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce)
 
     case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
       do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
         dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
-        Ps(:,:,co) = phase_mechanical_getP(co,ce)
+        Ps(:,:,co) = phase_P(co,ce)
       enddo
       call isostrain_averageStressAndItsTangent(&
         homogenization_P(1:3,1:3,ce), &
@@ -167,7 +172,7 @@ module subroutine mechanical_homogenize(dt,ce)
     case (HOMOGENIZATION_RGC_ID) chosenHomogenization
       do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
         dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
-        Ps(:,:,co) = phase_mechanical_getP(co,ce)
+        Ps(:,:,co) = phase_P(co,ce)
       enddo
       call RGC_averageStressAndItsTangent(&
         homogenization_P(1:3,1:3,ce), &
@@ -203,8 +208,8 @@ module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy)
   if (homogenization_type(material_homogenizationID(ce)) == HOMOGENIZATION_RGC_ID) then
       do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
         dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce)
-        Fs(:,:,co)        = phase_mechanical_getF(co,ce)
+        Fs(:,:,co)        = phase_F(co,ce)
-        Ps(:,:,co)        = phase_mechanical_getP(co,ce)
+        Ps(:,:,co)        = phase_P(co,ce)
       enddo
       doneAndHappy = RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce)
   else
@@ -244,4 +249,38 @@ module subroutine mechanical_results(group_base,ho)
 end subroutine mechanical_results
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief parses the homogenization part from the material configuration
+!--------------------------------------------------------------------------------------------------
+subroutine material_parseHomogenization2()
+
+  class(tNode), pointer :: &
+    material_homogenization, &
+    homog, &
+    homogMech
+
+  integer :: h
+
+  material_homogenization => config_material%get('homogenization')
+
+  allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
+
+  do h=1, size(material_name_homogenization)
+    homog => material_homogenization%get(h)
+    homogMech => homog%get('mechanical')
+    select case (homogMech%get_asString('type'))
+      case('pass')
+        homogenization_type(h) = HOMOGENIZATION_NONE_ID
+      case('isostrain')
+        homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
+      case('RGC')
+        homogenization_type(h) = HOMOGENIZATION_RGC_ID
+      case default
+        call IO_error(500,ext_msg=homogMech%get_asString('type'))
+    end select
+  enddo
+
+end subroutine material_parseHomogenization2
+
+
 end submodule mechanical

src/homogenization_thermal.f90

@@ -45,7 +45,7 @@ module subroutine thermal_init()
 
 
   print'(/,a)', ' <<<+-  homogenization:thermal init  -+>>>'
-  print'(/,a)', ' <<<+-  homogenization:thermal:isotemperature   init  -+>>>'
+  print'(/,a)', ' <<<+-  homogenization:thermal:pass init  -+>>>'
 
 
 
@@ -128,10 +128,20 @@ module function thermal_conduction_getConductivity(ce) result(K)
 end function thermal_conduction_getConductivity
 
 
+module function homogenization_thermal_mu_T(ce) result(mu_T)
+  
+  integer, intent(in) :: ce
+  real(pReal) :: mu_T
+
+  mu_T = c_P(ce) * rho(ce)
+
+end function homogenization_thermal_mu_T
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief returns homogenized specific heat capacity
 !--------------------------------------------------------------------------------------------------
-module function thermal_conduction_getSpecificHeat(ce) result(c_P)
+function c_P(ce)
 
   integer, intent(in) :: ce
   real(pReal) :: c_P
@@ -139,21 +149,20 @@ module function thermal_conduction_getSpecificHeat(ce) result(c_P)
   integer :: co
 
 
-  c_P = 0.0_pReal
+  c_P = lattice_c_p(material_phaseID(1,ce))
-
+  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-  do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
     c_P = c_P + lattice_c_p(material_phaseID(co,ce))
   enddo
 
   c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
-end function thermal_conduction_getSpecificHeat
+end function c_P
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief returns homogenized mass density
 !--------------------------------------------------------------------------------------------------
-module function thermal_conduction_getMassDensity(ce) result(rho)
+function rho(ce)
 
   integer, intent(in) :: ce
   real(pReal) :: rho
@@ -161,15 +170,14 @@ module function thermal_conduction_getMassDensity(ce) result(rho)
   integer :: co
 
 
-  rho = 0.0_pReal
+  rho = lattice_rho(material_phaseID(1,ce))
-
+  do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
-  do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
     rho = rho + lattice_rho(material_phaseID(co,ce))
   enddo
 
   rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
 
-end function thermal_conduction_getMassDensity
+end function rho
 
 
 
@@ -193,7 +201,7 @@ end subroutine homogenization_thermal_setField
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------
-module subroutine thermal_conduction_results(ho,group)
+module subroutine thermal_results(ho,group)
 
   integer,          intent(in) :: ho
   character(len=*), intent(in) :: group
@@ -209,7 +217,7 @@ module subroutine thermal_conduction_results(ho,group)
     enddo outputsLoop
   end associate
 
-end subroutine thermal_conduction_results
+end subroutine thermal_results
 
 
 module function homogenization_thermal_T(ce) result(T)

src/phase.f90

@@ -145,20 +145,20 @@ module phase
       real(pReal), dimension(3,3) :: L_p
     end function mechanical_L_p
 
-    module function phase_mechanical_getF(co,ce) result(F)
+    module function phase_F(co,ce) result(F)
       integer, intent(in) :: co, ce
       real(pReal), dimension(3,3) :: F
-    end function phase_mechanical_getF
+    end function phase_F
 
     module function mechanical_F_e(ph,me) result(F_e)
       integer, intent(in) :: ph,me
       real(pReal), dimension(3,3) :: F_e
     end function mechanical_F_e
 
-    module function phase_mechanical_getP(co,ce) result(P)
+    module function phase_P(co,ce) result(P)
       integer, intent(in) :: co, ce
       real(pReal), dimension(3,3) :: P
-    end function phase_mechanical_getP
+    end function phase_P
 
     module function phase_damage_get_phi(co,ip,el) result(phi)
       integer, intent(in) :: co, ip, el
@@ -181,10 +181,10 @@ module phase
     end function damage_phi
 
 
-    module subroutine phase_mechanical_setF(F,co,ce)
+    module subroutine phase_set_F(F,co,ce)
       real(pReal), dimension(3,3), intent(in) :: F
       integer, intent(in) :: co, ce
-    end subroutine phase_mechanical_setF
+    end subroutine phase_set_F
 
     module subroutine phase_thermal_setField(T,dot_T, co,ce)
       real(pReal), intent(in) :: T, dot_T
@@ -227,14 +227,13 @@ module phase
     end function phase_homogenizedC
 
 
-    module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
+    module function phase_damage_phi_dot(phi,co,ce) result(phi_dot)
-      integer, intent(in) :: ce
+      integer, intent(in) :: ce,co
       real(pReal), intent(in) :: &
         phi                                                                                         !< damage parameter
-      real(pReal), intent(inout) :: &
+      real(pReal) :: &
-        phiDot, &
+        phi_dot
-        dPhiDot_dPhi
+    end function phase_damage_phi_dot
-    end subroutine phase_damage_getRateAndItsTangents
 
     module subroutine phase_thermal_getRate(TDot, ph,me)
       integer, intent(in) :: ph, me
@@ -301,7 +300,7 @@ module phase
   public :: &
     phase_init, &
     phase_homogenizedC, &
-    phase_damage_getRateAndItsTangents, &
+    phase_damage_phi_dot, &
     phase_thermal_getRate, &
     phase_results, &
     phase_allocateState, &
@@ -321,9 +320,9 @@ module phase
     phase_thermal_setField, &
     phase_damage_set_phi, &
     phase_damage_get_phi, &
-    phase_mechanical_getP, &
+    phase_P, &
-    phase_mechanical_setF, &
+    phase_set_F, &
-    phase_mechanical_getF
+    phase_F
 
 contains
 
@@ -333,8 +332,7 @@ contains
 subroutine phase_init
 
   integer :: &
-    ph, &                                                                                            !< counter in phase loop
+    ph
-    so                                                                                               !< counter in source loop
   class (tNode), pointer :: &
     debug_constitutive, &
     materials, &
@@ -476,7 +474,6 @@ subroutine crystallite_init()
     co, &                                                                                           !< counter in integration point component loop
     ip, &                                                                                           !< counter in integration point loop
     el, &                                                                                           !< counter in element loop
-    so, &
     cMax, &                                                                                         !< maximum number of  integration point components
     iMax, &                                                                                         !< maximum number of integration points
     eMax                                                                                            !< maximum number of elements
@@ -591,7 +588,7 @@ function crystallite_push33ToRef(co,ce, tensor33)
 
   ph = material_phaseID(co,ce)
   en = material_phaseEntry(co,ce)
-  T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct?
+  T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct?
 
   crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
 

src/phase_damage.f90

@@ -65,43 +65,6 @@ submodule(phase) damage
       integer, intent(in) :: ph,me
     end subroutine isoductile_dotState
 
-
-    module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
-      integer, intent(in) :: ph,me
-      real(pReal),  intent(in) :: &
-        phi                                                                                           !< damage parameter
-      real(pReal),  intent(out) :: &
-        localphiDot, &
-        dLocalphiDot_dPhi
-    end subroutine anisobrittle_getRateAndItsTangent
-
-    module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
-      integer, intent(in) :: ph,me
-      real(pReal),  intent(in) :: &
-        phi                                                                                           !< damage parameter
-      real(pReal),  intent(out) :: &
-        localphiDot, &
-        dLocalphiDot_dPhi
-    end subroutine anisoductile_getRateAndItsTangent
-
-    module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
-      integer, intent(in) :: ph,me
-      real(pReal),  intent(in) :: &
-        phi                                                                                           !< damage parameter
-      real(pReal),  intent(out) :: &
-        localphiDot, &
-        dLocalphiDot_dPhi
-    end subroutine isobrittle_getRateAndItsTangent
-
-    module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
-      integer, intent(in) :: ph,me
-      real(pReal),  intent(in) :: &
-        phi                                                                                           !< damage parameter
-      real(pReal),  intent(out) :: &
-        localphiDot, &
-        dLocalphiDot_dPhi
-    end subroutine isoductile_getRateAndItsTangent
-
     module subroutine anisobrittle_results(phase,group)
       integer,          intent(in) :: phase
       character(len=*), intent(in) :: group
@@ -179,53 +142,30 @@ end subroutine damage_init
 !----------------------------------------------------------------------------------------------
 !< @brief returns local part of nonlocal damage driving force
 !----------------------------------------------------------------------------------------------
-module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
+module function phase_damage_phi_dot(phi,co,ce) result(phi_dot)
 
-  integer, intent(in) :: ce
+  integer, intent(in) :: ce,co
   real(pReal), intent(in) :: &
     phi                                                                                             !< damage parameter
-  real(pReal), intent(inout) :: &
-    phiDot, &
-    dPhiDot_dPhi
-
   real(pReal) :: &
-    localphiDot, &
+    phi_dot
-    dLocalphiDot_dPhi
+
   integer :: &
     ph, &
-    co, &
+    en
-    me
 
-   phiDot = 0.0_pReal
+  ph = material_phaseID(co,ce)
-   dPhiDot_dPhi = 0.0_pReal
+  en = material_phaseEntry(co,ce)
 
-   do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
+  select case(phase_source(ph))
-     ph = material_phaseID(co,ce)
+    case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID,DAMAGE_ANISODUCTILE_ID)
-     me = material_phaseEntry(co,ce)
+      phi_dot = 1.0_pReal &
+              - phi*damageState(ph)%state(1,en)
+    case default
+      phi_dot = 0.0_pReal
+  end select
 
-       select case(phase_source(ph))
+end function phase_damage_phi_dot
-         case (DAMAGE_ISOBRITTLE_ID)
-           call isobrittle_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
-
-         case (DAMAGE_ISODUCTILE_ID)
-           call isoductile_getRateAndItsTangent  (localphiDot, dLocalphiDot_dPhi, phi, ph, me)
-
-         case (DAMAGE_ANISOBRITTLE_ID)
-           call anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
-
-         case (DAMAGE_ANISODUCTILE_ID)
-           call anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
-
-         case default
-         localphiDot = 0.0_pReal
-         dLocalphiDot_dPhi = 0.0_pReal
-
-      end select
-      phiDot = phiDot + localphiDot
-      dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
-  enddo
-
-end subroutine phase_damage_getRateAndItsTangents
 
 
 

src/phase_damage_anisobrittle.f90

@@ -120,9 +120,6 @@ module subroutine anisobrittle_dotState(S, ph,me)
     S
 
   integer :: &
-    sourceOffset, &
-    damageOffset, &
-    homog, &
     i
   real(pReal) :: &
     traction_d, traction_t, traction_n, traction_crit
@@ -148,29 +145,6 @@ module subroutine anisobrittle_dotState(S, ph,me)
 end subroutine anisobrittle_dotState
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns local part of nonlocal damage driving force
-!--------------------------------------------------------------------------------------------------
-module subroutine anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
-
-  integer, intent(in) :: &
-    ph, &
-    me
-  real(pReal),  intent(in) :: &
-    phi
-  real(pReal),  intent(out) :: &
-    localphiDot, &
-    dLocalphiDot_dPhi
-
-
-  dLocalphiDot_dPhi = -damageState(ph)%state(1,me)
-
-  localphiDot = 1.0_pReal &
-              + dLocalphiDot_dPhi*phi
-
-end subroutine anisobrittle_getRateAndItsTangent
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------

src/phase_damage_anisoductile.f90

@@ -113,29 +113,6 @@ module subroutine anisoductile_dotState(ph,me)
 end subroutine anisoductile_dotState
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns local part of nonlocal damage driving force
-!--------------------------------------------------------------------------------------------------
-module subroutine anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph,me)
-
-  integer, intent(in) :: &
-    ph, &
-    me
-  real(pReal),  intent(in) :: &
-    phi
-  real(pReal),  intent(out) :: &
-    localphiDot, &
-    dLocalphiDot_dPhi
-
-
-  dLocalphiDot_dPhi = -damageState(ph)%state(1,me)
-
-  localphiDot = 1.0_pReal &
-              + dLocalphiDot_dPhi*phi
-
-end subroutine anisoductile_getRateAndItsTangent
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------

src/phase_damage_isobrittle.f90

@@ -113,29 +113,6 @@ module subroutine isobrittle_deltaState(C, Fe, ph,me)
 end subroutine isobrittle_deltaState
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns local part of nonlocal damage driving force
-!--------------------------------------------------------------------------------------------------
-module subroutine isobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi,  ph, me)
-
-  integer, intent(in) :: &
-    ph, me
-  real(pReal),  intent(in) :: &
-    phi
-  real(pReal),  intent(out) :: &
-    localphiDot, &
-    dLocalphiDot_dPhi
-
-
-  associate(prm => param(ph))
-    localphiDot = 1.0_pReal &
-                - phi*damageState(ph)%state(1,me)
-    dLocalphiDot_dPhi = - damageState(ph)%state(1,me)
-  end associate
-
-end subroutine isobrittle_getRateAndItsTangent
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------

src/phase_damage_isoductile.f90

@@ -103,29 +103,6 @@ module subroutine isoductile_dotState(ph, me)
 end subroutine isoductile_dotState
 
 
-!--------------------------------------------------------------------------------------------------
-!> @brief returns local part of nonlocal damage driving force
-!--------------------------------------------------------------------------------------------------
-module subroutine isoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ph, me)
-
-  integer, intent(in) :: &
-    ph, &
-    me
-  real(pReal),  intent(in) :: &
-    phi
-  real(pReal),  intent(out) :: &
-    localphiDot, &
-    dLocalphiDot_dPhi
-
-
-  dLocalphiDot_dPhi = -damageState(ph)%state(1,me)
-
-  localphiDot = 1.0_pReal &
-              + dLocalphiDot_dPhi*phi
-
-end subroutine isoductile_getRateAndItsTangent
-
-
 !--------------------------------------------------------------------------------------------------
 !> @brief writes results to HDF5 output file
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical.f90

@@ -1421,20 +1421,6 @@ module function mechanical_L_p(ph,me) result(L_p)
 end function mechanical_L_p
 
 
-!----------------------------------------------------------------------------------------------
-!< @brief Get deformation gradient (for use by homogenization)
-!----------------------------------------------------------------------------------------------
-module function phase_mechanical_getF(co,ce) result(F)
-
-  integer, intent(in) :: co, ce
-  real(pReal), dimension(3,3) :: F
-
-
-  F = phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce))
-
-end function phase_mechanical_getF
-
-
 !----------------------------------------------------------------------------------------------
 !< @brief Get elastic deformation gradient (for use by non-mech physics)
 !----------------------------------------------------------------------------------------------
@@ -1449,11 +1435,10 @@ module function mechanical_F_e(ph,me) result(F_e)
 end function mechanical_F_e
 
 
-
 !----------------------------------------------------------------------------------------------
 !< @brief Get second Piola-Kichhoff stress (for use by homogenization)
 !----------------------------------------------------------------------------------------------
-module function phase_mechanical_getP(co,ce) result(P)
+module function phase_P(co,ce) result(P)
 
   integer, intent(in) :: co, ce
   real(pReal), dimension(3,3) :: P
@@ -1461,11 +1446,27 @@ module function phase_mechanical_getP(co,ce) result(P)
 
   P = phase_mechanical_P(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce))
 
-end function phase_mechanical_getP
+end function phase_P
 
 
-! setter for homogenization
+!----------------------------------------------------------------------------------------------
-module subroutine phase_mechanical_setF(F,co,ce)
+!< @brief Get deformation gradient (for use by homogenization)
+!----------------------------------------------------------------------------------------------
+module function phase_F(co,ce) result(F)
+
+  integer, intent(in) :: co, ce
+  real(pReal), dimension(3,3) :: F
+
+
+  F = phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce))
+
+end function phase_F
+
+
+!----------------------------------------------------------------------------------------------
+!< @brief Set deformation gradient (for use by homogenization)
+!----------------------------------------------------------------------------------------------
+module subroutine phase_set_F(F,co,ce)
 
   real(pReal), dimension(3,3), intent(in) :: F
   integer, intent(in) :: co, ce
@@ -1473,7 +1474,7 @@ module subroutine phase_mechanical_setF(F,co,ce)
 
   phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) = F
 
-end subroutine phase_mechanical_setF
+end subroutine phase_set_F
 
 
 end submodule mechanical

src/phase_mechanical_eigen.f90

@@ -46,7 +46,6 @@ module subroutine eigendeformation_init(phases)
   class(tNode), pointer :: &
     phase, &
     kinematics, &
-    damage, &
     mechanics
 
   print'(/,a)', ' <<<+-  phase:mechanical:eigen init  -+>>>'