added brittle(elastic energy release rate)/ductile(vacancy condensation) phase field damage model coupled to vacancy concentration

simplified vacancy_generation

code/Makefile

@@ -326,7 +326,7 @@ COMPILE         =$(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATI
 COMPILE_MAXOPTI =$(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90))
 ###################################################################################################
 DAMAGE_FILES = \
-   damage_none.o damage_isoBrittle.o damage_isoDuctile.o damage_gurson.o damage_anisoBrittle.o damage_anisoDuctile.o
+   damage_none.o damage_isoBrittle.o damage_isoDuctile.o damage_gurson.o damage_anisoBrittle.o damage_anisoDuctile.o damage_phaseField.o
 
 THERMAL_FILES = \
    thermal_isothermal.o thermal_adiabatic.o
@@ -496,6 +496,9 @@ damage_anisoBrittle.o:         damage_anisoBrittle.f90 \
 damage_anisoDuctile.o:         damage_anisoDuctile.f90 \
                                lattice.o
                                
+damage_phaseField.o:           damage_phaseField.f90 \
+                               lattice.o
+                               
 damage_gurson.o:               damage_gurson.f90 \
                                lattice.o
 

code/constitutive.f90

@@ -47,6 +47,7 @@ module constitutive
    constitutive_putLocalVacancyConcentration, &
    constitutive_getVacancyConcentration, &
    constitutive_getVacancyDiffusion33, &
+   constitutive_getVacancyMobility33, &
    constitutive_postResults
  
  private :: &
@@ -131,6 +132,7 @@ subroutine constitutive_init(temperature_init)
    LOCAL_DAMAGE_anisoBrittle_ID, &
    LOCAL_DAMAGE_anisoDuctile_ID, &
    LOCAL_DAMAGE_gurson_ID, &
+   LOCAL_DAMAGE_phaseField_ID, &
    LOCAL_THERMAL_isothermal_ID, &
    LOCAL_THERMAL_adiabatic_ID, &
    LOCAL_VACANCY_constant_ID, &
@@ -141,6 +143,7 @@ subroutine constitutive_init(temperature_init)
    LOCAL_DAMAGE_anisoBrittle_LABEL, &
    LOCAL_DAMAGE_anisoDuctile_LABEL, &
    LOCAL_DAMAGE_gurson_LABEL, &
+   LOCAL_DAMAGE_phaseField_label, &
    LOCAL_THERMAL_isothermal_label, &
    LOCAL_THERMAL_adiabatic_label, &
    LOCAL_VACANCY_constant_label, &
@@ -165,6 +168,7 @@ subroutine constitutive_init(temperature_init)
  use damage_anisoDuctile
  use damage_anisoBrittle
  use damage_gurson
+ use damage_phaseField
  use thermal_isothermal
  use thermal_adiabatic
  use vacancy_constant
@@ -211,6 +215,7 @@ subroutine constitutive_init(temperature_init)
  if (any(phase_damage == LOCAL_DAMAGE_anisoBrittle_ID))     call damage_anisoBrittle_init(FILEUNIT)
  if (any(phase_damage == LOCAL_DAMAGE_anisoductile_ID))     call damage_anisoDuctile_init(FILEUNIT)
  if (any(phase_damage == LOCAL_DAMAGE_gurson_ID))           call damage_gurson_init(FILEUNIT)
+ if (any(phase_damage == LOCAL_DAMAGE_phaseField_ID))       call damage_phaseField_init(FILEUNIT)
  close(FILEUNIT)
  
 !--------------------------------------------------------------------------------------------------
@@ -324,6 +329,11 @@ subroutine constitutive_init(temperature_init)
        thisNoutput => damage_gurson_Noutput
        thisOutput => damage_gurson_output
        thisSize   => damage_gurson_sizePostResult
+     case (LOCAL_DAMAGE_phaseField_ID)
+       outputName = LOCAL_DAMAGE_phaseField_label
+       thisNoutput => damage_phaseField_Noutput
+       thisOutput => damage_phaseField_output
+       thisSize   => damage_phaseField_sizePostResult
      case default
        knownDamage = .false.
    end select   
@@ -509,9 +519,12 @@ function constitutive_damagedC(ipc,ip,el)
  use material, only: &
    material_phase, &
    LOCAL_DAMAGE_isoBrittle_ID, &
+   LOCAL_DAMAGE_phaseField_ID, &   
    phase_damage
  use damage_isoBrittle, only: &
    damage_isoBrittle_getDamagedC66  
+ use damage_phaseField, only: &
+   damage_phaseField_getDamagedC66  
 
  implicit none
  real(pReal), dimension(6,6) :: constitutive_damagedC
@@ -524,6 +537,9 @@ function constitutive_damagedC(ipc,ip,el)
    case (LOCAL_DAMAGE_isoBrittle_ID)
      constitutive_damagedC = damage_isoBrittle_getDamagedC66(constitutive_homogenizedC(ipc,ip,el), &
                                                              ipc,ip,el)
+   case (LOCAL_DAMAGE_phaseField_ID)
+     constitutive_damagedC = damage_phaseField_getDamagedC66(constitutive_homogenizedC(ipc,ip,el), &
+                                                             ipc,ip,el)
    case default
      constitutive_damagedC = constitutive_homogenizedC(ipc,ip,el)
 
@@ -547,7 +563,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el)
    PLASTICITY_nonlocal_ID, &
    LOCAL_DAMAGE_isoBrittle_ID, &
    LOCAL_DAMAGE_isoDuctile_ID, &
-   LOCAL_DAMAGE_gurson_ID
+   LOCAL_DAMAGE_gurson_ID, &
+   LOCAL_DAMAGE_phaseField_ID
 
  use constitutive_titanmod, only: &
    constitutive_titanmod_microstructure
@@ -559,6 +576,10 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el)
    constitutive_dislokmc_microstructure
  use damage_gurson, only: &
    damage_gurson_microstructure
+ use damage_gurson, only: &
+   damage_gurson_microstructure
+ use damage_phaseField, only: &
+   damage_phaseField_microstructure  
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -587,6 +608,10 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, ipc, ip, el)
  select case (phase_damage(material_phase(ipc,ip,el)))
    case (LOCAL_DAMAGE_gurson_ID)
      call damage_gurson_microstructure(ipc, ip, el)
+   case (LOCAL_DAMAGE_phaseField_ID)
+     call damage_phaseField_microstructure(constitutive_homogenizedC(ipc,ip,el), Fe, &
+                                           constitutive_getVacancyConcentration(ipc, ip, el), &
+                                           ipc, ip, el)
 
  end select
 
@@ -1037,6 +1062,7 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
    LOCAL_DAMAGE_anisoDuctile_ID, &
    LOCAL_DAMAGE_anisoBrittle_ID, &
    LOCAL_DAMAGE_gurson_ID, &
+   LOCAL_DAMAGE_phaseField_ID, &
    LOCAL_VACANCY_generation_ID
  use constitutive_j2, only:  &
    constitutive_j2_dotState
@@ -1060,6 +1086,8 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
    damage_anisoDuctile_dotState
  use damage_gurson, only: &
    damage_gurson_dotState
+ use damage_phaseField, only: &
+   damage_phaseField_dotState
  use vacancy_generation, only: &
    vacancy_generation_dotState
 
@@ -1121,6 +1149,8 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su
      call damage_anisoDuctile_dotState(nSlip, accumulatedSlip, ipc, ip, el)
    case (LOCAL_DAMAGE_gurson_ID)
      call damage_gurson_dotState(Tstar_v, Lp, ipc, ip, el)
+   case (LOCAL_DAMAGE_phaseField_ID)
+     call damage_phaseField_dotState(constitutive_getVacancyConcentration(ipc, ip, el), ipc, ip, el)
  end select
 
  select case (phase_vacancy(material_phase(ipc,ip,el)))
@@ -1216,6 +1246,7 @@ function constitutive_getLocalDamage(ipc, ip, el)
    LOCAL_DAMAGE_anisoBrittle_ID, &
    LOCAL_DAMAGE_anisoDuctile_ID, &
    LOCAL_DAMAGE_gurson_ID, &
+   LOCAL_DAMAGE_phaseField_ID, &
    phase_damage
  use damage_isoBrittle, only: &
    damage_isoBrittle_getLocalDamage
@@ -1227,6 +1258,8 @@ function constitutive_getLocalDamage(ipc, ip, el)
    damage_anisoDuctile_getLocalDamage
  use damage_gurson, only: &
    damage_gurson_getLocalDamage
+ use damage_phaseField, only: &
+   damage_phaseField_getLocalDamage
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -1249,12 +1282,14 @@ function constitutive_getLocalDamage(ipc, ip, el)
      constitutive_getLocalDamage = damage_anisoBrittle_getLocalDamage(ipc, ip, el)
 
    case (LOCAL_DAMAGE_anisoDuctile_ID)     
-     
      constitutive_getLocalDamage = damage_anisoDuctile_getLocalDamage(ipc, ip, el)
      
    case (LOCAL_DAMAGE_gurson_ID)
      constitutive_getLocalDamage = damage_gurson_getLocalDamage(ipc, ip, el)
 
+   case (LOCAL_DAMAGE_phaseField_ID)
+     constitutive_getLocalDamage = damage_phaseField_getLocalDamage(ipc, ip, el)
+
  end select
 
 end function constitutive_getLocalDamage
@@ -1272,6 +1307,7 @@ subroutine constitutive_putLocalDamage(ipc, ip, el, localDamage)
    LOCAL_DAMAGE_anisoBrittle_ID, &
    LOCAL_DAMAGE_anisoDuctile_ID, &
    LOCAL_DAMAGE_gurson_ID, &
+   LOCAL_DAMAGE_phaseField_ID, &
    phase_damage
  use damage_isoBrittle, only: &
    damage_isoBrittle_putLocalDamage
@@ -1283,6 +1319,8 @@ subroutine constitutive_putLocalDamage(ipc, ip, el, localDamage)
    damage_anisoDuctile_putLocalDamage
  use damage_gurson, only: &
    damage_gurson_putLocalDamage
+ use damage_phaseField, only: &
+   damage_phaseField_putLocalDamage
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -1308,6 +1346,9 @@ subroutine constitutive_putLocalDamage(ipc, ip, el, localDamage)
    case (LOCAL_DAMAGE_gurson_ID)
      call damage_gurson_putLocalDamage(ipc, ip, el, localDamage)
  
+   case (LOCAL_DAMAGE_phaseField_ID)
+     call damage_phaseField_putLocalDamage(ipc, ip, el, localDamage)
+
  end select
 
 end subroutine constitutive_putLocalDamage
@@ -1326,6 +1367,7 @@ function constitutive_getDamage(ipc, ip, el)
    LOCAL_DAMAGE_anisoBrittle_ID, &
    LOCAL_DAMAGE_anisoDuctile_ID, &
    LOCAL_DAMAGE_gurson_ID, &
+   LOCAL_DAMAGE_phaseField_ID, &
    phase_damage
  use damage_isoBrittle, only: &
    damage_isoBrittle_getDamage
@@ -1337,6 +1379,8 @@ function constitutive_getDamage(ipc, ip, el)
    damage_anisoDuctile_getDamage
  use damage_gurson, only: &
    damage_gurson_getDamage
+ use damage_phaseField, only: &
+   damage_phaseField_getDamage
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -1364,6 +1408,9 @@ function constitutive_getDamage(ipc, ip, el)
    case (LOCAL_DAMAGE_gurson_ID)
      constitutive_getDamage = damage_gurson_getDamage(ipc, ip, el)
 
+   case (LOCAL_DAMAGE_phaseField_ID)
+     constitutive_getDamage = damage_phaseField_getDamage(ipc, ip, el)
+
  end select
 
 end function constitutive_getDamage
@@ -1424,14 +1471,13 @@ function constitutive_getDamageDiffusion33(ipc, ip, el)
    lattice_DamageDiffusion33
  use material, only: &
    material_phase, &
-   LOCAL_DAMAGE_none_ID, &
+   phase_damage, &
    LOCAL_DAMAGE_isoBrittle_ID, &
-   LOCAL_DAMAGE_isoDuctile_ID, &
+   LOCAL_DAMAGE_phaseField_ID
-   LOCAL_DAMAGE_anisoBrittle_ID, &
-   LOCAL_DAMAGE_gurson_ID, &
-   phase_damage
  use damage_isoBrittle, only: &
    damage_isoBrittle_getDamageDiffusion33
+ use damage_phaseField, only: &
+   damage_phaseField_getDamageDiffusion33
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -1445,6 +1491,8 @@ function constitutive_getDamageDiffusion33(ipc, ip, el)
  select case(phase_damage(material_phase(ipc,ip,el)))                                                   
    case (LOCAL_DAMAGE_isoBrittle_ID)
     constitutive_getDamageDiffusion33 = damage_isoBrittle_getDamageDiffusion33(ipc, ip, el)
+   case (LOCAL_DAMAGE_phaseField_ID)
+    constitutive_getDamageDiffusion33 = damage_phaseField_getDamageDiffusion33(ipc, ip, el)
     
  end select
 
@@ -1559,7 +1607,7 @@ function constitutive_getLocalVacancyConcentration(ipc, ip, el)
    LOCAL_VACANCY_generation_ID, &
    phase_vacancy
  use vacancy_generation, only: &
-   vacancy_generation_getConcentration
+   vacancy_generation_getLocalConcentration
  use lattice, only: &
    lattice_equilibriumVacancyConcentration
 
@@ -1576,7 +1624,7 @@ function constitutive_getLocalVacancyConcentration(ipc, ip, el)
        lattice_equilibriumVacancyConcentration(material_phase(ipc,ip,el))
      
    case (LOCAL_VACANCY_generation_ID)
-     constitutive_getLocalVacancyConcentration = vacancy_generation_getConcentration(ipc, ip, el)
+     constitutive_getLocalVacancyConcentration = vacancy_generation_getLocalConcentration(ipc, ip, el)
  end select
 
 end function constitutive_getLocalVacancyConcentration
@@ -1592,7 +1640,7 @@ subroutine constitutive_putLocalVacancyConcentration(ipc, ip, el, localVacancyCo
    LOCAL_VACANCY_generation_ID, &
    phase_vacancy
  use vacancy_generation, only: &
-   vacancy_generation_putConcentration
+   vacancy_generation_putLocalConcentration
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -1604,7 +1652,7 @@ subroutine constitutive_putLocalVacancyConcentration(ipc, ip, el, localVacancyCo
  
  select case (phase_vacancy(material_phase(ipc,ip,el)))
    case (LOCAL_VACANCY_generation_ID)
-     call vacancy_generation_putConcentration(ipc, ip, el, localVacancyConcentration)
+     call vacancy_generation_putLocalConcentration(ipc, ip, el, localVacancyConcentration)
 
  end select
 
@@ -1617,13 +1665,12 @@ function constitutive_getVacancyConcentration(ipc, ip, el)
  use prec, only: &
    pReal
  use material, only: &
-   mappingHomogenization, &
    material_phase, &
-   fieldVacancy, &
+   LOCAL_VACANCY_constant_ID, &
-   field_vacancy_type, &
+   LOCAL_VACANCY_generation_ID, &
-   FIELD_VACANCY_local_ID, &
+   phase_vacancy
-   FIELD_VACANCY_nonlocal_ID, &
+ use vacancy_generation, only: &
-   material_homog
+   vacancy_generation_getConcentration
  use lattice, only: &
    lattice_equilibriumVacancyConcentration
  implicit none
@@ -1634,13 +1681,13 @@ function constitutive_getVacancyConcentration(ipc, ip, el)
    el                                                                                               !< element number
  real(pReal) :: constitutive_getVacancyConcentration
  
-   select case(field_vacancy_type(material_homog(ip,el)))                                                   
+ select case (phase_vacancy(material_phase(ipc,ip,el)))
-     case (FIELD_VACANCY_local_ID)
+   case (LOCAL_VACANCY_constant_ID)
-      constitutive_getVacancyConcentration = constitutive_getLocalVacancyConcentration(ipc, ip, el)      
+     constitutive_getVacancyConcentration = &
+       lattice_equilibriumVacancyConcentration(material_phase(ipc,ip,el))
      
-     case (FIELD_VACANCY_nonlocal_ID)
+   case (LOCAL_VACANCY_generation_ID)
-      constitutive_getVacancyConcentration = fieldVacancy(material_homog(ip,el))% &
+     constitutive_getVacancyConcentration = vacancy_generation_getConcentration(ipc, ip, el)
-        field(1,mappingHomogenization(1,ip,el))                                                     ! Taylor type 
  end select
 
 end function constitutive_getVacancyConcentration
@@ -1667,22 +1714,53 @@ function constitutive_getVacancyDiffusion33(ipc, ip, el)
    el                                                                                               !< element number
  real(pReal), dimension(3,3) :: &
    constitutive_getVacancyDiffusion33
+ 
+ select case(phase_vacancy(material_phase(ipc,ip,el)))                                                   
+   case (LOCAL_VACANCY_generation_ID)
+    constitutive_getVacancyDiffusion33 = &
+      vacancy_generation_getVacancyDiffusion33(ipc,ip,el)
+    
+ end select
+
+end function constitutive_getVacancyDiffusion33
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns vacancy diffusion tensor
+!--------------------------------------------------------------------------------------------------
+function constitutive_getVacancyMobility33(ipc, ip, el)
+ use prec, only: &
+   pReal
+ use lattice, only: &
+   lattice_VacancyDiffusion33
+ use material, only: &
+   material_phase, &
+   LOCAL_VACANCY_generation_ID, &
+   phase_vacancy
+ use vacancy_generation, only: &
+   vacancy_generation_getVacancyMobility33
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), dimension(3,3) :: &
+   constitutive_getVacancyMobility33
  real(pReal), dimension(:), allocatable :: &
    accumulatedSlip
  integer(pInt) :: &
    nSlip
  
- constitutive_getVacancyDiffusion33 = lattice_VacancyDiffusion33(1:3,1:3,material_phase(ipc,ip,el))
  select case(phase_vacancy(material_phase(ipc,ip,el)))                                                   
    case (LOCAL_VACANCY_generation_ID)
     call constitutive_getAccumulatedSlip(nSlip,accumulatedSlip,ipc,ip,el)
-    constitutive_getVacancyDiffusion33 = &
+    constitutive_getVacancyMobility33 = &
-      vacancy_generation_getVacancyDiffusion33(nSlip,accumulatedSlip,constitutive_getTemperature(ipc,ip,el), &
+      vacancy_generation_getVacancyMobility33(nSlip,accumulatedSlip,constitutive_getTemperature(ipc,ip,el), &
                                                ipc,ip,el)
     
  end select
 
-end function constitutive_getVacancyDiffusion33
+end function constitutive_getVacancyMobility33
 
 !--------------------------------------------------------------------------------------------------
 !> @brief returns accumulated slip on each system defined

code/damage_phaseField.f90

@@ -0,0 +1,520 @@
+!--------------------------------------------------------------------------------------------------
+! $Id: damage_phaseField.f90 3736 2014-11-21 13:12:54Z MPIE\l.sharma $
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine incoprorating isotropic brittle damage
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module damage_phaseField
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,          public, protected :: &
+   damage_phaseField_sizePostResults                                                           !< cumulative size of post results
+
+ integer(pInt),                       dimension(:,:),         allocatable, target,  public :: &
+   damage_phaseField_sizePostResult                                                            !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target,  public :: &
+   damage_phaseField_output                                                                    !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target,  public :: &
+   damage_phaseField_Noutput                                                                   !< number of outputs per instance of this damage 
+
+ real(pReal),                         dimension(:),           allocatable,         private :: &
+   damage_phaseField_aTol, &
+   damage_phaseField_surfaceEnergy, &
+   damage_phaseField_vacancyFormationEnergy
+
+ enum, bind(c) 
+   enumerator :: undefined_ID, &
+                 local_damage_ID
+ end enum                                                 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 ToDo
+ 
+ integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
+   damage_phaseField_outputID                                                                  !< ID of each post result output
+
+
+ public :: &
+   damage_phaseField_init, &
+   damage_phaseField_stateInit, &
+   damage_phaseField_aTolState, &
+   damage_phaseField_dotState, &
+   damage_phaseField_microstructure, &
+   damage_phaseField_getDamage, &
+   damage_phaseField_putLocalDamage, &
+   damage_phaseField_getLocalDamage, &
+   damage_phaseField_getDamageDiffusion33, &
+   damage_phaseField_getDamagedC66, &
+   damage_phaseField_postResults
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use mesh, only: &
+   mesh_maxNips, &
+   mesh_NcpElems
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   homogenization_maxNgrains, &
+   phase_damage, &
+   phase_damageInstance, &
+   phase_Noutput, &
+   LOCAL_damage_phaseField_label, &
+   LOCAL_damage_phaseField_ID, &
+   material_phase, &  
+   damageState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   worldrank, &
+   numerics_integrator
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
+ integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
+ integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,o
+ integer(pInt) :: sizeState, sizeDotState
+ integer(pInt) :: NofMyPhase   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  damage_'//LOCAL_damage_phaseField_label//' init  -+>>>'
+   write(6,'(a)')     ' $Id: damage_phaseField.f90 3736 2014-11-21 13:12:54Z MPIE\l.sharma $'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_damage == LOCAL_damage_phaseField_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(damage_phaseField_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(damage_phaseField_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(damage_phaseField_output(maxval(phase_Noutput),maxNinstance))
+          damage_phaseField_output = ''
+ allocate(damage_phaseField_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
+ allocate(damage_phaseField_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(damage_phaseField_surfaceEnergy(maxNinstance),                       source=0.0_pReal) 
+ allocate(damage_phaseField_vacancyFormationEnergy(maxNinstance),              source=0.0_pReal) 
+ allocate(damage_phaseField_aTol(maxNinstance),                                source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (phase_damage(phase) == LOCAL_damage_phaseField_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = phase_damageInstance(phase)                                                     ! which instance of my damage is present phase
+     positions = IO_stringPos(line,MAXNCHUNKS)
+     tag = IO_lc(IO_stringValue(line,positions,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('(output)')
+         select case(IO_lc(IO_stringValue(line,positions,2_pInt)))
+           case ('local_damage')
+             damage_phaseField_Noutput(instance) = damage_phaseField_Noutput(instance) + 1_pInt
+             damage_phaseField_outputID(damage_phaseField_Noutput(instance),instance) = local_damage_ID
+             damage_phaseField_output(damage_phaseField_Noutput(instance),instance) = &
+                                                       IO_lc(IO_stringValue(line,positions,2_pInt))
+          end select
+
+       case ('surfaceenergy')
+         damage_phaseField_surfaceEnergy(instance) = IO_floatValue(line,positions,2_pInt)
+
+       case ('vacancyformationenergy')
+         damage_phaseField_vacancyFormationEnergy(instance) = IO_floatValue(line,positions,2_pInt)
+
+       case ('atol_damage')
+         damage_phaseField_aTol(instance) = IO_floatValue(line,positions,2_pInt)
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+
+ sanityChecks: do phase = 1_pInt, size(phase_damage)   
+   myPhase: if (phase_damage(phase) == LOCAL_damage_phaseField_ID) then
+     NofMyPhase=count(material_phase==phase)
+     instance = phase_damageInstance(phase)
+!  sanity checks
+     if (damage_phaseField_aTol(instance) < 0.0_pReal) &
+       damage_phaseField_aTol(instance) = 1.0e-3_pReal                                              ! default absolute tolerance 1e-3
+     if (damage_phaseField_surfaceEnergy(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='surfaceEnergy ('//LOCAL_damage_phaseField_LABEL//')')
+   endif myPhase
+ enddo sanityChecks
+
+ initializeInstances: do phase = 1_pInt, size(phase_damage)
+   if (phase_damage(phase) == LOCAL_damage_phaseField_ID) then
+     NofMyPhase=count(material_phase==phase)
+     instance = phase_damageInstance(phase)
+!--------------------------------------------------------------------------------------------------
+!  Determine size of postResults array
+     outputsLoop: do o = 1_pInt,damage_phaseField_Noutput(instance)
+       select case(damage_phaseField_outputID(o,instance))
+         case(local_damage_ID)
+           mySize = 1_pInt
+       end select
+ 
+       if (mySize > 0_pInt) then  ! any meaningful output found
+          damage_phaseField_sizePostResult(o,instance) = mySize
+          damage_phaseField_sizePostResults(instance)  = damage_phaseField_sizePostResults(instance) + mySize
+       endif
+     enddo outputsLoop
+! Determine size of state array
+     sizeDotState              =   1_pInt
+     sizeState                 =   2_pInt
+                
+     damageState(phase)%sizeState = sizeState
+     damageState(phase)%sizeDotState = sizeDotState
+     damageState(phase)%sizePostResults = damage_phaseField_sizePostResults(instance)
+     allocate(damageState(phase)%aTolState           (sizeState),                source=0.0_pReal)
+     allocate(damageState(phase)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(damageState(phase)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(damageState(phase)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(damageState(phase)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(damageState(phase)%state_backup        (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(damageState(phase)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(damageState(phase)%deltaState          (sizeDotState,NofMyPhase),     source=0.0_pReal)
+     allocate(damageState(phase)%dotState_backup     (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(damageState(phase)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(damageState(phase)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(damageState(phase)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(damageState(phase)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)
+
+     call damage_phaseField_stateInit(phase)
+     call damage_phaseField_aTolState(phase,instance)
+   endif
+ 
+ enddo initializeInstances
+end subroutine damage_phaseField_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief sets the relevant  NEW state values for a given instance of this damage
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_stateInit(phase)
+ use material, only: &
+   damageState
+ 
+ implicit none
+ integer(pInt),              intent(in) :: phase                                                    !< number specifying the phase of the damage
+
+ real(pReal), dimension(damageState(phase)%sizeState) :: tempState
+
+ tempState = 1.0_pReal
+ damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:)))
+ damageState(phase)%state0 = damageState(phase)%state
+ damageState(phase)%partionedState0 = damageState(phase)%state
+end subroutine damage_phaseField_stateInit
+
+!--------------------------------------------------------------------------------------------------
+!> @brief sets the relevant state values for a given instance of this damage
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_aTolState(phase,instance)
+ use material, only: &
+  damageState
+
+ implicit none
+ integer(pInt), intent(in) ::  &
+   phase, &
+   instance                                                                                         ! number specifying the current instance of the damage
+ real(pReal), dimension(damageState(phase)%sizeState) :: tempTol
+
+ tempTol = damage_phaseField_aTol(instance)
+ damageState(phase)%aTolState = tempTol
+end subroutine damage_phaseField_aTolState
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_dotState(Cv, ipc, ip, el)
+ use material, only: &
+   mappingConstitutive, &
+   damageState
+ use lattice, only: &
+   lattice_DamageMobility
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal),  intent(in) :: &
+   Cv
+ integer(pInt) :: &
+   phase, constituent
+
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+
+ damageState(phase)%dotState(1,constituent) = &
+   (damageState(phase)%state(2,constituent)*(1.0_pReal - Cv) - &
+    damageState(phase)%state(1,constituent))/ &
+   lattice_DamageMobility(phase)
+
+end subroutine damage_phaseField_dotState
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_microstructure(C, Fe, Cv, ipc, ip, el)
+ use material, only: &
+   mappingConstitutive, &
+   phase_damageInstance, &
+   damageState
+ use math, only : &
+   math_mul33x33, &
+   math_mul66x6, &
+   math_Mandel33to6, &
+   math_transpose33, &
+   math_I3
+ use lattice, only: &
+   lattice_DamageMobility
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal),  intent(in), dimension(3,3) :: &
+   Fe
+ real(pReal),  intent(in), dimension(6,6) :: &
+   C
+ real(pReal),  intent(in) :: &
+   Cv
+ integer(pInt) :: &
+   phase, constituent, instance
+ real(pReal) :: &
+   strain(6), &
+   stress(6)
+
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ instance = phase_damageInstance(phase)
+
+ strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
+ stress = math_mul66x6(C,strain) 
+ 
+ damageState(phase)%state(2,constituent) = &
+   min(damageState(phase)%state0(2,constituent), &
+       damage_phaseField_surfaceEnergy(instance)/ &
+       (2.0_pReal*(sum(abs(stress*strain)) + Cv*damage_phaseField_vacancyFormationEnergy(instance))))
+ 
+end subroutine damage_phaseField_microstructure
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns damage
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_getDamage(ipc, ip, el)
+ use material, only: &
+   material_homog, &
+   mappingHomogenization, &
+   fieldDamage, &
+   field_damage_type, &
+   FIELD_DAMAGE_LOCAL_ID, &
+   FIELD_DAMAGE_NONLOCAL_ID
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal) :: damage_phaseField_getDamage
+ 
+ select case(field_damage_type(material_homog(ip,el)))                                                   
+   case (FIELD_DAMAGE_LOCAL_ID)
+    damage_phaseField_getDamage = damage_phaseField_getLocalDamage(ipc, ip, el)
+    
+   case (FIELD_DAMAGE_NONLOCAL_ID)
+    damage_phaseField_getDamage = fieldDamage(material_homog(ip,el))% &
+      field(1,mappingHomogenization(1,ip,el))                                                     ! Taylor type 
+
+ end select
+
+end function damage_phaseField_getDamage
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns temperature based on local damage model state layout 
+!--------------------------------------------------------------------------------------------------
+subroutine damage_phaseField_putLocalDamage(ipc, ip, el, localDamage)
+ use material, only: &
+   mappingConstitutive, &
+   damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),   intent(in) :: localDamage
+ 
+ damageState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)) = &
+   localDamage
+ 
+end subroutine damage_phaseField_putLocalDamage
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local damage
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_getLocalDamage(ipc, ip, el)
+ use material, only: &
+   mappingConstitutive, &
+   damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal) :: damage_phaseField_getLocalDamage
+ 
+ damage_phaseField_getLocalDamage = &
+   damageState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el))
+ 
+end function damage_phaseField_getLocalDamage
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns brittle damage diffusion tensor 
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_getDamageDiffusion33(ipc, ip, el)
+ use lattice, only: &
+   lattice_DamageDiffusion33
+ use material, only: &
+   mappingConstitutive, &
+   damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), dimension(3,3) :: &
+   damage_phaseField_getDamageDiffusion33
+ integer(pInt) :: &
+   phase, constituent
+ 
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ damage_phaseField_getDamageDiffusion33 = &
+   damageState(phase)%state(2,constituent)* &
+   lattice_DamageDiffusion33(1:3,1:3,phase)
+    
+end function damage_phaseField_getDamageDiffusion33
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns brittle damaged stiffness tensor 
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_getDamagedC66(C, ipc, ip, el)
+ use material, only: &
+   mappingConstitutive, &
+   damageState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),  intent(in), dimension(6,6) :: &
+   C
+ real(pReal),              dimension(6,6) :: &
+   damage_phaseField_getDamagedC66
+ integer(pInt) :: &
+   phase, constituent
+ real(pReal) :: &
+   damage  
+ 
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ damage = damage_phaseField_getDamage(ipc, ip, el)
+ 
+ damage_phaseField_getDamagedC66 = &
+   damage*damage*C
+    
+end function damage_phaseField_getDamagedC66
+
+!--------------------------------------------------------------------------------------------------
+!> @brief return array of constitutive results
+!--------------------------------------------------------------------------------------------------
+function damage_phaseField_postResults(ipc,ip,el)
+ use material, only: &
+   mappingConstitutive, &
+   phase_damageInstance,& 
+   damageState
+
+ implicit none
+ integer(pInt),              intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal), dimension(damage_phaseField_sizePostResults(phase_damageInstance(mappingConstitutive(2,ipc,ip,el)))) :: &
+   damage_phaseField_postResults
+
+ integer(pInt) :: &
+   instance, phase, constituent, o, c
+   
+ phase = mappingConstitutive(2,ipc,ip,el)
+ constituent = mappingConstitutive(1,ipc,ip,el)
+ instance = phase_damageInstance(phase)
+
+ c = 0_pInt
+ damage_phaseField_postResults = 0.0_pReal
+
+ do o = 1_pInt,damage_phaseField_Noutput(instance)
+    select case(damage_phaseField_outputID(o,instance))
+      case (local_damage_ID)
+        damage_phaseField_postResults(c+1_pInt) = damageState(phase)%state(2,constituent)
+        c = c + 1
+
+    end select
+ enddo
+end function damage_phaseField_postResults
+
+end module damage_phaseField

code/material.f90

@@ -33,6 +33,7 @@ module material
    LOCAL_DAMAGE_anisoBrittle_LABEL= 'anisobrittle', &
    LOCAL_DAMAGE_anisoDuctile_LABEL= 'anisoductile', &
    LOCAL_DAMAGE_gurson_LABEL      = 'gurson', &   
+   LOCAL_DAMAGE_phaseField_LABEL  = 'phasefield', &   
    LOCAL_THERMAL_isothermal_label = 'isothermal', &
    LOCAL_THERMAL_adiabatic_label  = 'adiabatic', &
    LOCAL_VACANCY_constant_label   = 'constant', &
@@ -69,7 +70,8 @@ module material
                  LOCAL_DAMAGE_isoDuctile_ID, &
                  LOCAL_DAMAGE_anisoBrittle_ID, &
                  LOCAL_DAMAGE_anisoDuctile_ID, &
-                 LOCAL_DAMAGE_gurson_ID
+                 LOCAL_DAMAGE_gurson_ID, &
+                 LOCAL_DAMAGE_phaseField_ID
  end enum
  enum, bind(c)
    enumerator :: LOCAL_THERMAL_isothermal_ID, &
@@ -244,6 +246,7 @@ module material
    LOCAL_DAMAGE_anisoBrittle_ID, &
    LOCAL_DAMAGE_anisoDuctile_ID, &
    LOCAL_DAMAGE_gurson_ID, &
+   LOCAL_DAMAGE_phaseField_ID, &
    LOCAL_THERMAL_isothermal_ID, &
    LOCAL_THERMAL_adiabatic_ID, &
    LOCAL_VACANCY_constant_ID, &
@@ -859,6 +862,8 @@ subroutine material_parsePhase(fileUnit,myPart)
              phase_damage(section) = LOCAL_DAMAGE_anisoDuctile_ID
            case (LOCAL_DAMAGE_gurson_label)
              phase_damage(section) = LOCAL_DAMAGE_gurson_ID
+           case (LOCAL_DAMAGE_phaseField_label)
+             phase_damage(section) = LOCAL_DAMAGE_phaseField_ID
            case default
              call IO_error(200_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt)))
          end select

code/vacancy_generation.f90

@@ -28,21 +28,11 @@ module vacancy_generation
    vacancy_generation_aTol, &
    vacancy_generation_freq, &
    vacancy_generation_formationEnergy, &
-   vacancy_generation_diffusionEnergy, &
+   vacancy_generation_migrationEnergy, &
-   vacancy_generation_diffusionCoeff0, &                                                      !< the temperature-independent pre-exponential of diffusion coefficient D_0
+   vacancy_generation_diffusionCoeff0, &                                                        !< the temperature-independent diffusion coefficient D_0
-   vacancy_generation_stressCoeff, &
+   vacancy_generation_atomicVol, &
-   vacancy_generation_jogHeight, &                                                              !< the height of jogs in Burgers vectors
+   vacancy_generation_surfaceEnergy, &
-   vacancy_generation_jogSeparation, &                                                          !< the jog seperation
+   vacancy_generation_plasticityCoeff
-   vacancy_generation_nLatticeSites, &                                                          !< the number of lattice sites per unit volume
-   vacancy_generation_burgersVec, &                                                             !< the Burgers vector
-   vacancy_generation_dislocationCoeff, &
-   vacancy_generation_equilibConcentration                                                      !< the equilibrium concentration of vacancy
-
- real(pReal),                         dimension(:),           allocatable,         public :: &
-   pore_nucleation_surfaceEnergy, &                                                             !< surface energy of metal which controls the necleation of pores
-   pore_nucleation_atomVolume, &                                                                !< the volume of atom
-   pore_nucleation_shellThickness, &                                                            !< the thickness of spherical shell surrounding the pore
-   pore_nucleation_concentrationCoeff0                                                          !< the pre-exponential of equilibrium concentration of critical pore
 
  real(pReal),                                                 parameter,           private :: &
    kB = 1.38e-23_pReal                                                                          !< Boltzmann constant in J/Kelvin
@@ -60,9 +50,11 @@ module vacancy_generation
    vacancy_generation_stateInit, &
    vacancy_generation_aTolState, &
    vacancy_generation_dotState, &
+   vacancy_generation_getLocalConcentration, &
+   vacancy_generation_putLocalConcentration, &
    vacancy_generation_getConcentration, &
-   vacancy_generation_putConcentration, &
    vacancy_generation_getVacancyDiffusion33, &
+   vacancy_generation_getVacancyMobility33, &
    vacancy_generation_postResults
 
 contains
@@ -141,21 +133,10 @@ subroutine vacancy_generation_init(fileUnit)
  allocate(vacancy_generation_aTol(maxNinstance),                                source=0.0_pReal) 
  allocate(vacancy_generation_freq(maxNinstance),                                source=0.0_pReal) 
  allocate(vacancy_generation_formationEnergy(maxNinstance),                     source=0.0_pReal) 
- allocate(vacancy_generation_diffusionEnergy(maxNinstance),                     source=0.0_pReal) 
+ allocate(vacancy_generation_migrationEnergy(maxNinstance),                     source=0.0_pReal) 
- allocate(vacancy_generation_stressCoeff(maxNinstance),                         source=0.0_pReal) 
+ allocate(vacancy_generation_atomicVol(maxNinstance),                           source=0.0_pReal) 
- allocate(vacancy_generation_jogHeight(maxNinstance),                           source=0.0_pReal)
+ allocate(vacancy_generation_surfaceEnergy(maxNinstance),                       source=0.0_pReal)
- allocate(vacancy_generation_jogSeparation(maxNinstance),                       source=0.0_pReal) 
+ allocate(vacancy_generation_plasticityCoeff(maxNinstance),                     source=0.0_pReal)
- allocate(vacancy_generation_nLatticeSites(maxNinstance),                       source=0.0_pReal)
- allocate(vacancy_generation_burgersVec(maxNinstance),                          source=0.0_pReal)
- allocate(vacancy_generation_diffusionCoeff0(maxNinstance),                     source=0.0_pReal)
- allocate(vacancy_generation_equilibConcentration(maxNinstance),                source=0.0_pReal)
- 
- allocate(vacancy_generation_dislocationCoeff(maxNinstance),                    source=0.0_pReal)
- 
- allocate(pore_nucleation_surfaceEnergy(maxNinstance),                          source=0.0_pReal)
- allocate(pore_nucleation_atomVolume(maxNinstance),                             source=0.0_pReal)
- allocate(pore_nucleation_shellThickness(maxNinstance),                         source=0.0_pReal)
- allocate(pore_nucleation_concentrationCoeff0(maxNinstance),                    source=0.0_pReal)
 
  rewind(fileUnit)
  phase = 0_pInt
@@ -190,50 +171,29 @@ subroutine vacancy_generation_init(fileUnit)
                                                        IO_lc(IO_stringValue(line,positions,2_pInt))
           end select
 
-       case ('atol_vacancygeneration')
+       case ('atolvacancygeneration')
          vacancy_generation_aTol(instance) = IO_floatValue(line,positions,2_pInt)
 
-       case ('vacancy_frequency')
+       case ('debyefrequency')
          vacancy_generation_freq(instance) = IO_floatValue(line,positions,2_pInt)
 
-       case ('vacancy_formationenergy')
+       case ('vacancyformationenergy')
          vacancy_generation_formationEnergy(instance) = IO_floatValue(line,positions,2_pInt)
 
-       case ('vacancy_equilibconcentration')
+       case ('vacancymigrationenergy')
-         vacancy_generation_equilibConcentration(instance) = IO_floatValue(line,positions,2_pInt)
+         vacancy_generation_migrationEnergy(instance) = IO_floatValue(line,positions,2_pInt)
 
-       case ('vacancy_diffusionenergy')
+       case ('vacancydiffusioncoeff0')
-         vacancy_generation_diffusionEnergy(instance) = IO_floatValue(line,positions,2_pInt)
-
-       case ('vacancy_diffusioncoeff0')
          vacancy_generation_diffusionCoeff0(instance) = IO_floatValue(line,positions,2_pInt)
 
-       case ('vacancy_stresscoeff')
+       case ('atomicvolume')
-         vacancy_generation_stressCoeff(instance) = IO_floatValue(line,positions,2_pInt)
+         vacancy_generation_atomicVol(instance) = IO_floatValue(line,positions,2_pInt)
          
-       case ('vacancy_jogheight')
+       case ('surfaceenergy')
-         vacancy_generation_jogHeight(instance) = IO_floatValue(line,positions,2_pInt)
+         vacancy_generation_surfaceEnergy(instance) = IO_floatValue(line,positions,2_pInt)
 
-       case ('vacancy_jogseparation')
+       case ('vacancyplasticitycoeff')
-         vacancy_generation_jogSeparation(instance) = IO_floatValue(line,positions,2_pInt)
+         vacancy_generation_plasticityCoeff(instance) = IO_floatValue(line,positions,2_pInt)
-
-       case ('vacancy_nlatticesites')
-         vacancy_generation_nLatticeSites(instance) = IO_floatValue(line,positions,2_pInt)
-
-       case ('vacancy_burgersvec')
-         vacancy_generation_burgersVec(instance) = IO_floatValue(line,positions,2_pInt)
-
-       case ('pore_surfacefnergy')
-         pore_nucleation_surfaceEnergy(instance) = IO_floatValue(line,positions,2_pInt)
-
-       case ('pore_atomvolume')
-         pore_nucleation_atomVolume(instance) = IO_floatValue(line,positions,2_pInt)
-
-       case ('pore_shellthickness')
-         pore_nucleation_shellThickness(instance) = IO_floatValue(line,positions,2_pInt)
-
-       case ('pore_concentrationcoeff0')
-         pore_nucleation_concentrationCoeff0(instance) = IO_floatValue(line,positions,2_pInt)
 
      end select
    endif; endif
@@ -244,14 +204,6 @@ subroutine vacancy_generation_init(fileUnit)
      NofMyPhase=count(material_phase==phase)
      instance = phase_vacancyInstance(phase)
 
-!--------------------------------------------------------------------------------------------------
-!  Calculate the coefficient for dislocation motion induced vacancy generation
-     vacancy_generation_dislocationCoeff(instance) = vacancy_generation_jogHeight(instance)/     &
-                                                     vacancy_generation_jogSeparation(instance)/ &
-                                                     vacancy_generation_nLatticeSites(instance)/ &
-                                                     vacancy_generation_burgersVec(instance)/    &
-                                                     vacancy_generation_burgersVec(instance)
-
 !--------------------------------------------------------------------------------------------------
 !  Determine size of postResults array
      outputsLoop: do o = 1_pInt,vacancy_generation_Noutput(instance)
@@ -337,17 +289,13 @@ end subroutine vacancy_generation_aTolState
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
 subroutine vacancy_generation_dotState(nSlip, accumulatedSlip, Tstar_v, Temperature, ipc, ip, el)
- use lattice, only: &
-   lattice_massDensity, &
-   lattice_specificHeat
  use material, only: &
    mappingConstitutive, &
    phase_vacancyInstance, &
    vacancyState
  use math, only: &
    math_Mandel6to33, &
-   math_trace33, &
+   math_trace33
-   pi
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -362,70 +310,29 @@ subroutine vacancy_generation_dotState(nSlip, accumulatedSlip, Tstar_v, Temperat
  real(pReal),  intent(in) :: &
    Temperature                                                                                      !< 2nd Piola Kirchhoff stress tensor (Mandel)
  real(pReal) :: &
-   pressure                                                                                         !< 2nd Piola Kirchhoff stress tensor (Mandel)
+   pressure, &
+   energyBarrier                                                                     
  integer(pInt) :: &
    instance, phase, constituent 
- real(pReal) :: &
-   vacancyConcentration, &                                                                          !< current vacancy concentration
-   vacancyDiffusion, &                                                                              !< the diffusion coefficient D_v
-   poleZeldovichCoeff, &                                                                            !< Zeldovich factor of pore nucleation
-   vacancyAbsorpRateCoeff, &                                                                        !< vacancy absorption rate
-   chemicalPotential, &                                                                             !< the chemical potential due to vacancy concentration
-   criticalRadius, &                                                                                !< the critical pore radius
-   Gibbs4Pore, &                                                                                    !< the Gibbs free energy for generating a critical pore
-   equilibPoreConcentration, &                                                                      !< the equilibrium pore concentration
-   nucleationRatePore, &                                                                            !< the nucleation rate of pore
-   ratioCvCve                                                                                       !< the ratio of Cv with respect to Cve
- real(pReal) :: &
-   threshold4ratioCvCve = 2.0_pReal                                                                 !< the threshold value for Cv/Cve
 
  phase = mappingConstitutive(2,ipc,ip,el)
  constituent = mappingConstitutive(1,ipc,ip,el)
  instance = phase_vacancyInstance(phase)
  pressure = math_trace33(math_Mandel6to33(Tstar_v))
+ energyBarrier = vacancy_generation_formationEnergy(instance) - &
+                 pressure*vacancy_generation_atomicVol(instance) - &
+                 sum(accumulatedSlip)*vacancy_generation_plasticityCoeff(instance)
 
-!--------------------------------------------------------------------------------------------------
- vacancyConcentration = vacancy_generation_getConcentration(ipc, ip, el)
- ratioCvCve = vacancyConcentration/vacancy_generation_equilibConcentration(instance)
-
- if(ratioCvCve < threshold4ratioCvCve) then
-   nucleationRatePore = 0.0_pReal
- else
-!  Calculate nucleation rate of pore
-   vacancyDiffusion = vacancy_generation_diffusionCoeff0(instance)* &
-                      exp( -vacancy_generation_diffusionEnergy(instance)/(kB*temperature) )
-   chemicalPotential = kB*Temperature * log(vacancyConcentration/ &
-                       vacancy_generation_equilibConcentration(instance))
-   criticalRadius = 2_pReal/chemicalPotential* &
-                    pore_nucleation_surfaceEnergy(instance) * pore_nucleation_atomVolume(instance)
-   Gibbs4Pore = 4_pReal/3_pReal * pi * pore_nucleation_surfaceEnergy(instance)* &
-                criticalRadius * criticalRadius
-   equilibPoreConcentration = pore_nucleation_concentrationCoeff0(instance)* &
-                              exp( -Gibbs4Pore/(kB*temperature) )
-
-   vacancyAbsorpRateCoeff = 2_pReal/pore_nucleation_shellThickness(instance) * &
-                            vacancyDiffusion * vacancyConcentration
-   poleZeldovichCoeff = pore_nucleation_atomVolume(instance)* &
-                        sqrt( pore_nucleation_surfaceEnergy(instance)/(kB*temperature) )
-   nucleationRatePore = poleZeldovichCoeff * vacancyAbsorpRateCoeff* equilibPoreConcentration
- endif
-
-!--------------------------------------------------------------------------------------------------
-!  the net generating rate vacancy                            
  vacancyState(phase)%dotState(1,constituent) = &
     vacancy_generation_freq(instance)* &
-    exp(-(vacancy_generation_formationEnergy(instance) - vacancy_generation_stressCoeff(instance)*pressure)/ &
+    exp(-energyBarrier/(kB*Temperature))
-         (kB*Temperature)) + &
-    sum(accumulatedSlip) * vacancy_generation_dislocationCoeff(instance)- &                         !< Induced by dislocation motion
-    nucleationRatePore * (4_pReal/3_pReal * pi * criticalRadius**3_pReal)/ &                        !< Reduced by the formation of pore
-    pore_nucleation_atomVolume(instance)   
 
 end subroutine vacancy_generation_dotState
 
 !--------------------------------------------------------------------------------------------------
 !> @brief returns vacancy concentration based on state layout 
 !--------------------------------------------------------------------------------------------------
-function vacancy_generation_getConcentration(ipc, ip, el)
+function vacancy_generation_getLocalConcentration(ipc, ip, el)
  use material, only: &
    mappingConstitutive, &
    vacancyState
@@ -435,17 +342,17 @@ function vacancy_generation_getConcentration(ipc, ip, el)
    ipc, &                                                                                           !< grain number
    ip, &                                                                                            !< integration point number
    el                                                                                               !< element number
- real(pReal) :: vacancy_generation_getConcentration
+ real(pReal) :: vacancy_generation_getLocalConcentration
  
- vacancy_generation_getConcentration = &
+ vacancy_generation_getLocalConcentration = &
    vacancyState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el))
  
-end function vacancy_generation_getConcentration
+end function vacancy_generation_getLocalConcentration
  
 !--------------------------------------------------------------------------------------------------
 !> @brief returns temperature based on local damage model state layout 
 !--------------------------------------------------------------------------------------------------
-subroutine vacancy_generation_putConcentration(ipc, ip, el, localVacancyConcentration)
+subroutine vacancy_generation_putLocalConcentration(ipc, ip, el, localVacancyConcentration)
  use material, only: &
    mappingConstitutive, &
    vacancyState
@@ -461,18 +368,68 @@ subroutine vacancy_generation_putConcentration(ipc, ip, el, localVacancyConcentr
  vacancyState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el))= &
    localVacancyConcentration
  
-end subroutine vacancy_generation_putConcentration
+end subroutine vacancy_generation_putLocalConcentration
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns vacancy concentration based on state layout 
+!--------------------------------------------------------------------------------------------------
+function vacancy_generation_getConcentration(ipc, ip, el)
+ use material, only: &
+   mappingHomogenization, &
+   material_phase, &
+   fieldVacancy, &
+   field_vacancy_type, &
+   FIELD_VACANCY_local_ID, &
+   FIELD_VACANCY_nonlocal_ID, &
+   material_homog
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal) :: vacancy_generation_getConcentration
+ 
+ select case(field_vacancy_type(material_homog(ip,el)))                                                   
+   case (FIELD_VACANCY_local_ID)
+    vacancy_generation_getConcentration = vacancy_generation_getLocalConcentration(ipc, ip, el)      
+    
+   case (FIELD_VACANCY_nonlocal_ID)
+    vacancy_generation_getConcentration = fieldVacancy(material_homog(ip,el))% &
+      field(1,mappingHomogenization(1,ip,el))                                                     ! Taylor type 
+ end select
+
+end function vacancy_generation_getConcentration
  
 !--------------------------------------------------------------------------------------------------
 !> @brief returns generation vacancy diffusion tensor 
 !--------------------------------------------------------------------------------------------------
-function vacancy_generation_getVacancyDiffusion33(nSlip,accumulatedSlip,temperature,ipc,ip,el)
+function vacancy_generation_getVacancyDiffusion33(ipc,ip,el)
  use lattice, only: &
    lattice_VacancyDiffusion33
+ use material, only: &
+   mappingConstitutive
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), dimension(3,3) :: &
+   vacancy_generation_getVacancyDiffusion33
+ 
+ vacancy_generation_getVacancyDiffusion33 = &
+   lattice_VacancyDiffusion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el))
+    
+end function vacancy_generation_getVacancyDiffusion33
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns generation vacancy mobility tensor 
+!--------------------------------------------------------------------------------------------------
+function vacancy_generation_getVacancyMobility33(nSlip,accumulatedSlip,temperature,ipc,ip,el)
  use material, only: &
    mappingConstitutive, &
-   phase_vacancyInstance, &
+   phase_vacancyInstance
-   vacancyState
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -481,7 +438,7 @@ function vacancy_generation_getVacancyDiffusion33(nSlip,accumulatedSlip,temperat
    ip, &                                                                                            !< integration point number
    el                                                                                               !< element number
  real(pReal), dimension(3,3) :: &
-   vacancy_generation_getVacancyDiffusion33
+   vacancy_generation_getVacancyMobility33
  real(pReal), dimension(nSlip) :: &
    accumulatedSlip
  real(pReal) :: &
@@ -493,11 +450,13 @@ function vacancy_generation_getVacancyDiffusion33(nSlip,accumulatedSlip,temperat
  constituent = mappingConstitutive(1,ipc,ip,el)
  instance = phase_vacancyInstance(phase)
 
- vacancy_generation_getVacancyDiffusion33 = &
+ vacancy_generation_getVacancyMobility33 = &
-   lattice_VacancyDiffusion33(1:3,1:3,phase)* &
+   vacancy_generation_surfaceEnergy(instance)* &
-   exp(-vacancy_generation_diffusionEnergy(instance)/(kB*temperature))
+   vacancy_generation_diffusionCoeff0(instance)* &
+   exp(-vacancy_generation_migrationEnergy(instance)/(kB*temperature))/ &
+   (kB*temperature)
     
-end function vacancy_generation_getVacancyDiffusion33
+end function vacancy_generation_getVacancyMobility33
 
 !--------------------------------------------------------------------------------------------------
 !> @brief return array of constitutive results