added dislokmc (dislotwin modification)

code/DAMASK_abaqus_exp.f

@@ -95,6 +95,7 @@ end module DAMASK_interface
 #include "constitutive_phenopowerlaw.f90"
 #include "constitutive_titanmod.f90"
 #include "constitutive_dislotwin.f90"
+#include "constitutive_dislokmc.f90"
 #include "constitutive_nonlocal.f90"
 #include "constitutive.f90"
 #include "crystallite.f90"

code/DAMASK_abaqus_std.f

@@ -95,6 +95,7 @@ end module DAMASK_interface
 #include "constitutive_phenopowerlaw.f90"
 #include "constitutive_titanmod.f90"
 #include "constitutive_dislotwin.f90"
+#include "constitutive_dislokmc.f90"
 #include "constitutive_nonlocal.f90"
 #include "constitutive.f90"
 #include "crystallite.f90"

code/DAMASK_marc.f90

@@ -124,6 +124,7 @@ end module DAMASK_interface
 #include "constitutive_phenopowerlaw.f90"
 #include "constitutive_titanmod.f90"
 #include "constitutive_dislotwin.f90"
+#include "constitutive_dislokmc.f90"
 #include "constitutive_nonlocal.f90"
 #include "constitutive.f90"
 #include "crystallite.f90"

code/Makefile

@@ -349,7 +349,7 @@ COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHEC
 COMPILED_FILES = prec.o DAMASK_spectral_interface.o IO.o libs.o numerics.o debug.o math.o \
                  FEsolving.o mesh.o material.o lattice.o \
                  damage_none.o damage_local.o damage_gradient.o thermal_none.o thermal_conduction.o \
-                 constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o \
+                 constitutive_dislotwin.o constitutive_dislokmc.o constitutive_j2.o constitutive_phenopowerlaw.o \
                  constitutive_titanmod.o constitutive_nonlocal.o constitutive_none.o \
                  constitutive_damage.o constitutive_thermal.o constitutive.o crystallite.o \
                  homogenization_RGC.o homogenization_isostrain.o homogenization_none.o homogenization.o CPFEM.o \
@@ -411,6 +411,7 @@ constitutive.o:                constitutive.f90 \
                                constitutive_nonlocal.o \
                                constitutive_titanmod.o \
                                constitutive_dislotwin.o \
+                               constitutive_dislokmc.o \
                                constitutive_phenopowerlaw.o \
                                constitutive_j2.o \
                                constitutive_none.o
@@ -421,6 +422,9 @@ constitutive_nonlocal.o:       constitutive_nonlocal.f90 \
 constitutive_titanmod.o:       constitutive_titanmod.f90 \
                                lattice.o
 
+constitutive_dislokmc.o:       constitutive_dislokmc.f90 \
+                               lattice.o
+
 constitutive_dislotwin.o:      constitutive_dislotwin.f90 \
                                lattice.o
 

code/constitutive.f90

@@ -81,6 +81,7 @@ subroutine constitutive_init
    PLASTICITY_J2_ID, &
    PLASTICITY_PHENOPOWERLAW_ID, &
    PLASTICITY_DISLOTWIN_ID, &
+   PLASTICITY_DISLOKMC_ID, &
    PLASTICITY_TITANMOD_ID, &
    PLASTICITY_NONLOCAL_ID ,&
    ELASTICITY_HOOKE_label, &
@@ -88,15 +89,18 @@ subroutine constitutive_init
    PLASTICITY_J2_label, &
    PLASTICITY_PHENOPOWERLAW_label, &
    PLASTICITY_DISLOTWIN_label, &
+   PLASTICITY_DISLOKMC_label, &
    PLASTICITY_TITANMOD_label, &
+   PLASTICITY_NONLOCAL_label, &
    plasticState, &
-   mappingConstitutive, &
+   mappingConstitutive
+ 
 
-   PLASTICITY_NONLOCAL_label
  use constitutive_none
  use constitutive_j2
  use constitutive_phenopowerlaw
  use constitutive_dislotwin
+ use constitutive_dislokmc
  use constitutive_titanmod
  use constitutive_nonlocal
  implicit none
@@ -120,6 +124,7 @@ subroutine constitutive_init
  if (any(phase_plasticity == PLASTICITY_J2_ID))            call constitutive_j2_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call constitutive_phenopowerlaw_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID))     call constitutive_dislotwin_init(FILEUNIT)
+ if (any(phase_plasticity == PLASTICITY_DISLOKMC_ID))      call constitutive_dislokmc_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_TITANMOD_ID))      call constitutive_titanmod_init(FILEUNIT)
  if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then
   call constitutive_nonlocal_init(FILEUNIT)
@@ -155,6 +160,10 @@ subroutine constitutive_init
        outputName = PLASTICITY_DISLOTWIN_label
        thisOutput => constitutive_dislotwin_output
        thisSize   => constitutive_dislotwin_sizePostResult
+     case (PLASTICITY_DISLOKMC_ID)
+       outputName = PLASTICITY_DISLOKMC_label
+       thisOutput => constitutive_dislokmc_output
+       thisSize   => constitutive_dislokmc_sizePostResult
      case (PLASTICITY_TITANMOD_ID)
        outputName = PLASTICITY_TITANMOD_label
        thisOutput => constitutive_titanmod_output
@@ -244,13 +253,17 @@ function constitutive_homogenizedC(ipc,ip,el)
    phase_plasticity, &
    material_phase, &
    PLASTICITY_TITANMOD_ID, &
+   PLASTICITY_DISLOTWIN_ID, &
+   PLASTICITY_DISLOKMC_ID, &
    plasticState,&
-   mappingConstitutive, &
+   mappingConstitutive
-   PLASTICITY_DISLOTWIN_ID
+
  use constitutive_titanmod, only: &
    constitutive_titanmod_homogenizedC
  use constitutive_dislotwin, only: &
    constitutive_dislotwin_homogenizedC
+ use constitutive_dislokmc, only: &
+   constitutive_dislokmc_homogenizedC
  use lattice, only: &
    lattice_C66
 
@@ -265,6 +278,8 @@ function constitutive_homogenizedC(ipc,ip,el)
 
    case (PLASTICITY_DISLOTWIN_ID)
      constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(ipc,ip,el)
+   case (PLASTICITY_DISLOKMC_ID)
+     constitutive_homogenizedC = constitutive_dislokmc_homogenizedC(ipc,ip,el) 
    case (PLASTICITY_TITANMOD_ID)
      constitutive_homogenizedC = constitutive_titanmod_homogenizedC (ipc,ip,el)
    case default
@@ -285,16 +300,20 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el)
    phase_plasticity, &
    material_phase, &
    PLASTICITY_DISLOTWIN_ID, &
-   plasticState, &
+   PLASTICITY_DISLOKMC_ID, &
-   mappingConstitutive, &
    PLASTICITY_TITANMOD_ID, &
-   PLASTICITY_NONLOCAL_ID
+   PLASTICITY_NONLOCAL_ID, &
+   plasticState, &
+   mappingConstitutive
+
  use constitutive_titanmod, only: &
    constitutive_titanmod_microstructure
  use constitutive_nonlocal, only: &
    constitutive_nonlocal_microstructure
  use constitutive_dislotwin, only: &
    constitutive_dislotwin_microstructure
+ use constitutive_dislokmc, only: &
+   constitutive_dislokmc_microstructure
 
  implicit none
  integer(pInt), intent(in) :: &
@@ -311,6 +330,8 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el)
        
    case (PLASTICITY_DISLOTWIN_ID)
      call constitutive_dislotwin_microstructure(temperature,ipc,ip,el)
+   case (PLASTICITY_DISLOKMC_ID)
+     call constitutive_dislokmc_microstructure(temperature,ipc,ip,el)
    case (PLASTICITY_TITANMOD_ID)
      call constitutive_titanmod_microstructure (temperature,ipc,ip,el)
    case (PLASTICITY_NONLOCAL_ID)
@@ -338,6 +359,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
    PLASTICITY_J2_ID, &
    PLASTICITY_PHENOPOWERLAW_ID, &
    PLASTICITY_DISLOTWIN_ID, &
+   PLASTICITY_DISLOKMC_ID, &
    PLASTICITY_TITANMOD_ID, &
    PLASTICITY_NONLOCAL_ID
  use constitutive_j2, only: &
@@ -346,6 +368,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
    constitutive_phenopowerlaw_LpAndItsTangent
  use constitutive_dislotwin, only: &
    constitutive_dislotwin_LpAndItsTangent
+ use constitutive_dislokmc, only: &
+   constitutive_dislokmc_LpAndItsTangent
  use constitutive_titanmod, only: &
    constitutive_titanmod_LpAndItsTangent
  use constitutive_nonlocal, only: &
@@ -378,6 +402,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
      call constitutive_nonlocal_LpAndItsTangent     (Lp,dLp_dTstar,Tstar_v,temperature,    ip,el)
    case (PLASTICITY_DISLOTWIN_ID)
      call constitutive_dislotwin_LpAndItsTangent    (Lp,dLp_dTstar,Tstar_v,temperature,ipc,ip,el)
+   case (PLASTICITY_DISLOKMC_ID)
+     call constitutive_dislokmc_LpAndItsTangent    (Lp,dLp_dTstar,Tstar_v,temperature,ipc,ip,el)
    case (PLASTICITY_TITANMOD_ID)
      call constitutive_titanmod_LpAndItsTangent     (Lp,dLp_dTstar,Tstar_v,temperature,ipc,ip,el)
 
@@ -526,6 +552,7 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature,
    PLASTICITY_J2_ID, &
    PLASTICITY_PHENOPOWERLAW_ID, &
    PLASTICITY_DISLOTWIN_ID, &
+   PLASTICITY_DISLOKMC_ID, &
    PLASTICITY_TITANMOD_ID, &
    PLASTICITY_NONLOCAL_ID
  use constitutive_j2, only:  &
@@ -534,6 +561,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature,
    constitutive_phenopowerlaw_dotState
  use constitutive_dislotwin, only: &
    constitutive_dislotwin_dotState
+ use constitutive_dislokmc, only: &
+   constitutive_dislokmc_dotState
  use constitutive_titanmod, only: &
    constitutive_titanmod_dotState
  use constitutive_nonlocal, only: &
@@ -569,6 +598,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature,
      call constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
    case (PLASTICITY_DISLOTWIN_ID)
      call constitutive_dislotwin_dotState    (Tstar_v,Temperature,ipc,ip,el)
+   case (PLASTICITY_DISLOKMC_ID)
+     call constitutive_dislokmc_dotState    (Tstar_v,Temperature,ipc,ip,el)
    case (PLASTICITY_TITANMOD_ID)
      call constitutive_titanmod_dotState     (Tstar_v,Temperature,ipc,ip,el)
    case (PLASTICITY_NONLOCAL_ID)
@@ -669,6 +700,7 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
    PLASTICITY_J2_ID, &
    PLASTICITY_PHENOPOWERLAW_ID, &
    PLASTICITY_DISLOTWIN_ID, &
+   PLASTICITY_DISLOKMC_ID, &
    PLASTICITY_TITANMOD_ID, &
    PLASTICITY_NONLOCAL_ID
  use constitutive_j2, only: &
@@ -680,6 +712,8 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
    constitutive_phenopowerlaw_postResults
  use constitutive_dislotwin, only: &
    constitutive_dislotwin_postResults
+ use constitutive_dislokmc, only: &
+   constitutive_dislokmc_postResults
  use constitutive_titanmod, only: &
    constitutive_titanmod_postResults
  use constitutive_nonlocal, only: &
@@ -709,6 +743,8 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
      constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
    case (PLASTICITY_DISLOTWIN_ID)
      constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
+   case (PLASTICITY_DISLOKMC_ID)
+     constitutive_postResults = constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
    case (PLASTICITY_NONLOCAL_ID)
      constitutive_postResults = constitutive_nonlocal_postResults (Tstar_v,FeArray,        ip,el)
  end select

code/constitutive_dislokmc.f90

@@ -232,10 +232,6 @@ subroutine constitutive_dislokmc_init(fileUnit)
  if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
    write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
 
- allocate(constitutive_dislotkmc_sizeDotState(maxNinstance),                       source=0_pInt)
- allocate(constitutive_dislotkmc_sizeState(maxNinstance),                          source=0_pInt)
- allocate(constitutive_dislokmc_sizeDotState(maxNinstance),                        source=0_pInt)
- allocate(constitutive_dislokmc_sizeState(maxNinstance),                           source=0_pInt)
  allocate(constitutive_dislokmc_sizePostResults(maxNinstance),                     source=0_pInt)
  allocate(constitutive_dislokmc_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
  allocate(constitutive_dislokmc_output(maxval(phase_Noutput),maxNinstance))
@@ -1316,7 +1312,7 @@ subroutine constitutive_dislokmc_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperatu
  
         !* Shear rates due to slip                                                                                                                                              
         vel_slip(j) = exp(-BoltzmannRatio*(1-StressRatio_p) ** constitutive_dislokmc_qPerSlipFamily(f,instance)) &
-                     * (1-constitutive_dislokmc_sPerSlipFamily(f,instance) &
+                     * (1.0_pReal-constitutive_dislokmc_sPerSlipFamily(f,instance) &
                      * exp(-BoltzmannRatio*(1-StressRatio_p) ** constitutive_dislokmc_qPerSlipFamily(f,instance)))
                        
         gdot_slip(j) = (1.0_pReal - sumf) * DotGamma0 &
@@ -1552,10 +1548,10 @@ subroutine constitutive_dislokmc_dotState(Tstar_v,Temperature,ipc,ip,el)
       (constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
          **(constitutive_dislokmc_pPerSlipFamily(f,instance)-1.0_pReal)
 
-        StressRatio_u = ((abs(tau_slip(j))-state%p(6*ns+4*nt+j, of))/&
+        StressRatio_u = ((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/&
          (constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
        **constitutive_dislokmc_uPerSlipFamily(f,instance)
-        StressRatio_uminus1 = ((abs(tau_slip(j))-state%p(6*ns+4*nt+j,of))/&
+        StressRatio_uminus1 = ((abs(tau_slip(j))-plasticState(ph)%state(6*ns+4*nt+j,of))/&
    (constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
            **(constitutive_dislokmc_uPerSlipFamily(f,instance)-1.0_pReal)
 
@@ -1803,7 +1799,7 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
                 constitutive_dislokmc_postResults(c+j) = &
                   DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
                                constitutive_dislokmc_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau) & 
-                               * (1-constitutive_dislokmc_sPerSlipFamily(f,instance) &
+                               * (1.0_pReal-constitutive_dislokmc_sPerSlipFamily(f,instance) &
                                * exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislokmc_qPerSlipFamily(f,instance)))
               else
                 constitutive_dislokmc_postResults(c+j) = 0.0_pReal
@@ -1844,8 +1840,6 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
                 (16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph))))
               constitutive_dislokmc_postResults(c+j)=min(constitutive_dislokmc_postResults(c+j),&
                                                             plasticState(ph)%state(5*ns+3*nt+j, of))
- !            constitutive_dislokmc_postResults(c+j)=max(constitutive_dislokmc_postResults(c+j),&
- !                                                            plasticState(ph)%state(4*ns+2*nt+j, of))
         enddo; enddo
         c = c + ns
        case (resolved_stress_shearband_ID)
@@ -1904,10 +1898,10 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
                                         constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
                    **(constitutive_dislokmc_pPerSlipFamily(f,instance)-1.0_pReal)
 
-                 StressRatio_u = ((abs(tau)-state%p(6*ns+4*nt+j, of))/&
+                 StressRatio_u = ((abs(tau)-plasticState(ph)%state(6*ns+4*nt+j, of))/&
      (constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
             **constitutive_dislokmc_uPerSlipFamily(f,instance)
-                  StressRatio_uminus1 = ((abs(tau)-state%p(6*ns+4*nt+j, of))/&
+                  StressRatio_uminus1 = ((abs(tau)-plasticState(ph)%state(6*ns+4*nt+j, of))/&
    (constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
                **(constitutive_dislokmc_uPerSlipFamily(f,instance)-1.0_pReal)
 
@@ -1921,7 +1915,7 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
                !* Shear rates due to slip                                                                                                                                                      
                  gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**&
                                 constitutive_dislokmc_qPerSlipFamily(f,instance))*sign(1.0_pReal,tau) &
-                     * (1_pInt-constitutive_dislokmc_sPerSlipFamily(f,instance) &
+                     * (1.0_pReal-constitutive_dislokmc_sPerSlipFamily(f,instance) &
                 * exp(-BoltzmannRatio*(1_pInt-StressRatio_p) ** constitutive_dislokmc_qPerSlipFamily(f,instance)) ) &
                 * StressRatio_u
                else
@@ -2010,10 +2004,10 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
                                      (constitutive_dislokmc_SolidSolutionStrength(instance)+&
                                       constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
                        **(constitutive_dislokmc_pPerSlipFamily(f,instance)-1.0_pReal)
-               StressRatio_u = ((abs(tau)-state%p(6*ns+4*nt+j))/&
+               StressRatio_u = ((abs(tau)-plasticState(ph)%state(6*ns+4*nt+j,of))/&
            (constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
                              **constitutive_dislokmc_uPerSlipFamily(f,instance)
-               StressRatio_uminus1 = ((abs(tau)-state%p(6*ns+4*nt+j))/&
+               StressRatio_uminus1 = ((abs(tau)-plasticState(ph)%state(6*ns+4*nt+j,of))/&
        (constitutive_dislokmc_SolidSolutionStrength(instance)+constitutive_dislokmc_tau_peierlsPerSlipFamily(f,instance)))&
                          **(constitutive_dislokmc_uPerSlipFamily(f,instance)-1.0_pReal)
 
@@ -2026,7 +2020,7 @@ function constitutive_dislokmc_postResults(Tstar_v,Temperature,ipc,ip,el)
 
              !* Shear rates due to slip                                                                                                                      
                vel_slip(j) = exp(-BoltzmannRatio*(1-StressRatio_p) ** constitutive_dislokmc_qPerSlipFamily(f,instance)) &
-                     * (1-constitutive_dislokmc_sPerSlipFamily(f,instance) &
+                     * (1.0_pReal-constitutive_dislokmc_sPerSlipFamily(f,instance) &
                      * exp(-BoltzmannRatio*(1-StressRatio_p) ** constitutive_dislokmc_qPerSlipFamily(f,instance)))
 
                gdot_slip(j) = (1.0_pReal - sumf) * DotGamma0 &

code/material.f90

@@ -23,6 +23,7 @@ module material
    PLASTICITY_J2_label            = 'j2', &
    PLASTICITY_PHENOPOWERLAW_label = 'phenopowerlaw', &
    PLASTICITY_DISLOTWIN_label     = 'dislotwin', &
+   PLASTICITY_DISLOKMC_label      = 'dislokmc', &
    PLASTICITY_TITANMOD_label      = 'titanmod', &
    PLASTICITY_NONLOCAL_label      = 'nonlocal', &
    DAMAGE_NONE_label              = 'none', &
@@ -46,6 +47,7 @@ module material
                  PLASTICITY_J2_ID, &
                  PLASTICITY_phenopowerlaw_ID, &
                  PLASTICITY_dislotwin_ID, &
+                 PLASTICITY_dislokmc_ID, &
                  PLASTICITY_titanmod_ID, &
                  PLASTICITY_nonlocal_ID
  end enum
@@ -186,6 +188,7 @@ module material
    PLASTICITY_J2_ID, &
    PLASTICITY_phenopowerlaw_ID, &
    PLASTICITY_dislotwin_ID, &
+   PLASTICITY_dislokmc_ID, &
    PLASTICITY_titanmod_ID, &
    PLASTICITY_nonlocal_ID, &
    DAMAGE_none_ID, &
@@ -670,6 +673,8 @@ subroutine material_parsePhase(fileUnit,myPart)
              phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID
            case (PLASTICITY_DISLOTWIN_label)
              phase_plasticity(section) = PLASTICITY_DISLOTWIN_ID
+           case (PLASTICITY_DISLOKMC_label)
+             phase_plasticity(section) = PLASTICITY_DISLOKMC_ID
            case (PLASTICITY_TITANMOD_label)
              phase_plasticity(section) = PLASTICITY_TITANMOD_ID
            case (PLASTICITY_NONLOCAL_label)