new variables for damage

PRIVATE

@@ -1 +1 @@
-Subproject commit 9e8011876b16896a18902737790eba8018e94f85
+Subproject commit 5fd23c60c721aada27bcd3a4ba96301753b7bd0d

src/commercialFEM_fileList.f90

@@ -52,4 +52,5 @@
 #include "homogenization_mech_isostrain.f90"
 #include "homogenization_mech_RGC.f90"
 #include "homogenization_thermal.f90"
+#include "homogenization_damage.f90"
 #include "CPFEM.f90"

src/constitutive.f90

@@ -188,6 +188,11 @@ module constitutive
       real(pReal), dimension(3,3) :: P
     end function constitutive_mech_getP
 
+    module function constitutive_damage_get_phi(co,ip,el) result(phi)
+      integer, intent(in) :: co, ip, el
+      real(pReal) :: phi
+    end function constitutive_damage_get_phi
+
     module function thermal_T(ph,me) result(T)
       integer, intent(in) :: ph,me
       real(pReal) :: T
@@ -203,6 +208,11 @@ module constitutive
       real(pReal), intent(in) :: T
       integer, intent(in) :: co, ce
     end subroutine constitutive_thermal_setT
+    
+    module subroutine constitutive_damage_set_phi(phi,co,ce)
+      real(pReal), intent(in) :: phi
+      integer, intent(in) :: co, ce
+    end subroutine constitutive_damage_set_phi
 
 ! == cleaned:end ===================================================================================
 
@@ -353,6 +363,8 @@ module constitutive
     constitutive_restartRead, &
     integrateDamageState, &
     constitutive_thermal_setT, &
+    constitutive_damage_set_phi, &
+    constitutive_damage_get_phi, &
     constitutive_mech_getP, &
     constitutive_mech_setF, &
     constitutive_mech_getF, &

src/constitutive_damage.f90

@@ -10,9 +10,16 @@ submodule(constitutive) constitutive_damage
     DAMAGE_ANISODUCTILE_ID
   end enum
 
+
+  type :: tDataContainer
+    real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi
+  end type tDataContainer
+
   integer(kind(DAMAGE_UNDEFINED_ID)),     dimension(:,:), allocatable :: &
     phase_source                                                                                !< active sources mechanisms of each phase
 
+  type(tDataContainer), dimension(:), allocatable :: current
+
   interface
 
   module function source_damage_anisoBrittle_init(source_length) result(mySources)
@@ -152,7 +159,8 @@ contains
 module subroutine damage_init
 
   integer :: &
-    ph                                                                                              !< counter in phase loop
+    ph, &                                                                                           !< counter in phase loop
+    Nconstituents
   class(tNode), pointer :: &
    phases, &
    phase, &
@@ -161,10 +169,18 @@ module subroutine damage_init
 
   phases => config_material%get('phase')
 
+  allocate(current(phases%length))
+
   allocate(damageState (phases%length))
   allocate(phase_Nsources(phases%length),source = 0)           ! same for kinematics
 
   do ph = 1,phases%length
+
+    Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
+
+    allocate(current(ph)%phi(Nconstituents),source=1.0_pReal)
+    allocate(current(ph)%d_phi_d_dot_phi(Nconstituents),source=0.0_pReal)
+
     phase => phases%get(ph)
     sources => phase%get('source',defaultVal=emptyList)
     phase_Nsources(ph) = sources%length
@@ -511,4 +527,28 @@ function source_active(source_label,src_length)  result(active_source)
 end function source_active
 
 
+!----------------------------------------------------------------------------------------------
+!< @brief Set damage parameter
+!----------------------------------------------------------------------------------------------
+module subroutine constitutive_damage_set_phi(phi,co,ce)
+
+  real(pReal), intent(in) :: phi
+  integer, intent(in) :: ce, co
+
+
+  current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi
+
+end subroutine constitutive_damage_set_phi
+
+
+module function constitutive_damage_get_phi(co,ip,el) result(phi)
+  
+  integer, intent(in) :: co, ip, el
+  real(pReal) :: phi
+
+  phi = current(material_phaseAt(co,el))%phi(material_phaseMemberAt(co,ip,el))
+
+end function constitutive_damage_get_phi
+
+
 end submodule constitutive_damage

src/grid/grid_damage_spectral.f90

@@ -17,8 +17,9 @@ module grid_damage_spectral
   use discretization_grid
   use damage_nonlocal
   use YAML_types
+  use homogenization
   use config
- 
+
   implicit none
   private
 
@@ -43,13 +44,13 @@ module grid_damage_spectral
     phi_current, &                                                                                  !< field of current damage
     phi_lastInc, &                                                                                  !< field of previous damage
     phi_stagInc                                                                                     !< field of staggered damage
- 
+
 !--------------------------------------------------------------------------------------------------
-! reference diffusion tensor, mobility etc. 
+! reference diffusion tensor, mobility etc.
   integer                     :: totalIter = 0                                                      !< total iteration in current increment
   real(pReal), dimension(3,3) :: K_ref
   real(pReal)                 :: mu_ref
-  
+
   public :: &
     grid_damage_spectral_init, &
     grid_damage_spectral_solution, &
@@ -62,8 +63,8 @@ contains
 ! ToDo: Restart not implemented
 !--------------------------------------------------------------------------------------------------
 subroutine grid_damage_spectral_init
-    
+
-  PetscInt, dimension(0:worldsize-1) :: localK  
+  PetscInt, dimension(0:worldsize-1) :: localK
   DM :: damage_grid
   Vec :: uBound, lBound
   PetscErrorCode :: ierr
@@ -72,22 +73,22 @@ subroutine grid_damage_spectral_init
     num_generic
   character(len=pStringLen) :: &
     snes_type
- 
+
   print'(/,a)', ' <<<+-  grid_spectral_damage init  -+>>>'
 
   print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019'
   print*, 'https://doi.org/10.1007/978-981-10-6855-3_80'
- 
+
 !-------------------------------------------------------------------------------------------------
 ! read numerical parameters and do sanity checks
   num_grid => config_numerics%get('grid',defaultVal=emptyDict)
   num%itmax           = num_grid%get_asInt   ('itmax',defaultVal=250)
   num%eps_damage_atol = num_grid%get_asFloat ('eps_damage_atol',defaultVal=1.0e-2_pReal)
   num%eps_damage_rtol = num_grid%get_asFloat ('eps_damage_rtol',defaultVal=1.0e-6_pReal)
- 
+
   num_generic => config_numerics%get('generic',defaultVal=emptyDict)
   num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal)
- 
+
   if (num%residualStiffness < 0.0_pReal)   call IO_error(301,ext_msg='residualStiffness')
   if (num%itmax <= 1)                      call IO_error(301,ext_msg='itmax')
   if (num%eps_damage_atol <= 0.0_pReal)    call IO_error(301,ext_msg='eps_damage_atol')
@@ -104,7 +105,7 @@ subroutine grid_damage_spectral_init
 !--------------------------------------------------------------------------------------------------
 ! initialize solver specific parts of PETSc
   call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr)
-  call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) 
+  call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr)
   localK            = 0
   localK(worldrank) = grid3
   call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
@@ -122,7 +123,7 @@ subroutine grid_damage_spectral_init
   call DMsetUp(damage_grid,ierr); CHKERRQ(ierr)
   call DMCreateGlobalVector(damage_grid,solution_vec,ierr); CHKERRQ(ierr)                           ! global solution vector (grid x 1, i.e. every def grad tensor)
   call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr)        ! residual vector of same shape as solution vector
-  CHKERRQ(ierr) 
+  CHKERRQ(ierr)
   call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr)                                          ! pull it all together with additional CLI arguments
   call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr)
   if (trim(snes_type) == 'vinewtonrsls' .or. &
@@ -137,12 +138,12 @@ subroutine grid_damage_spectral_init
   endif
 
 !--------------------------------------------------------------------------------------------------
-! init fields             
+! init fields
   call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr)
   CHKERRQ(ierr)
   xend = xstart + xend - 1
   yend = ystart + yend - 1
-  zend = zstart + zend - 1 
+  zend = zstart + zend - 1
   allocate(phi_current(grid(1),grid(2),grid3), source=1.0_pReal)
   allocate(phi_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
   allocate(phi_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
@@ -157,26 +158,26 @@ end subroutine grid_damage_spectral_init
 !> @brief solution for the spectral damage scheme with internal iterations
 !--------------------------------------------------------------------------------------------------
 function grid_damage_spectral_solution(timeinc) result(solution)
- 
+
   real(pReal), intent(in) :: &
     timeinc                                                                                         !< increment in time for current solution
   integer :: i, j, k, cell
   type(tSolutionState) :: solution
   PetscInt  :: devNull
   PetscReal :: phi_min, phi_max, stagNorm, solnNorm
-  
+
   PetscErrorCode :: ierr
   SNESConvergedReason :: reason
 
   solution%converged =.false.
- 
+
 !--------------------------------------------------------------------------------------------------
-! set module wide availabe data 
+! set module wide availabe data
   params%timeinc = timeinc
- 
+
   call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
   call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)
- 
+
   if (reason < 1) then
     solution%converged = .false.
     solution%iterationsNeeded = num%itmax
@@ -192,20 +193,21 @@ function grid_damage_spectral_solution(timeinc) result(solution)
   solution%stagConverged = stagNorm < max(num%eps_damage_atol, num%eps_damage_rtol*solnNorm)
 
 !--------------------------------------------------------------------------------------------------
-! updating damage state 
+! updating damage state
   cell = 0
   do k = 1, grid3;  do j = 1, grid(2);  do i = 1,grid(1)
     cell = cell + 1
     call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
+    homogenization_phi(cell) = phi_current(i,j,k)
   enddo; enddo; enddo
- 
+
   call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr)
   call VecMax(solution_vec,devNull,phi_max,ierr); CHKERRQ(ierr)
   if (solution%converged) &
     print'(/,a)', ' ... nonlocal damage converged .....................................'
   print'(/,a,f8.6,2x,f8.6,2x,e11.4)', ' Minimum|Maximum|Delta Damage      = ', phi_min, phi_max, stagNorm
   print'(/,a)', ' ==========================================================================='
-  flush(IO_STDOUT) 
+  flush(IO_STDOUT)
 
 end function grid_damage_spectral_solution
 
@@ -214,31 +216,32 @@ end function grid_damage_spectral_solution
 !> @brief spectral damage forwarding routine
 !--------------------------------------------------------------------------------------------------
 subroutine grid_damage_spectral_forward(cutBack)
-   
+
   logical, intent(in) :: cutBack
   integer                               :: i, j, k, cell
   DM :: dm_local
   PetscScalar,  dimension(:,:,:), pointer     :: x_scal
   PetscErrorCode                              :: ierr
 
-  if (cutBack) then 
+  if (cutBack) then
     phi_current = phi_lastInc
     phi_stagInc = phi_lastInc
 !--------------------------------------------------------------------------------------------------
-! reverting damage field state 
+! reverting damage field state
     cell = 0
     call SNESGetDM(damage_snes,dm_local,ierr); CHKERRQ(ierr)
     call DMDAVecGetArrayF90(dm_local,solution_vec,x_scal,ierr); CHKERRQ(ierr)                       !< get the data out of PETSc to work with
     x_scal(xstart:xend,ystart:yend,zstart:zend) = phi_current
     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)
-      cell = cell + 1                                                                           
+      cell = cell + 1
       call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell)
+      homogenization_phi(cell) = phi_current(i,j,k)
     enddo; enddo; enddo
   else
     phi_lastInc = phi_current
     call updateReference
-  endif  
+  endif
 
 end subroutine grid_damage_spectral_forward
 
@@ -247,7 +250,7 @@ end subroutine grid_damage_spectral_forward
 !> @brief forms the spectral damage residual vector
 !--------------------------------------------------------------------------------------------------
 subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
- 
+
   DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
     in
   PetscScalar, dimension( &
@@ -261,11 +264,11 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
   integer :: i, j, k, cell
   real(pReal)   :: phiDot, dPhiDot_dPhi, mobility
 
-  phi_current = x_scal 
+  phi_current = x_scal
 !--------------------------------------------------------------------------------------------------
 ! evaluate polarization field
   scalarField_real = 0.0_pReal
-  scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current 
+  scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_current
   call utilities_FFTscalarForward
   call utilities_fourierScalarGradient                                                              !< calculate gradient of damage field
   call utilities_FFTvectorBackward
@@ -297,7 +300,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
         scalarField_real(1:grid(1),1:grid(2),1:grid3) = phi_lastInc
   where(scalarField_real(1:grid(1),1:grid(2),1:grid3) < num%residualStiffness) &
         scalarField_real(1:grid(1),1:grid(2),1:grid3) = num%residualStiffness
- 
+
 !--------------------------------------------------------------------------------------------------
 ! constructing residual
   f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - phi_current
@@ -311,7 +314,7 @@ end subroutine formResidual
 subroutine updateReference
 
   integer :: i,j,k,cell,ierr
-  
+
   cell = 0
   K_ref = 0.0_pReal
   mu_ref = 0.0_pReal

src/homogenization.f90

@@ -64,6 +64,9 @@ module homogenization
     module subroutine thermal_init
     end subroutine thermal_init
 
+    module subroutine damage_init
+    end subroutine damage_init
+
     module subroutine mech_partition(subF,ip,el)
       real(pReal), intent(in), dimension(3,3) :: &
         subF
@@ -77,6 +80,11 @@ module homogenization
       integer,     intent(in) :: ce
     end subroutine thermal_partition
 
+    module subroutine damage_partition(phi,ce)
+      real(pReal), intent(in) :: phi
+      integer,     intent(in) :: ce
+    end subroutine damage_partition
+
     module subroutine thermal_homogenize(ip,el)
       integer, intent(in) :: ip,el
     end subroutine thermal_homogenize
@@ -120,7 +128,7 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief module initialization
 !--------------------------------------------------------------------------------------------------
-subroutine homogenization_init
+subroutine homogenization_init()
 
   class (tNode) , pointer :: &
     num_homog, &
@@ -144,6 +152,7 @@ subroutine homogenization_init
 
   call mech_init(num_homog)
   call thermal_init()
+  call damage_init()
 
   if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T)
   if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T)

src/homogenization_damage.f90

@@ -0,0 +1,40 @@
+!--------------------------------------------------------------------------------------------------
+!> @author Martin Diehl, KU Leuven
+!--------------------------------------------------------------------------------------------------
+submodule(homogenization) homogenization_damage
+
+
+contains
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Allocate variables and set parameters.
+!--------------------------------------------------------------------------------------------------
+module subroutine damage_init()
+
+
+  print'(/,a)',   ' <<<+-  homogenization_damage init  -+>>>'
+
+  allocate(homogenization_phi(discretization_nIPs*discretization_Nelems))
+  allocate(homogenization_dot_phi(discretization_nIPs*discretization_Nelems))
+
+end subroutine damage_init
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief Partition temperature onto the individual constituents.
+!--------------------------------------------------------------------------------------------------
+module subroutine damage_partition(phi,ce)
+
+  real(pReal), intent(in) :: phi
+  integer,     intent(in) :: ce
+
+  integer :: co
+
+  do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce))
+    call constitutive_damage_set_phi(phi,co,ce)
+  enddo
+
+end subroutine damage_partition
+
+
+end submodule homogenization_damage