Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development

VERSION

@@ -1 +1 @@
-v2.0.1-1035-gd80a255
+v2.0.1-1054-gb834b2a

lib/damask/util.py

@@ -39,7 +39,7 @@ def srepr(arg,glue = '\n'):
   if (not hasattr(arg, "strip") and
           hasattr(arg, "__getitem__") or
           hasattr(arg, "__iter__")):
-     return glue.join(srepr(x) for x in arg)
+     return glue.join(str(x) for x in arg)
   return arg if isinstance(arg,str) else repr(arg)
 
 # -----------------------------
@@ -233,6 +233,7 @@ def leastsqBound(func, x0, args=(), bounds=None, Dfun=None, full_output=0,
  
   def _check_func(checker, argname, thefunc, x0, args, numinputs,
                 output_shape=None):
+    from numpy import shape
     """The same as that of minpack.py"""
     res = np.atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     if (output_shape is not None) and (shape(res) != output_shape):

src/DAMASK_spectral.f90

@@ -78,7 +78,8 @@ program DAMASK_spectral
    FIELD_UNDEFINED_ID, &
    FIELD_MECH_ID, &
    FIELD_THERMAL_ID, &
-   FIELD_DAMAGE_ID
+   FIELD_DAMAGE_ID, &
+   utilities_calcPlasticity
  use spectral_mech_Basic
  use spectral_mech_AL
  use spectral_mech_Polarisation
@@ -156,6 +157,19 @@ program DAMASK_spectral
    MPI_finalize, &
    MPI_allreduce, &
    PETScFinalize
+!--------------------------------------------------------------------------------------------------
+! variables related to stop criterion for yielding
+ real(pReal) :: plasticWorkOld, plasticWorkNew, &                                                   ! plastic work
+ eqTotalStrainOld,  eqTotalStrainNew, &                                                             ! total equivalent strain
+ eqPlasticStrainOld, eqPlasticStrainNew, &                                                          ! total equivalent plastic strain
+ eqStressOld, eqStressNew , &                                                                       ! equivalent stress
+ yieldStopValue
+ real(pReal), dimension(3,3) :: yieldStress,yieldStressOld,yieldStressNew, &
+   plasticStrainOld, plasticStrainNew, plasticStrainRate
+ integer(pInt) :: yieldResUnit = 0_pInt
+ integer(pInt) :: stressstrainUnit = 0_pInt
+ character(len=13)  :: stopFlag
+ logical :: yieldStop, yieldStopSatisfied
 
 !--------------------------------------------------------------------------------------------------
 ! init DAMASK (all modules)
@@ -213,6 +227,8 @@ program DAMASK_spectral
 
 !--------------------------------------------------------------------------------------------------
 ! reading the load case and assign values to the allocated data structure
+ yieldStop = .False.
+ yieldStopSatisfied = .False.
  rewind(FILEUNIT)
  do
    line = IO_read(FILEUNIT)
@@ -287,10 +303,30 @@ program DAMASK_spectral
            temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j)
          enddo
          loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector)
+       case('totalstrain')
+         yieldStop = .True.
+         stopFlag = 'totalStrain'
+         yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt)
+       case('plasticstrain')
+         yieldStop = .True.
+         stopFlag = 'plasticStrain'
+         yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt)
+       case('plasticwork')
+         yieldStop = .True.
+         stopFlag = 'plasticWork'
+         yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt)
      end select
  enddo; enddo
  close(FILEUNIT)
-
+ 
+ if(yieldStop) then                                                                                 ! initialize variables related to yield stop
+   yieldStressNew = 0.0_pReal
+   plasticStrainNew = 0.0_pReal
+   eqStressNew = 0.0_pReal
+   eqTotalStrainNew = 0.0_pReal
+   eqPlasticStrainNew = 0.0_pReal
+   plasticWorkNew = 0.0_pReal
+ endif
 !--------------------------------------------------------------------------------------------------
 ! consistency checks and output of load case
  loadCases(1)%followFormerTrajectory = .false.                                                      ! cannot guess along trajectory for first inc of first currentLoadCase
@@ -661,10 +697,76 @@ program DAMASK_spectral
        time = time + timeinc
        guess = .true.
      endif forwarding
+     
+     yieldCheck: if(yieldStop) then                                                                 ! check if it yields or satisfies the certain stop condition
+       yieldStressOld = yieldStressNew
+       plasticStrainOld = plasticStrainNew
+       eqStressOld = eqStressNew
+       eqTotalStrainOld = eqTotalStrainNew
+       eqPlasticStrainOld = eqPlasticStrainNew
+       plasticWorkOld = plasticWorkNew
+       
+       call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, &
+                                     eqPlasticStrainNew, plasticWorkNew, loadCases(currentLoadCase)%rotation)
+       
+       if (worldrank == 0) then                                                                     ! output the stress-strain curve to file if yield stop criterion is used
+         if ((currentLoadCase == 1_pInt) .and. (inc == 1_pInt)) then
+           open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
+                                  '.stressstrain',form='FORMATTED',status='REPLACE')
+           write(stressstrainUnit,*) 0.0_pReal, 0.0_pReal
+           write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew
+           close(stressstrainUnit)
+         else
+           open(newunit=stressstrainUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
+                                 '.stressstrain',form='FORMATTED', position='APPEND', status='OLD')
+           write(stressstrainUnit,*) eqTotalStrainNew, eqStressNew
+           close(stressstrainUnit)
+         endif
+       endif
+       
+       if(stopFlag == 'totalStrain') then
+         if(eqTotalStrainNew > yieldStopValue) then
+           yieldStress = yieldStressOld * (eqTotalStrainNew - yieldStopValue)/(eqTotalStrainNew - eqTotalStrainOld) &   ! linear interpolation of stress values
+                       + yieldStressNew * (yieldStopValue - eqTotalStrainOld)/(eqTotalStrainNew - eqTotalStrainOld)
+           plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0)                                     ! calculate plastic strain rate
+           yieldStopSatisfied = .True.
+         endif
+       elseif(stopFlag == 'plasticStrain') then
+         if(eqPlasticStrainNew > yieldStopValue) then
+           yieldStress = yieldStressOld * (eqPlasticStrainNew - yieldStopValue)/(eqPlasticStrainNew - eqPlasticStrainOld) &
+                       + yieldStressNew * (yieldStopValue - eqPlasticStrainOld)/(eqPlasticStrainNew - eqPlasticStrainOld)
+           plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0)
+           yieldStopSatisfied = .True.
+         endif
+       elseif(stopFlag == 'plasticWork') then
+         if(plasticWorkNew > yieldStopValue) then
+           yieldStress = yieldStressOld * (plasticWorkNew - yieldStopValue)/(plasticWorkNew - plasticWorkOld) &
+                       + yieldStressNew * (yieldStopValue - plasticWorkOld)/(plasticWorkNew - plasticWorkOld)
+           plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0)
+           yieldStopSatisfied = .True.
+         endif
+       endif
+     endif yieldCheck
+
+     if (yieldStopSatisfied) then                                                                   ! when yield, write the yield stress and strain rate to file and quit the job
+       if (worldrank == 0) then
+         open(newunit=yieldResUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//&
+                                  '.yield',form='FORMATTED',status='REPLACE')
+         do i = 1_pInt,3_pInt
+           write(yieldResUnit,*) (yieldStress(i,j), j=1,3)
+         enddo
+         do i = 1_pInt,3_pInt
+           write(yieldResUnit,*) (plasticStrainRate(i,j), j=1,3)
+         enddo
+         close(yieldResUnit)
+         call quit(0_pInt)
+       endif
+     endif
 
     enddo incLooping
  enddo loadCaseLooping
-
+ 
+ 
 !--------------------------------------------------------------------------------------------------
 ! report summary of whole calculation
  write(6,'(/,a)') ' ###########################################################################'

src/crystallite.f90

@@ -55,14 +55,14 @@ module crystallite
    crystallite_Li0, &                                                                               !< intermediate velocitiy grad at start of FE inc
    crystallite_partionedLi0,&                                                                       !< intermediate velocity grad at start of homog inc
    crystallite_Fe, &                                                                                !< current "elastic" def grad (end of converged time step)
-   crystallite_P                                                                                    !< 1st Piola-Kirchhoff stress per grain
+   crystallite_P,  &                                                                                !< 1st Piola-Kirchhoff stress per grain
+   crystallite_subF                                                                                 !< def grad to be reached at end of crystallite inc
  real(pReal),                dimension(:,:,:,:,:),    allocatable, private :: &
    crystallite_subFe0,&                                                                             !< "elastic" def grad at start of crystallite inc
    crystallite_invFp, &                                                                             !< inverse of current plastic def grad (end of converged time step)
    crystallite_subFp0,&                                                                             !< plastic def grad at start of crystallite inc
    crystallite_invFi, &                                                                             !< inverse of current intermediate def grad (end of converged time step)
    crystallite_subFi0,&                                                                             !< intermediate def grad at start of crystallite inc
-   crystallite_subF,  &                                                                             !< def grad to be reached at end of crystallite inc
    crystallite_subF0, &                                                                             !< def grad at start of crystallite inc
    crystallite_subLp0,&                                                                             !< plastic velocity grad at start of crystallite inc
    crystallite_subLi0,&                                                                             !< intermediate velocity grad at start of crystallite inc

src/math.f90

@@ -144,6 +144,7 @@ module math
    math_sampleGaussVar, &
    math_symmetricEulers, &
    math_eigenvectorBasisSym33, &
+   math_eigenvectorBasisSym33_log, &
    math_eigenvectorBasisSym, &
    math_eigenValuesVectorsSym33, &
    math_eigenValuesVectorsSym, &
@@ -2117,6 +2118,70 @@ function math_eigenvectorBasisSym33(m)
 end function math_eigenvectorBasisSym33
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief logarithm eigenvector basis of symmetric 33 matrix m
+!--------------------------------------------------------------------------------------------------
+function math_eigenvectorBasisSym33_log(m)
+
+ implicit none
+ real(pReal), dimension(3,3)              :: math_eigenvectorBasisSym33_log
+ real(pReal), dimension(3)                :: invariants, values
+ real(pReal), dimension(3,3), intent(in) :: m
+ real(pReal) :: P, Q, rho, phi
+ real(pReal), parameter :: TOL=1.e-14_pReal
+ real(pReal), dimension(3,3,3) :: N, EB
+
+ invariants = math_invariantsSym33(m)
+ EB = 0.0_pReal
+
+ P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal
+ Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3)
+
+ threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then
+   values = invariants(1)/3.0_pReal
+!   this is not really correct, but at least the basis is correct
+   EB(1,1,1)=1.0_pReal
+   EB(2,2,2)=1.0_pReal
+   EB(3,3,3)=1.0_pReal
+ else threeSimilarEigenvalues
+   rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal
+   phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal))
+   values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* &
+                             [cos(phi/3.0_pReal), &
+                              cos((phi+2.0_pReal*PI)/3.0_pReal), &
+                              cos((phi+4.0_pReal*PI)/3.0_pReal) &
+                             ] + invariants(1)/3.0_pReal
+   N(1:3,1:3,1) = m-values(1)*math_I3
+   N(1:3,1:3,2) = m-values(2)*math_I3
+   N(1:3,1:3,3) = m-values(3)*math_I3
+   twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then 
+     EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
+                                               ((values(3)-values(1))*(values(3)-values(2)))
+     EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
+   elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
+     EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
+                                               ((values(1)-values(2))*(values(1)-values(3)))
+     EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
+   elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues 
+     EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
+                                               ((values(2)-values(1))*(values(2)-values(3)))
+     EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
+   else twoSimilarEigenvalues
+     EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ &
+                                               ((values(1)-values(2))*(values(1)-values(3)))
+     EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ &
+                                               ((values(2)-values(1))*(values(2)-values(3)))
+     EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ &
+                                               ((values(3)-values(1))*(values(3)-values(2)))
+   endif twoSimilarEigenvalues
+ endif threeSimilarEigenvalues
+
+ math_eigenvectorBasisSym33_log = log(sqrt(values(1))) * EB(1:3,1:3,1) &
+                            + log(sqrt(values(2))) * EB(1:3,1:3,2) &
+                            + log(sqrt(values(3))) * EB(1:3,1:3,3)
+
+end function math_eigenvectorBasisSym33_log
+
 !--------------------------------------------------------------------------------------------------
 !> @brief rotational part from polar decomposition of 33 tensor m
 !--------------------------------------------------------------------------------------------------

src/spectral_utilities.f90

@@ -145,7 +145,8 @@ module spectral_utilities
    FIELD_UNDEFINED_ID, &
    FIELD_MECH_ID, &
    FIELD_THERMAL_ID, &
-   FIELD_DAMAGE_ID
+   FIELD_DAMAGE_ID, &
+   utilities_calcPlasticity
  private :: &
    utilities_getFreqDerivative
 
@@ -1047,6 +1048,125 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
 
 end subroutine utilities_constitutiveResponse
 
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates yield stress, plastic strain, total strain and their equivalent values
+!--------------------------------------------------------------------------------------------------
+subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTotalStrain, &
+                                    eqPlasticStrain, plasticWork, rotation_BC)
+ use crystallite, only: &
+   crystallite_Fe, &
+   crystallite_P, &
+   crystallite_subF
+ use material, only: &
+   homogenization_maxNgrains
+ use mesh, only: &
+   mesh_maxNips,&
+   mesh_NcpElems
+ use math, only: &
+   math_det33, &
+   math_inv33, &
+   math_mul33x33, &
+   math_trace33, &
+   math_transpose33, &
+   math_equivStrain33, &
+   math_equivStress33, &
+   math_rotate_forward33, &
+   math_identity2nd, &
+   math_crossproduct, &
+   math_eigenvectorBasisSym, &
+   math_eigenvectorBasisSym33, &
+   math_eigenvectorBasisSym33_log, &
+   math_eigenValuesVectorsSym33
+   
+ implicit none
+
+ real(pReal), intent(inout) :: eqStress, eqPlasticStrain, plasticWork
+ real(pReal), intent(out) :: eqTotalStrain
+ real(pReal), dimension(3,3),intent(out) :: yieldStress, plasticStrain
+ real(pReal), intent(in), dimension(3,3) :: rotation_BC                                             !< rotation of load frame
+ real(pReal), dimension(3,3) :: cauchy, P_av, F_av, Ve_av                                           !< average
+ real(pReal), dimension(3)   :: Values, S
+ real(pReal), dimension(3,3) :: Vectors, diag
+ real(pReal), dimension(3,3) :: &
+   Vp, F_temp, U, VT, R, V, V_total
+ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+   Be, Ve, Fe
+ real(pReal), dimension(15)   :: WORK                                                               !< previous deformation gradient
+ integer(pInt) :: INFO, i, j, k, l, ierr
+ real(pReal) :: wgtm
+ real(pReal) :: eqStressOld, eqPlasticStrainOld, plasticWorkOld
+ 
+ external :: dgesvd
+ 
+ eqStressOld = eqStress
+ eqPlasticStrainOld = eqPlasticStrain
+ plasticWorkOld = plasticWork
+ wgtm = 1.0_pReal/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal)
+ diag = 0.0_pReal
+ 
+ P_av = sum(sum(sum(crystallite_P,dim=5),dim=4),dim=3) * wgtm
+ call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ P_av = math_rotate_forward33(P_av,rotation_BC)
+ 
+ F_av = sum(sum(sum(crystallite_subF,dim=5),dim=4),dim=3) * wgtm
+ call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ F_av = math_rotate_forward33(F_av,rotation_BC)
+ 
+ cauchy = 1.0_pReal/math_det33(F_av)*math_mul33x33(P_av,transpose(F_av))
+ yieldStress = cauchy
+ eqStress = math_equivStress33(cauchy)
+ 
+ F_temp = F_av
+ call dgesvd ('A', 'A', 3, 3, F_temp, 3, S, U, 3, VT, 3, WORK, 15, INFO)                            ! singular value decomposition
+ 
+ R = math_mul33x33(U, VT)                                                                           ! rotation of polar decomposition
+ V = math_mul33x33(F_av,math_inv33(R))
+ 
+ call math_eigenValuesVectorsSym33(V,Values,Vectors)
+ do l = 1_pInt, 3_pInt
+   if (Values(l) < 0.0_pReal) then
+     Values(l) = -Values(l)
+     Vectors(1:3, l) = -Vectors(1:3, l)
+   endif
+   Values(l) = log(Values(l))
+   diag(l,l) = Values(l)
+ enddo
+ if (dot_product(Vectors(1:3,1),Vectors(1:3,2)) /= 0) then
+   Vectors(1:3,2) = math_crossproduct(Vectors(1:3,3), Vectors(1:3,1))
+   Vectors(1:3,2) = Vectors(1:3,2)/sqrt(dot_product(Vectors(1:3,2),Vectors(1:3,2)))
+ endif
+ if (dot_product(Vectors(1:3,2),Vectors(1:3,3)) /= 0) then
+   Vectors(1:3,3) = math_crossproduct(Vectors(1:3,1), Vectors(1:3,2))
+   Vectors(1:3,3) = Vectors(1:3,3)/sqrt(dot_product(Vectors(1:3,3),Vectors(1:3,3)))
+ endif
+ if (dot_product(Vectors(1:3,3),Vectors(1:3,1)) /= 0) then
+   Vectors(1:3,1) = math_crossproduct(Vectors(1:3,2), Vectors(1:3,3))
+   Vectors(1:3,1) = Vectors(1:3,1)/sqrt(dot_product(Vectors(1:3,1),Vectors(1:3,1)))
+ endif
+ 
+ V_total = REAL(math_mul33x33(Vectors, math_mul33x33(diag, transpose(Vectors))))
+ eqTotalStrain = math_equivStrain33(V_total)
+ 
+ do k = 1_pInt, mesh_NcpElems;  do j = 1_pInt, mesh_maxNips;  do i = 1_pInt,homogenization_maxNgrains
+   Fe(1:3,1:3,i,j,k) = crystallite_Fe(1:3,1:3,i,j,k)
+   Fe(1:3,1:3,i,j,k) = math_rotate_forward33(Fe(1:3,1:3,i,j,k),rotation_BC)
+   Be(1:3,1:3,i,j,k) = math_mul33x33(Fe(1:3,1:3,i,j,k),math_transpose33(Fe(1:3,1:3,i,j,k)))         ! elastic part of left Cauchy–Green deformation tensor
+   Ve(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Be(1:3,1:3,i,j,k))
+ enddo; enddo; enddo
+ 
+ Ve_av = sum(sum(sum(Ve,dim=5),dim=4),dim=3) * wgtm
+ call MPI_Allreduce(MPI_IN_PLACE,Ve_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
+ 
+ Vp = V_total - Ve_av
+ 
+ eqPlasticStrain = math_equivStrain33(Vp)
+ 
+ plasticStrain = Vp
+ 
+ plasticWork = plasticWorkOld + 0.5*(eqStressOld + eqStress) * (eqPlasticStrain - eqPlasticStrainOld)
+ 
+end subroutine utilities_calcPlasticity
+
 
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates forward rate, either guessing or just add delta/timeinc