isoductile is not working

PRIVATE

@@ -1 +1 @@
-Subproject commit fe88ce67a7c3deefe10a3a8a7eeab2215464bc76
+Subproject commit 122022609581777fdb323fc4b1b97d593f22bd58

src/commercialFEM_fileList.f90

@@ -32,14 +32,12 @@
 #include "phase_mechanical_plastic_nonlocal.f90"
 #include "phase_mechanical_eigen.f90"
 #include "phase_mechanical_eigen_cleavageopening.f90"
-#include "phase_mechanical_eigen_slipplaneopening.f90"
 #include "phase_mechanical_eigen_thermalexpansion.f90"
 #include "phase_thermal.f90"
 #include "phase_thermal_dissipation.f90"
 #include "phase_thermal_externalheat.f90"
 #include "phase_damage.f90"
 #include "phase_damage_isobrittle.f90"
-#include "phase_damage_isoductile.f90"
 #include "phase_damage_anisobrittle.f90"
 #include "homogenization.f90"
 #include "homogenization_mechanical.f90"

src/phase.f90

@@ -280,16 +280,6 @@ module phase
         dLd_dTstar                                                                                  !< derivative of Ld with respect to Tstar (4th-order tensor)
     end subroutine damage_anisobrittle_LiAndItsTangent
 
-    module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
-      integer, intent(in) :: ph, me
-      real(pReal),   intent(in),  dimension(3,3) :: &
-        S
-      real(pReal),   intent(out), dimension(3,3) :: &
-        Ld                                                                                          !< damage velocity gradient
-      real(pReal),   intent(out), dimension(3,3,3,3) :: &
-        dLd_dTstar                                                                                  !< derivative of Ld with respect to Tstar (4th-order tensor)
-    end subroutine damage_isoductile_LiAndItsTangent
-
   end interface
 
 

src/phase_damage.f90

@@ -11,7 +11,6 @@ submodule(phase) damage
   enum, bind(c); enumerator :: &
     DAMAGE_UNDEFINED_ID, &
     DAMAGE_ISOBRITTLE_ID, &
-    DAMAGE_ISODUCTILE_ID, &
     DAMAGE_ANISOBRITTLE_ID
   end enum
 
@@ -39,10 +38,6 @@ submodule(phase) damage
       logical, dimension(:), allocatable :: mySources
     end function isobrittle_init
 
-    module function isoductile_init() result(mySources)
-      logical, dimension(:), allocatable :: mySources
-    end function isoductile_init
-
 
     module subroutine isobrittle_deltaState(C, Fe, ph, me)
       integer, intent(in) :: ph,me
@@ -59,10 +54,6 @@ submodule(phase) damage
         S
     end subroutine anisobrittle_dotState
 
-    module subroutine isoductile_dotState(ph,me)
-      integer, intent(in) :: ph,me
-    end subroutine isoductile_dotState
-
     module subroutine anisobrittle_results(phase,group)
       integer,          intent(in) :: phase
       character(len=*), intent(in) :: group
@@ -73,11 +64,6 @@ submodule(phase) damage
       character(len=*), intent(in) :: group
     end subroutine isobrittle_results
 
-    module subroutine isoductile_results(phase,group)
-      integer,          intent(in) :: phase
-      character(len=*), intent(in) :: group
-    end subroutine isoductile_results
-
  end interface
 
 contains
@@ -131,7 +117,6 @@ module subroutine damage_init
 
   if (damage_active) then
     where(isobrittle_init()  ) phase_damage = DAMAGE_ISOBRITTLE_ID
-    where(isoductile_init()  ) phase_damage = DAMAGE_ISODUCTILE_ID
     where(anisobrittle_init()) phase_damage = DAMAGE_ANISOBRITTLE_ID
   endif
 
@@ -178,7 +163,7 @@ module function phase_f_phi(phi,co,ce) result(f)
   en = material_phaseEntry(co,ce)
 
   select case(phase_damage(ph))
-    case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID)
+    case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ANISOBRITTLE_ID)
       f = 1.0_pReal &
         - phi*damageState(ph)%state(1,en)
     case default
@@ -304,9 +289,6 @@ module subroutine damage_results(group,ph)
     case (DAMAGE_ISOBRITTLE_ID) sourceType
       call isobrittle_results(ph,group//'damage/')
 
-    case (DAMAGE_ISODUCTILE_ID) sourceType
-      call isoductile_results(ph,group//'damage/')
-
     case (DAMAGE_ANISOBRITTLE_ID) sourceType
       call anisobrittle_results(ph,group//'damage/')
 
@@ -332,9 +314,6 @@ function phase_damage_collectDotState(ph,me) result(broken)
 
     sourceType: select case (phase_damage(ph))
 
-      case (DAMAGE_ISODUCTILE_ID) sourceType
-        call isoductile_dotState(ph,me)
-
       case (DAMAGE_ANISOBRITTLE_ID) sourceType
         call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress?
 

src/phase_damage_isoductile.f90

@@ -1,127 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine incorporating isotropic ductile damage source mechanism
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-submodule(phase:damage) isoductile
-
-  type:: tParameters                                                                                !< container type for internal constitutive parameters
-    real(pReal) :: &
-      gamma_crit, &                                                                                 !< critical plastic strain
-      q
-    character(len=pStringLen), allocatable, dimension(:) :: &
-      output
-  end type tParameters
-
-  type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters (len Ninstances)
-
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-module function isoductile_init() result(mySources)
-
-  logical, dimension(:), allocatable :: mySources
-
-  class(tNode), pointer :: &
-    phases, &
-    phase, &
-    sources, &
-    src
-  integer :: Ninstances,Nmembers,ph
-  character(len=pStringLen) :: extmsg = ''
-
-
-  mySources = source_active('isoductile')
-  if(count(mySources) == 0) return
-
-  print'(/,a)', ' <<<+-  phase:damage:isoductile init  -+>>>'
-  print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
-
-
-  phases => config_material%get('phase')
-  allocate(param(phases%length))
-
-  do ph = 1, phases%length
-    if(mySources(ph)) then
-      phase => phases%get(ph)
-      sources => phase%get('damage')
-
-        associate(prm  => param(ph))
-        src => sources%get(1)
-
-        prm%q          = src%get_asFloat('q')
-        prm%gamma_crit = src%get_asFloat('gamma_crit')
-
-#if defined (__GFORTRAN__)
-        prm%output = output_as1dString(src)
-#else
-        prm%output = src%get_as1dString('output',defaultVal=emptyStringArray)
-#endif
-
-        ! sanity checks
-        if (prm%q          <= 0.0_pReal) extmsg = trim(extmsg)//' q'
-        if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
-
-        Nmembers=count(material_phaseID==ph)
-        call phase_allocateState(damageState(ph),Nmembers,1,1,0)
-        damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
-        if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
-
-        end associate
-
-!--------------------------------------------------------------------------------------------------
-!  exit if any parameter is out of range
-        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isoductile)')
-      endif
-  enddo
-
-
-end function isoductile_init
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief calculates derived quantities from state
-!--------------------------------------------------------------------------------------------------
-module subroutine isoductile_dotState(ph, me)
-
-  integer, intent(in) :: &
-    ph, &
-    me
-
-
-  associate(prm => param(ph))
-    damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)) &
-                                   / (prm%gamma_crit*damage_phi(ph,me)**prm%q)
-  end associate
-
-end subroutine isoductile_dotState
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief writes results to HDF5 output file
-!--------------------------------------------------------------------------------------------------
-module subroutine isoductile_results(phase,group)
-
-  integer,          intent(in) :: phase
-  character(len=*), intent(in) :: group
-
-  integer :: o
-
-  associate(prm => param(phase), stt => damageState(phase)%state)
-    outputsLoop: do o = 1,size(prm%output)
-      select case(trim(prm%output(o)))
-        case ('f_phi')
-          call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³')
-      end select
-    enddo outputsLoop
-  end associate
-
-end subroutine isoductile_results
-
-end submodule isoductile

src/phase_mechanical.f90

@@ -15,7 +15,6 @@ submodule(phase) mechanical
     PLASTICITY_NONLOCAL_ID, &
     KINEMATICS_UNDEFINED_ID, &
     KINEMATICS_CLEAVAGE_OPENING_ID, &
-    KINEMATICS_SLIPPLANE_OPENING_ID, &
     KINEMATICS_THERMAL_EXPANSION_ID
   end enum
 

src/phase_mechanical_eigen.f90

@@ -13,10 +13,6 @@ submodule(phase:mechanical) eigen
       logical, dimension(:), allocatable :: myKinematics
     end function damage_anisobrittle_init
 
-    module function damage_isoductile_init() result(myKinematics)
-      logical, dimension(:), allocatable :: myKinematics
-    end function damage_isoductile_init
-
     module function thermalexpansion_init(kinematics_length) result(myKinematics)
       integer, intent(in) :: kinematics_length
       logical, dimension(:,:), allocatable :: myKinematics
@@ -70,7 +66,6 @@ module subroutine eigendeformation_init(phases)
   allocate(model_damage(phases%length),  source = KINEMATICS_UNDEFINED_ID)
 
   where(damage_anisobrittle_init())  model_damage = KINEMATICS_cleavage_opening_ID
-  where(damage_isoductile_init())    model_damage = KINEMATICS_slipplane_opening_ID
 
 
 end subroutine eigendeformation_init
@@ -201,11 +196,6 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
       Li = Li + my_Li
       dLi_dS = dLi_dS + my_dLi_dS
       active = .true.
-    case (KINEMATICS_slipplane_opening_ID)
-      call damage_isoductile_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, en)
-      Li = Li + my_Li
-      dLi_dS = dLi_dS + my_dLi_dS
-      active = .true.
   end select
 
   if(.not. active) return

src/phase_mechanical_eigen_slipplaneopening.f90

@@ -1,184 +0,0 @@
-!--------------------------------------------------------------------------------------------------
-!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
-!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
-!> @brief material subroutine incorporating kinematics resulting from opening of slip planes
-!> @details to be done
-!--------------------------------------------------------------------------------------------------
-submodule(phase:eigen) slipplaneopening
-
-  integer, dimension(:), allocatable :: damage_isoductile_instance
-
-  type :: tParameters                                                                               !< container type for internal constitutive parameters
-    integer :: &
-      sum_N_sl                                                                                      !< total number of cleavage planes
-    real(pReal) :: &
-      dot_o, &                                                                                      !< opening rate of cleavage planes
-      q                                                                                             !< damage rate sensitivity
-    real(pReal), dimension(:),   allocatable :: &
-      g_crit
-    real(pReal), dimension(:,:,:), allocatable     :: &
-      P_d, &
-      P_t, &
-      P_n
-  end type tParameters
-
-  type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters (len Ninstances)
-
-
-contains
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief module initialization
-!> @details reads in material parameters, allocates arrays, and does sanity checks
-!--------------------------------------------------------------------------------------------------
-module function damage_isoductile_init() result(myKinematics)
-
-  logical, dimension(:), allocatable :: myKinematics
-
-  integer :: p,i
-  character(len=pStringLen) :: extmsg = ''
-  integer,     dimension(:),   allocatable :: N_sl
-  real(pReal), dimension(:,:), allocatable :: d,n,t
-  class(tNode), pointer :: &
-    phases, &
-    phase, &
-    mech, &
-    pl, &
-    kinematics, &
-    kinematic_type
-
-
-  myKinematics = kinematics_active2('isoductile')
-  if(count(myKinematics) == 0) return
-  print'(/,a)', ' <<<+-  phase:mechanical:eigen:slipplaneopening init  -+>>>'
-  print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT)
-
-
-  phases => config_material%get('phase')
-  allocate(param(phases%length))
-
-  do p = 1, phases%length
-    if(myKinematics(p)) then
-      phase => phases%get(p)
-      mech  => phase%get('mechanical')
-      pl    => mech%get('plastic')
-
-      kinematics => phase%get('damage')
-
-      associate(prm  => param(p))
-        kinematic_type => kinematics%get(1)
-
-        prm%dot_o    = kinematic_type%get_asFloat('dot_o')
-        prm%q        = kinematic_type%get_asFloat('q')
-        N_sl         = pl%get_as1dInt('N_sl')
-        prm%sum_N_sl = sum(abs(N_sl))
-
-        d = lattice_slip_direction (N_sl,phase%get_asString('lattice'),&
-                                    phase%get_asFloat('c/a',defaultVal=0.0_pReal))
-        t = lattice_slip_transverse(N_sl,phase%get_asString('lattice'),&
-                                    phase%get_asFloat('c/a',defaultVal=0.0_pReal))
-        n = lattice_slip_normal    (N_sl,phase%get_asString('lattice'),&
-                                    phase%get_asFloat('c/a',defaultVal=0.0_pReal))
-        allocate(prm%P_d(3,3,size(d,2)),prm%P_t(3,3,size(t,2)),prm%P_n(3,3,size(n,2)))
-
-        do i=1, size(n,2)
-          prm%P_d(1:3,1:3,i) = math_outer(d(1:3,i), n(1:3,i))
-          prm%P_t(1:3,1:3,i) = math_outer(t(1:3,i), n(1:3,i))
-          prm%P_n(1:3,1:3,i) = math_outer(n(1:3,i), n(1:3,i))
-        enddo
-
-        prm%g_crit = kinematic_type%get_as1dFloat('g_crit',requiredSize=size(N_sl))
-
-        ! expand: family => system
-        prm%g_crit = math_expand(prm%g_crit,N_sl)
-
-        ! sanity checks
-        if (prm%q          <= 0.0_pReal)  extmsg = trim(extmsg)//' anisoDuctile_n'
-        if (prm%dot_o      <= 0.0_pReal)  extmsg = trim(extmsg)//' anisoDuctile_sdot0'
-        if (any(prm%g_crit <  0.0_pReal)) extmsg = trim(extmsg)//' anisoDuctile_critLoad'
-
-!--------------------------------------------------------------------------------------------------
-!  exit if any parameter is out of range
-        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(slipplane_opening)')
-
-      end associate
-    endif
-  enddo
-
-
-end function damage_isoductile_init
-
-
-!--------------------------------------------------------------------------------------------------
-!> @brief  contains the constitutive equation for calculating the velocity gradient
-!--------------------------------------------------------------------------------------------------
-module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
-
-  integer, intent(in) :: &
-    ph, me
-  real(pReal),   intent(in),  dimension(3,3) :: &
-    S
-  real(pReal),   intent(out), dimension(3,3) :: &
-    Ld                                                                                              !< damage velocity gradient
-  real(pReal),   intent(out), dimension(3,3,3,3) :: &
-    dLd_dTstar                                                                                      !< derivative of Ld with respect to Tstar (4th-order tensor)
-
-  integer :: &
-    i, k, l, m, n
-  real(pReal) :: &
-    traction_d, traction_t, traction_n, traction_crit, &
-    udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
-
-
-  associate(prm => param(ph))
-  Ld = 0.0_pReal
-  dLd_dTstar = 0.0_pReal
-  do i = 1, prm%sum_N_sl
-
-    traction_d = math_tensordot(S,prm%P_d(1:3,1:3,i))
-    traction_t = math_tensordot(S,prm%P_t(1:3,1:3,i))
-    traction_n = math_tensordot(S,prm%P_n(1:3,1:3,i))
-
-    traction_crit = prm%g_crit(i)* damage_phi(ph,me)
-
-    udotd = sign(1.0_pReal,traction_d)* prm%dot_o* (  abs(traction_d)/traction_crit &
-                                                    - abs(traction_d)/prm%g_crit(i))**prm%q
-    udott = sign(1.0_pReal,traction_t)* prm%dot_o* (  abs(traction_t)/traction_crit &
-                                                    - abs(traction_t)/prm%g_crit(i))**prm%q
-    udotn = prm%dot_o* (  max(0.0_pReal,traction_n)/traction_crit &
-                        - max(0.0_pReal,traction_n)/prm%g_crit(i))**prm%q
-
-    if (dNeq0(traction_d)) then
-      dudotd_dt = udotd*prm%q/traction_d
-    else
-      dudotd_dt = 0.0_pReal
-    endif
-    if (dNeq0(traction_t)) then
-      dudott_dt = udott*prm%q/traction_t
-    else
-      dudott_dt = 0.0_pReal
-    endif
-    if (dNeq0(traction_n)) then
-      dudotn_dt = udotn*prm%q/traction_n
-    else
-      dudotn_dt = 0.0_pReal
-    endif
-
-    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      dLd_dTstar(k,l,m,n) = dLd_dTstar(k,l,m,n) &
-                          + dudotd_dt*prm%P_d(k,l,i)*prm%P_d(m,n,i) &
-                          + dudott_dt*prm%P_t(k,l,i)*prm%P_t(m,n,i) &
-                          + dudotn_dt*prm%P_n(k,l,i)*prm%P_n(m,n,i)
-
-    Ld = Ld &
-       + udotd*prm%P_d(1:3,1:3,i) &
-       + udott*prm%P_t(1:3,1:3,i) &
-       + udotn*prm%P_n(1:3,1:3,i)
-  enddo
-
-  end associate
-
-end subroutine damage_isoductile_LiAndItsTangent
-
-end submodule slipplaneopening

src/phase_mechanical_plastic_dislotungsten.f90

@@ -249,8 +249,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
     dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
-    ! global alias
-    plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
 
     allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers),         source=0.0_pReal)
     allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers),  source=0.0_pReal)

src/phase_mechanical_plastic_dislotwin.f90

@@ -438,8 +438,6 @@ module function plastic_dislotwin_init() result(myPlasticity)
     dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
-    ! global alias
-    plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_tw

src/phase_mechanical_plastic_isotropic.f90

@@ -139,8 +139,6 @@ module function plastic_isotropic_init() result(myPlasticity)
     dot%gamma  => plasticState(ph)%dotState(2,:)
     plasticState(ph)%atol(2) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
     if (plasticState(ph)%atol(2) < 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma'
-    ! global alias
-    plasticState(ph)%slipRate => plasticState(ph)%dotState(2:2,:)
 
     end associate
 

src/phase_mechanical_plastic_kinehardening.f90

@@ -194,8 +194,6 @@ module function plastic_kinehardening_init() result(myPlasticity)
     dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
-    ! global alias
-    plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
 
     o = plasticState(ph)%offsetDeltaState
     startIndex = endIndex + 1

src/phase_mechanical_plastic_nonlocal.f90

@@ -491,8 +491,6 @@ module function plastic_nonlocal_init() result(myPlasticity)
     plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl )  = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-6_pReal)
     if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
       extmsg = trim(extmsg)//' atol_gamma'
-    ! global alias
-    plasticState(ph)%slipRate => plasticState(ph)%dotState   (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nmembers)
 
     stt%rho_forest => plasticState(ph)%state                 (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nmembers)
     stt%v          => plasticState(ph)%state                 (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nmembers)

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -254,8 +254,6 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
-    ! global alias
-    plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_tw

src/prec.f90

@@ -51,10 +51,7 @@ module prec
   end type
 
   type, extends(tState) :: tPlasticState
-    logical :: &
+    logical :: nonlocal = .false.
-      nonlocal = .false.
-    real(pReal), pointer, dimension(:,:) :: &
-      slipRate                                                                                      !< slip rate
   end type
 
   type :: tSourceState