diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index aea276d4f..8eebf7d7d 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -300,7 +300,11 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el #endif material_phase use constitutive, only: & - constitutive_state0,constitutive_state + constitutive_state0, & +#ifdef NEWSTATE + mappingConstitutive,& +#endif + constitutive_state use crystallite, only: & crystallite_partionedF,& crystallite_F0, & @@ -389,17 +393,29 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress + + + forall ( i = 1:homogenization_maxNgrains, & j = 1:mesh_maxNips, & k = 1:mesh_NcpElems ) & constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites + #ifdef NEWSTATE -!(:) needed? A component cannot be an array if the encompassing structure is an array -!plasticState(:)%state0(:,:)=plasticState(:)%state(:,:) - forall ( i = 1:size(plasticState)) & - plasticState(i)%state0= plasticState(i)%state ! microstructure of crystallites + 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 + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + 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)),/)') '<< CPFEM >> aged state of elFE ip grain',& + debug_e, debug_i, 1, & + plasticState(mappingConstitutive(2,1,debug_i,debug_e))%state(:,mappingConstitutive(1,1,debug_i,debug_e)) + endif + !$OMP END CRITICAL (write2out) + endif #endif + if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then !$OMP CRITICAL (write2out) write(6,'(a)') '<< CPFEM >> aging states' @@ -409,6 +425,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el endif !$OMP END CRITICAL (write2out) endif + !$OMP PARALLEL DO do k = 1,mesh_NcpElems do j = 1,mesh_maxNips diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 0beddd718..1d477a6d5 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -147,8 +147,8 @@ subroutine constitutive_init character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready logical :: knownPlasticity, nonlocalConstitutionPresent #if defined(HDF) || defined(NEWSTATE) - allocate(mappingConstitutive(homogenization_maxngrains,mesh_maxNips,mesh_ncpelems,2),source=0_pInt) - allocate(mappingCrystallite (homogenization_maxngrains,mesh_ncpelems,2),source=0_pInt) + allocate(mappingConstitutive(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems),source=0_pInt) + allocate(mappingCrystallite (2,homogenization_maxngrains,mesh_ncpelems),source=0_pInt) allocate(ConstitutivePosition(material_nphase),source=0_pInt) allocate(CrystallitePosition(material_nphase),source=0_pInt) #endif @@ -221,7 +221,6 @@ subroutine constitutive_init iMax = mesh_maxNips eMax = mesh_NcpElems - ! lumped into new state allocate(constitutive_state0(cMax,iMax,eMax)) @@ -258,8 +257,8 @@ subroutine constitutive_init phase = material_phase(g,i,e) instance = phase_plasticityInstance(phase) #if defined(HDF) || defined(NEWSTATE) - ConstitutivePosition(phase) = ConstitutivePosition(phase)+1_pInt - mappingConstitutive(g,i,e,1:2) = [ConstitutivePosition(phase),phase] + ConstitutivePosition(phase) = ConstitutivePosition(phase)+1_pInt ! not distinguishing between instances of same phase + mappingConstitutive(1:2,g,i,e) = [ConstitutivePosition(phase),phase] #endif select case(phase_plasticity(material_phase(g,i,e))) case (PLASTICITY_NONE_ID) diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 0f0ea086b..225ead1f7 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -316,7 +316,7 @@ subroutine constitutive_j2_init(fileUnit) allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%state (sizeState,NofMyPhase),source=constitutive_j2_tau0(instance)) allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%aTolState (sizeState,NofMyPhase),constitutive_j2_aTolResistance(instance)) + allocate(plasticState(phase)%aTolState (sizeState,NofMyPhase),source=constitutive_j2_aTolResistance(instance)) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase),source=0.0_pReal) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 04613076e..841d75e83 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -529,6 +529,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) FE_cellType use material, only: & homogenization_Ngrains, & +#ifdef NEWSTATE + plasticState, & +#endif homogenization_maxNgrains use constitutive, only: & constitutive_sizeState, & @@ -539,6 +542,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) constitutive_partionedState0, & constitutive_dotState, & constitutive_dotState_backup, & +#ifdef NEWSTATE + mappingConstitutive, & +#endif constitutive_TandItsTangent implicit none @@ -921,6 +927,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) +#endif crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress if (crystallite_syncSubFrac(i,e)) then ! if we just did a synchronization of states, then we wind forward without any further time integration crystallite_syncSubFracCompleted(i,e) = .true. @@ -965,6 +975,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) !$OMP FLUSH(crystallite_invFp) crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e)) +#endif crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e) ! ...2nd PK stress ! cant restore dotState here, since not yet calculated in first cutback after initialization crystallite_todo(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair) @@ -1180,6 +1194,13 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) constitutive_state(g,i,e)%p(1:constitutive_sizeState(g,i,e)) ! remember unperturbed, converged state, ... constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) = & constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) ! ... dotStates, ... +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) + plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) +#endif + F_backup(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ... and kinematics Fp_backup(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) InvFp_backup(1:3,1:3,g,i,e) = crystallite_invFp(1:3,1:3,g,i,e) @@ -1208,7 +1229,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) ! --- INITIALIZE UNPERTURBED STATE --- select case(numerics_integrator(numerics_integrationMode)) - case(1_pInt) ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state + case(1_pInt) ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains) @@ -1216,6 +1237,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) constitutive_state_backup(g,i,e)%p(1:constitutive_sizeState(g,i,e)) constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) = & constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e)) + plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) +#endif crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e) crystallite_invFp(1:3,1:3,g,i,e) = InvFp_backup(1:3,1:3,g,i,e) crystallite_Fe(1:3,1:3,g,i,e) = Fe_backup(1:3,1:3,g,i,e) @@ -1223,8 +1250,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) crystallite_Tstar_v(1:6,g,i,e) = Tstar_v_backup(1:6,g,i,e) endforall enddo - case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step - case(4_pInt,5_pInt) ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress + case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step + case(4_pInt,5_pInt) ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains) @@ -1232,6 +1259,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) constitutive_subState0(g,i,e)%p(1:constitutive_sizeState(g,i,e)) constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) = & constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e)) + plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) +#endif crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) crystallite_Fe(1:3,1:3,g,i,e) = crystallite_subFe0(1:3,1:3,g,i,e) crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) @@ -1318,6 +1351,13 @@ subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity) constitutive_state_backup(g,i,e)%p(1:constitutive_sizeState(g,i,e)) constitutive_dotState(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) = & constitutive_dotState_backup(g,i,e)%p(1:constitutive_sizeDotState(g,i,e)) +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e)) + plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) +#endif + crystallite_subF(1:3,1:3,g,i,e) = F_backup(1:3,1:3,g,i,e) crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e) crystallite_invFp(1:3,1:3,g,i,e) = InvFp_backup(1:3,1:3,g,i,e) @@ -1371,6 +1411,9 @@ subroutine crystallite_integrateStateRK4() mesh_maxNips use material, only: & homogenization_Ngrains, & +#ifdef NEWSTATE + plasticState, & +#endif homogenization_maxNgrains use constitutive, only: & constitutive_sizeDotState, & @@ -1381,6 +1424,9 @@ subroutine crystallite_integrateStateRK4() constitutive_collectDotState, & constitutive_deltaState, & constitutive_collectDeltaState, & +#ifdef NEWSTATE + mappingConstitutive, & +#endif constitutive_microstructure implicit none @@ -1413,6 +1459,9 @@ subroutine crystallite_integrateStateRK4() !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains constitutive_RK4dotState(g,i,e)%p = 0.0_pReal ! initialize Runge-Kutta dotState +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,g,i,e))%RK4dotState(:,mappingConstitutive(1,g,i,e)) = 0.0_pReal +#endif if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, & crystallite_Fp, crystallite_temperature(i,e), & @@ -1424,7 +1473,21 @@ subroutine crystallite_integrateStateRK4() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then +#ifdef NEWSTATE + if ( any(plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e)) & + /= plasticState(mappingConstitutive(2,g,i,e))%dotState(:,mappingConstitutive(1,g,i,e))) )then ! NaN occured in dotState + if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL (checkTodo) + crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped + !$OMP END CRITICAL (checkTodo) + else ! if broken local... + crystallite_todo(g,i,e) = .false. ! ... skip this one next time + endif + endif + +#endif if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) )then ! NaN occured in dotState + if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped @@ -3547,4 +3610,11 @@ function crystallite_postResults(ipc, ip, el) end function crystallite_postResults +!subroutine gradients +!use DAMASK_spectral_utilities + +!implicit none + +!end subroutine gradients + end module crystallite diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 4a33bb9d1..2bde49f46 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -373,6 +373,10 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness crystallite_partionedF0(1:3,1:3,g,i,e) = crystallite_F0(1:3,1:3,g,i,e) ! ...def grads crystallite_partionedTstar0_v(1:6,g,i,e) = crystallite_Tstar0_v(1:6,g,i,e) ! ...2nd PK stress +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = & + plasticState(mappingConstitutive(2,g,i,e))%state0(:,mappingConstitutive(1,g,i,e)) +#endif endforall forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e)) materialpoint_subF0(1:3,1:3,i,e) = materialpoint_F0(1:3,1:3,i,e) ! ...def grad @@ -384,15 +388,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), homogenization_sizeState(i,e) > 0_pInt) & homogenization_subState0(i,e)%p = homogenization_state0(i,e)%p ! ...internal homogenization state enddo -#ifdef NEWSTATE - do i = FEsolving_execIP(1,e), FEsolving_execIP(2,e) - do g = 1, myNgrains - plasticState(mappingConstitutive(g,i,e,1))%partionedState0(:,mappingConstitutive(g,i,e,2)) = & - plasticState(mappingConstitutive(g,i,e,1))%state0(:,mappingConstitutive(g,i,e,2)) - enddo - enddo -#endif - NiterationHomog = 0_pInt cutBackLooping: do while (.not. terminallyIll .and. &