diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index f2464545b..281a643d9 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -210,7 +210,7 @@ subroutine CPFEM_init m = 0_pInt readInstances: do ph = 1_pInt, size(phase_plasticity) do k = 1_pInt, plasticState(ph)%sizeState - do l = 1, size(plasticState(ph)%state(1,:)) + do l = 1, size(plasticState(ph)%state0(1,:)) m = m+1_pInt read(777,rec=m) plasticState(ph)%state0(k,l) enddo; enddo @@ -224,6 +224,7 @@ subroutine CPFEM_init m = m+1_pInt read(777,rec=m) homogenization_state0(j,k)%p(l) enddo + enddo; enddo close (777) #if defined(Marc4DAMASK) || defined(Abaqus) @@ -403,7 +404,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) forall ( i = 1:size(plasticState)) plasticState(i)%state0= plasticState(i)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array forall ( i = 1:size(damageState)) damageState(i)%state0 = damageState(i)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array forall ( i = 1:size(thermalState)) thermalState(i)%state0= thermalState(i)%state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then write(6,'(a)') '<< CPFEM >> aging states' if (debug_e <= mesh_NcpElems .and. debug_i <= mesh_maxNips) then write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)),/)') & @@ -424,7 +425,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) ! * dump the last converged values of each essential variable to a binary file if (restartWrite) then - if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & + if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) & write(6,'(a)') '<< CPFEM >> writing state variables of last converged step to binary files' call IO_write_jobRealFile(777,'recordedPhase',size(material_phase)) @@ -455,7 +456,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) m = 0_pInt writeInstances: do ph = 1_pInt, size(phase_plasticity) do k = 1_pInt, plasticState(ph)%sizeState - do l = 1, size(plasticState(ph)%state(1,:)) + do l = 1, size(plasticState(ph)%state0(1,:)) m = m+1_pInt write(777,rec=m) plasticState(ph)%state0(k,l) enddo; enddo @@ -559,14 +560,12 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) endif !* map stress and stiffness (or return odd values if terminally ill) - +#if defined(Marc4DAMASK) || defined(Abaqus) if ( terminallyIll ) then -#if defined(Marc4DAMASK) || defined(Abaqus) call random_number(rnd) if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal CPFEM_cs(1:6,ip,elCP) = rnd * CPFEM_odd_stress CPFEM_dcsde(1:6,1:6,ip,elCP) = CPFEM_odd_jacobian * math_identity2nd(6) -#endif else if (microstructure_elemhomo(mesh_element(4,elCP)) .and. ip > 1_pInt) then ! me homogenous? --> copy from first ip materialpoint_P(1:3,1:3,ip,elCP) = materialpoint_P(1:3,1:3,1,elCP) @@ -575,7 +574,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) materialpoint_results(1:materialpoint_sizeResults,ip,elCP) = & materialpoint_results(1:materialpoint_sizeResults,1,elCP) endif -#if defined(Marc4DAMASK) || defined(Abaqus) ! translate from P to CS Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), math_transpose33(materialpoint_F(1:3,1:3,ip,elCP))) J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) @@ -595,8 +593,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) forall(i=1:3, j=1:3,k=1:3,l=1:3) & H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k)) CPFEM_dcsde(1:6,1:6,ip,elCP) = math_Mandel3333to66(J_inverse * H_sym) -#endif endif +#endif endif #if defined(Marc4DAMASK) || defined(Abaqus) @@ -637,7 +635,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature, dt, elFE, ip) #if defined(Marc4DAMASK) || defined(Abaqus) !*** warn if stiffness close to zero if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip) - !*** copy to output if required (FEM solver) + !*** copy to output if using commercial FEM solver cauchyStress = CPFEM_cs (1:6, ip,elCP) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP) #endif diff --git a/code/DAMASK_spectral_utilities.f90 b/code/DAMASK_spectral_utilities.f90 index 5035e87c5..8708e18c6 100644 --- a/code/DAMASK_spectral_utilities.f90 +++ b/code/DAMASK_spectral_utilities.f90 @@ -827,14 +827,13 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& integer(pInt) :: & calcMode, & !< CPFEM mode for calculation - collectMode, & !< CPFEM mode for collection j,k real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet write(6,'(/,a)') ' ... evaluating constitutive response ......................................' calcMode = CPFEM_CALCRESULTS - collectMode = CPFEM_COLLECT + if (forwardData) then ! aging results calcMode = ior(calcMode, CPFEM_AGERESULTS) materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid)]) @@ -843,7 +842,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,& calcMode = iand(calcMode, not(CPFEM_AGERESULTS)) endif - call CPFEM_general(collectMode,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), & ! collect mode handles Jacobian backup / restoration + call CPFEM_general(CPFEM_COLLECT,F_lastInc(1:3,1:3,1,1,1),F(1:3,1:3,1,1,1), & temperature,timeinc,1_pInt,1_pInt) crystallite_temperature = temperature diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index 27ef4d69d..d8c3f93d8 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -608,6 +608,7 @@ subroutine constitutive_phenopowerlaw_stateInit(ph,instance) i real(pReal), dimension(plasticState(ph)%sizeState) :: & tempState + tempState = 0.0_pReal do i = 1_pInt,lattice_maxNslipFamily tempState(1+sum(constitutive_phenopowerlaw_Nslip(1:i-1,instance)) : & @@ -624,8 +625,8 @@ subroutine constitutive_phenopowerlaw_stateInit(ph,instance) enddo plasticState(ph)%state0(:,:) = spread(tempState,2,size(plasticState(ph)%state0(1,:))) -end subroutine constitutive_phenopowerlaw_stateInit +end subroutine constitutive_phenopowerlaw_stateInit !-------------------------------------------------------------------------------------------------- diff --git a/code/crystallite.f90 b/code/crystallite.f90 index f7f25987f..7be9504c1 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -3052,9 +3052,9 @@ subroutine crystallite_integrateStateFPI() if (crystallite_todo(g,i,e)) then p = mappingConstitutive(2,g,i,e) c = mappingConstitutive(1,g,i,e) - if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or.& + if ( any(plasticState(p)%dotState(:,c) /= plasticState(p)%dotState(:,c)) .or.& any(damageState(p)%dotState(:,c) /= damageState(p)%dotState(:,c)) .or.& - any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c))) then !NaN occured in dotState + any(thermalState(p)%dotState(:,c) /= thermalState(p)%dotState(:,c))) then !NaN occured in dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... !$OMP CRITICAL (checkTodo)