From 18f60a94a9680878930f073fc27efad79f5268a5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 09:04:08 +0200 Subject: [PATCH 01/14] one variable suffices --- src/crystallite.f90 | 30 +++++++++++++----------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 8a5a89d62..ece4db020 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1014,15 +1014,13 @@ subroutine integrateStateFPI sizeDotState real(pReal) :: & zeta - real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & - residuum_plastic ! residuum for plastic state - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & - residuum_source ! residuum for source state + real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & + r ! state residuum logical :: & nonlocalBroken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1095,13 +1093,12 @@ subroutine integrateStateFPI plasticState(p)%previousDotState2(:,c)) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) - residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - - plasticState(p)%subState0(1:sizeDotState,c) & - - plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) + r(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & + - plasticState(p)%subState0(1:sizeDotState,c) & + - plasticState(p)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & - - residuum_plastic(1:sizeDotState) - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), & + - r(1:sizeDotState) + crystallite_converged(g,i,e) = converged(r(1:sizeDotState), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) @@ -1111,14 +1108,13 @@ subroutine integrateStateFPI sourceState(p)%p(s)%previousDotState2(:,c)) sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & + sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - zeta) - residuum_source(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) & - - sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - - sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) + r(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) & + - sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + - sourceState(p)%p(s)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) & - - residuum_source(1:sizeDotState) + - r(1:sizeDotState) crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState), & + crystallite_converged(g,i,e) .and. converged(r(1:sizeDotState), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo From 604bcd1229b2bdda91c89c998b457f3341471342 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 09:38:08 +0200 Subject: [PATCH 02/14] memory efficient FPI state integrator --- src/crystallite.f90 | 30 +++++++++++++----------------- src/material.f90 | 8 -------- 2 files changed, 13 insertions(+), 25 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index ece4db020..26fa1b3f8 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1016,11 +1016,13 @@ subroutine integrateStateFPI zeta real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & r ! state residuum + real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 + type(group_float), dimension(maxval(phase_Nsources)) :: source_dotState_p1, source_dotState_p2 logical :: & nonlocalBroken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState_p1, source_dotState_p2) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1045,24 +1047,22 @@ subroutine integrateStateFPI plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) + plastic_dotState_p2 = 0.0_pReal * plasticState(p)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) + source_dotState_p2(p)%p = 0.0_pReal * sourceState(p)%p(s)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer enddo iteration: do NiterationState = 1, num%nState - plasticState(p)%previousDotState2(:,c) = merge(plasticState(p)%previousDotState(:,c),& - 0.0_pReal,& - NiterationState > 1) - plasticState(p)%previousDotState (:,c) = plasticState(p)%dotState(:,c) + if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 + plastic_dotState_p1 = plasticState(p)%dotState(:,c) do s = 1, phase_Nsources(p) - sourceState(p)%p(s)%previousDotState2(:,c) = merge(sourceState(p)%p(s)%previousDotState(:,c),& - 0.0_pReal, & - NiterationState > 1) - sourceState(p)%p(s)%previousDotState (:,c) = sourceState(p)%p(s)%dotState(:,c) + if(nIterationState > 1) source_dotState_p2(p)%p = source_dotState_p1(p)%p + source_dotState_p1(p)%p = sourceState(p)%p(s)%dotState(:,c) enddo call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & @@ -1088,11 +1088,9 @@ subroutine integrateStateFPI if(.not. crystallite_todo(g,i,e)) exit iteration sizeDotState = plasticState(p)%sizeDotState - zeta = damper(plasticState(p)%dotState (:,c), & - plasticState(p)%previousDotState (:,c), & - plasticState(p)%previousDotState2(:,c)) + zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & - + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) + + plastic_dotState_p1 * (1.0_pReal - zeta) r(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - plasticState(p)%subState0(1:sizeDotState,c) & - plasticState(p)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) @@ -1103,11 +1101,9 @@ subroutine integrateStateFPI plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - zeta = damper(sourceState(p)%p(s)%dotState (:,c), & - sourceState(p)%p(s)%previousDotState (:,c), & - sourceState(p)%p(s)%previousDotState2(:,c)) + zeta = damper(sourceState(p)%p(s)%dotState(:,c),source_dotState_p1(p)%p,source_dotState_p2(p)%p) sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - zeta) + + source_dotState_p1(p)%p* (1.0_pReal - zeta) r(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) & - sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - sourceState(p)%p(s)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) diff --git a/src/material.f90 b/src/material.f90 index 4461ca958..14d4e4c19 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -724,10 +724,6 @@ subroutine material_allocatePlasticState(phase,NipcMyPhase,& allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 1) then - allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(plasticState(phase)%previousDotState2 (sizeDotState,NipcMyPhase),source=0.0_pReal) - endif if (numerics_integrator == 4) & allocate(plasticState(phase)%RK4dotState (4,sizeDotState,NipcMyPhase),source=0.0_pReal) if (numerics_integrator == 5) & @@ -762,10 +758,6 @@ subroutine material_allocateSourceState(phase,of,NipcMyPhase,& allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 1) then - allocate(sourceState(phase)%p(of)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - allocate(sourceState(phase)%p(of)%previousDotState2 (sizeDotState,NipcMyPhase),source=0.0_pReal) - endif if (numerics_integrator == 4) & allocate(sourceState(phase)%p(of)%RK4dotState (4,sizeDotState,NipcMyPhase),source=0.0_pReal) if (numerics_integrator == 5) & From 28dadb4422d3898d0ba2420aa946d8853dcf330d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 09:49:41 +0200 Subject: [PATCH 03/14] no need to check multiple times --- src/crystallite.f90 | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 26fa1b3f8..5c6510864 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1070,8 +1070,6 @@ subroutine integrateStateFPI g, i, e) crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & - nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit iteration call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1083,8 +1081,6 @@ subroutine integrateStateFPI do s = 1, phase_Nsources(p) crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) enddo - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & - nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit iteration sizeDotState = plasticState(p)%sizeDotState @@ -1117,12 +1113,12 @@ subroutine integrateStateFPI if(crystallite_converged(g,i,e)) then crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & - nonlocalBroken = .true. exit iteration endif enddo iteration + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + nonlocalBroken = .true. endif enddo; enddo; enddo From 2cc0c746d31f3ed784c8664cff3389e45da89344 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 10:07:01 +0200 Subject: [PATCH 04/14] no variable needed --- src/crystallite.f90 | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 5c6510864..9c4eeb326 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1540,15 +1540,6 @@ subroutine integrateStateRKCK45 logical :: & nonlocalBroken - real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & - homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - residuum_plastic - real(pReal), dimension(constitutive_source_maxSizeDotState, & - maxval(phase_Nsources), & - homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - residuum_source - - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -1629,13 +1620,12 @@ subroutine integrateStateRKCK45 sizeDotState = plasticState(p)%sizeDotState plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c) - residuum_plastic(1:sizeDotState,g,i,e) = matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) & - * crystallite_subdt(g,i,e) plasticState(p)%dotState(:,c) = matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & + crystallite_todo(g,i,e) = converged(matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) & + * crystallite_subdt(g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) @@ -1643,14 +1633,13 @@ subroutine integrateStateRKCK45 sizeDotState = sourceState(p)%p(s)%sizeDotState sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c) - residuum_source(1:sizeDotState,s,g,i,e) = matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) & - * crystallite_subdt(g,i,e) sourceState(p)%p(s)%dotState(:,c) = matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & - converged(residuum_source(1:sizeDotState,s,g,i,e), & + converged(matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) & + * crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo From 505c1432b1f61bc0718f7780110b722d7a4a2684 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 10:12:25 +0200 Subject: [PATCH 05/14] no need for temp storage --- src/crystallite.f90 | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 9c4eeb326..52d102536 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -804,7 +804,6 @@ logical function integrateStress(ipc,ip,el,timeFraction) residuumLi_old, & ! last residuum of intermediate velocity gradient deltaLi, & ! direction of next guess Fe, & ! elastic deformation gradient - Fe_new, & S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration A, & B, & @@ -982,17 +981,15 @@ logical function integrateStress(ipc,ip,el,timeFraction) invFp_new = matmul(invFp_current,B) call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error - Fp_new = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize - Fe_new = matmul(matmul(F,invFp_new),invFi_new) integrateStress = .true. crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) crystallite_S (1:3,1:3,ipc,ip,el) = S crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess crystallite_Li (1:3,1:3,ipc,ip,el) = Liguess - crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new + crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new - crystallite_Fe (1:3,1:3,ipc,ip,el) = Fe_new + crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) end function integrateStress From 5b29af847349b8df64e40c1a5975be7092c004f8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 10:16:14 +0200 Subject: [PATCH 06/14] clearer --- src/crystallite.f90 | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 52d102536..4df98c429 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -910,21 +910,19 @@ logical function integrateStress(ipc,ip,el,timeFraction) cycle LpLoop endif - !* calculate Jacobian for correction term - if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then + calculateJacobiLi: if (mod(jacoCounterLp, num%iJacoLpresiduum) == 0) then jacoCounterLp = jacoCounterLp + 1 do o=1,3; do p=1,3 - dFe_dLp(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFe_dLp(o,1:3,p,1:3) = - dt * A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) enddo; enddo - dFe_dLp = - dt * dFe_dLp dRLp_dLp = math_identity2nd(9) & - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) temp_9 = math_33to9(residuumLp) call dgesv(9,1,dRLp_dLp,9,devNull_9,temp_9,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp if (ierr /= 0) return ! error deltaLp = - math_9to33(temp_9) - endif + endif calculateJacobiLi Lpguess = Lpguess & + deltaLp * steplengthLp @@ -952,8 +950,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) cycle LiLoop endif - !* calculate Jacobian for correction term - if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then + calculateJacobiLp: if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then jacoCounterLi = jacoCounterLi + 1 temp_33 = matmul(matmul(A,B),invFi_current) @@ -972,7 +969,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) call dgesv(9,1,dRLi_dLi,9,devNull_9,temp_9,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li if (ierr /= 0) return ! error deltaLi = - math_9to33(temp_9) - endif + endif calculateJacobiLp Liguess = Liguess & + deltaLi * steplengthLi From 2de4b87c6197b8565f83e425e2cf0ee78117cf2c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 11:24:13 +0200 Subject: [PATCH 07/14] bugfix for FPI: loop over sources! no memory waste for adaptive Euler --- src/crystallite.f90 | 39 ++++++++++++++------------------------- 1 file changed, 14 insertions(+), 25 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 4df98c429..ae5091978 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1047,7 +1047,7 @@ subroutine integrateStateFPI sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - source_dotState_p2(p)%p = 0.0_pReal * sourceState(p)%p(s)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer + source_dotState_p2(s)%p = 0.0_pReal * sourceState(p)%p(s)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer enddo iteration: do NiterationState = 1, num%nState @@ -1055,8 +1055,8 @@ subroutine integrateStateFPI if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 plastic_dotState_p1 = plasticState(p)%dotState(:,c) do s = 1, phase_Nsources(p) - if(nIterationState > 1) source_dotState_p2(p)%p = source_dotState_p1(p)%p - source_dotState_p1(p)%p = sourceState(p)%p(s)%dotState(:,c) + if(nIterationState > 1) source_dotState_p2(s)%p = source_dotState_p1(s)%p + source_dotState_p1(s)%p = sourceState(p)%p(s)%dotState(:,c) enddo call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & @@ -1091,9 +1091,9 @@ subroutine integrateStateFPI plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - zeta = damper(sourceState(p)%p(s)%dotState(:,c),source_dotState_p1(p)%p,source_dotState_p2(p)%p) + zeta = damper(sourceState(p)%p(s)%dotState(:,c),source_dotState_p1(s)%p,source_dotState_p2(s)%p) sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState_p1(p)%p* (1.0_pReal - zeta) + + source_dotState_p1(s)%p* (1.0_pReal - zeta) r(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) & - sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - sourceState(p)%p(s)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) @@ -1236,14 +1236,8 @@ subroutine integrateStateAdaptiveEuler logical :: & nonlocalBroken - real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & - homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - residuum_plastic - real(pReal), dimension(constitutive_source_maxSizeDotState,& - maxval(phase_Nsources), & - homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & - residuum_source - + real(pReal), dimension(:), allocatable :: residuum_plastic + type(group_float), dimension(maxval(phase_Nsources)) :: residuum_source nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) @@ -1269,15 +1263,14 @@ subroutine integrateStateAdaptiveEuler sizeDotState = plasticState(p)%sizeDotState - residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) & - * (- 0.5_pReal * crystallite_subdt(g,i,e)) + residuum_plastic = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - residuum_source(1:sizeDotState,s,g,i,e) = sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * (- 0.5_pReal * crystallite_subdt(g,i,e)) + residuum_source(s)%p = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & + * 0.5_pReal * crystallite_subdt(g,i,e) sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo @@ -1312,21 +1305,17 @@ subroutine integrateStateAdaptiveEuler sizeDotState = plasticState(p)%sizeDotState - residuum_plastic(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) & - + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e) - - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & + crystallite_converged(g,i,e) = converged(residuum_plastic & + + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - residuum_source(1:sizeDotState,s,g,i,e) = & - residuum_source(1:sizeDotState,s,g,i,e) + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e) - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s,g,i,e), & + crystallite_converged(g,i,e) .and. converged(residuum_source(s)%p & + + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo From 00be291fa02ce5af06641efa6adc506509ff580c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 12:14:54 +0200 Subject: [PATCH 08/14] todo will be reset after state integration --- src/crystallite.f90 | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index ae5091978..0555e28a6 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1118,8 +1118,7 @@ subroutine integrateStateFPI enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false. - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck contains @@ -1214,8 +1213,7 @@ subroutine integrateStateEuler enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false. - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateEuler @@ -1324,7 +1322,7 @@ subroutine integrateStateAdaptiveEuler enddo; enddo; enddo !$OMP END PARALLEL DO - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateAdaptiveEuler @@ -1478,8 +1476,7 @@ subroutine integrateStateRK4 enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false. - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateRK4 @@ -1652,8 +1649,7 @@ subroutine integrateStateRKCK45 enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) where(.not. crystallite_localPlasticity) crystallite_todo = .false. - if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + if (nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateRKCK45 @@ -1664,8 +1660,7 @@ end subroutine integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- subroutine nonlocalConvergenceCheck - if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... - where( .not. crystallite_localPlasticity) crystallite_converged = .false. + where( .not. crystallite_localPlasticity) crystallite_converged = .false. end subroutine nonlocalConvergenceCheck From e818dfdb3e5c6d9a4bd44a05291b951ffa02cfee Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 20:42:51 +0200 Subject: [PATCH 09/14] not used anymore --- src/prec.f90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/prec.f90 b/src/prec.f90 index a8d7d9e33..974c63c2f 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -30,10 +30,6 @@ module prec real(pReal), dimension(:), pointer :: p end type group_float - type :: group_int - integer, dimension(:), pointer :: p - end type group_int - ! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array type :: tState integer :: & From 9c95ce36f4b6ef1a783fce150481501f889891e0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 20:57:09 +0200 Subject: [PATCH 10/14] automatic LHS (re)-allocation does not work for pointers group_float has pointers, not allocatables --- src/crystallite.f90 | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 5c53e9165..d547a6a01 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1011,7 +1011,7 @@ subroutine integrateStateFPI real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & r ! state residuum real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 - type(group_float), dimension(maxval(phase_Nsources)) :: source_dotState_p1, source_dotState_p2 + real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & nonlocalBroken @@ -1047,7 +1047,7 @@ subroutine integrateStateFPI sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - source_dotState_p2(s)%p = 0.0_pReal * sourceState(p)%p(s)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer + source_dotState(1:sizeDotState,2,s) = 0.0_pReal enddo iteration: do NiterationState = 1, num%nState @@ -1055,8 +1055,9 @@ subroutine integrateStateFPI if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 plastic_dotState_p1 = plasticState(p)%dotState(:,c) do s = 1, phase_Nsources(p) - if(nIterationState > 1) source_dotState_p2(s)%p = source_dotState_p1(s)%p - source_dotState_p1(s)%p = sourceState(p)%p(s)%dotState(:,c) + sizeDotState = sourceState(p)%p(s)%sizeDotState + if(nIterationState > 1) source_dotState(1:sizeDotState,2,s) = source_dotState(1:sizeDotState,1,s) + source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) enddo call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & @@ -1091,9 +1092,11 @@ subroutine integrateStateFPI plasticState(p)%atol(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - zeta = damper(sourceState(p)%p(s)%dotState(:,c),source_dotState_p1(s)%p,source_dotState_p2(s)%p) + zeta = damper(sourceState(p)%p(s)%dotState(:,c), & + source_dotState(1:sizeDotState,1,s),& + source_dotState(1:sizeDotState,2,s)) sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState_p1(s)%p* (1.0_pReal - zeta) + + source_dotState(1:sizeDotState,1,s)* (1.0_pReal - zeta) r(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) & - sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - sourceState(p)%p(s)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) From e46220cd8a471039438d967fae0dd698e97a727e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 31 Mar 2020 22:02:53 +0200 Subject: [PATCH 11/14] OMP bugfix for FPI integrator, memory-efficient RK4 --- src/crystallite.f90 | 31 ++++++++++++++++++------------- src/material.f90 | 4 ---- src/prec.f90 | 5 +---- 3 files changed, 19 insertions(+), 21 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d547a6a01..9617cc283 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1016,7 +1016,7 @@ subroutine integrateStateFPI nonlocalBroken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState_p1, source_dotState_p2) + !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1359,8 +1359,10 @@ subroutine integrateStateRK4 logical :: & nonlocalBroken + real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1382,20 +1384,23 @@ subroutine integrateStateRK4 if(.not. crystallite_todo(g,i,e)) cycle do stage = 1,3 - - plasticState(p)%RK4dotState(stage,:,c) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RK4dotState(1,:,c) + sizeDotState = plasticState(p)%sizeDotState + plastic_RK4dotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) + plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RK4dotState(1:sizeDotState,1) do s = 1, phase_Nsources(p) - sourceState(p)%p(s)%RK4dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RK4dotState(1,:,c) + sizeDotState = sourceState(p)%p(s)%sizeDotState + source_RK4dotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RK4dotState(1:sizeDotState,1,s) enddo do n = 2, stage + sizeDotState = plasticState(p)%sizeDotState plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + A(n,stage) * plasticState(p)%RK4dotState(n,:,c) + + A(n,stage) * plastic_RK4dotState(1:sizeDotState,n) do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & - + A(n,stage) * sourceState(p)%p(s)%RK4dotState(n,:,c) + + A(n,stage) * source_RK4dotState(1:sizeDotState,n,s) enddo enddo @@ -1438,9 +1443,9 @@ subroutine integrateStateRK4 sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%RK4dotState(4,:,c) = plasticState (p)%dotState(:,c) + plastic_RK4dotState(1:sizeDotState,4) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(B,plasticState(p)%RK4dotState(1:4,1:sizeDotState,c)) + plasticState(p)%dotState(:,c) = matmul(plastic_RK4dotState(1:sizeDotState,1:4),B) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) @@ -1448,9 +1453,9 @@ subroutine integrateStateRK4 do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%RK4dotState(4,:,c) = sourceState(p)%p(s)%dotState(:,c) + source_RK4dotState(1:sizeDotState,4,s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = matmul(B,sourceState(p)%p(s)%RK4dotState(1:4,1:sizeDotState,c)) + sourceState(p)%p(s)%dotState(:,c) = matmul(source_RK4dotState(1:sizeDotState,1:4,s),B) sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) diff --git a/src/material.f90 b/src/material.f90 index 14d4e4c19..fc55bb17c 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -724,8 +724,6 @@ subroutine material_allocatePlasticState(phase,NipcMyPhase,& allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 4) & - allocate(plasticState(phase)%RK4dotState (4,sizeDotState,NipcMyPhase),source=0.0_pReal) if (numerics_integrator == 5) & allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) @@ -758,8 +756,6 @@ subroutine material_allocateSourceState(phase,of,NipcMyPhase,& allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 4) & - allocate(sourceState(phase)%p(of)%RK4dotState (4,sizeDotState,NipcMyPhase),source=0.0_pReal) if (numerics_integrator == 5) & allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) diff --git a/src/prec.f90 b/src/prec.f90 index 974c63c2f..0afaeb536 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -46,11 +46,8 @@ module prec deltaState !< increment of state change real(pReal), allocatable, dimension(:,:) :: & partionedState0, & - subState0, & - previousDotState, & - previousDotState2 + subState0 real(pReal), allocatable, dimension(:,:,:) :: & - RK4dotState, & RKCK45dotState end type From 6bfa51c30777d6a9ccab1898af32a27ec32f6c20 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 07:02:23 +0200 Subject: [PATCH 12/14] LHS allocation does not work for pointers --- src/crystallite.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 9617cc283..5884e2cef 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1237,11 +1237,11 @@ subroutine integrateStateAdaptiveEuler logical :: & nonlocalBroken - real(pReal), dimension(:), allocatable :: residuum_plastic - type(group_float), dimension(maxval(phase_Nsources)) :: residuum_source + real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic + real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1264,14 +1264,14 @@ subroutine integrateStateAdaptiveEuler sizeDotState = plasticState(p)%sizeDotState - residuum_plastic = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) + residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - residuum_source(s)%p = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * 0.5_pReal * crystallite_subdt(g,i,e) + residuum_source(1:sizeDotState,s) = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & + * 0.5_pReal * crystallite_subdt(g,i,e) sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo @@ -1306,7 +1306,7 @@ subroutine integrateStateAdaptiveEuler sizeDotState = plasticState(p)%sizeDotState - crystallite_converged(g,i,e) = converged(residuum_plastic & + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) @@ -1315,7 +1315,7 @@ subroutine integrateStateAdaptiveEuler sizeDotState = sourceState(p)%p(s)%sizeDotState crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(s)%p & + crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) From 54e3455bd47a31a9e82b63b7b7e6be2bcb04e313 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 07:09:33 +0200 Subject: [PATCH 13/14] RKCK45 store intermediate state only per point --- src/crystallite.f90 | 33 +++++++++++++++++++-------------- src/material.f90 | 4 ---- src/prec.f90 | 2 -- 3 files changed, 19 insertions(+), 20 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 5884e2cef..fb12c439b 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1526,9 +1526,11 @@ subroutine integrateStateRKCK45 sizeDotState logical :: & nonlocalBroken + real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1550,20 +1552,23 @@ subroutine integrateStateRKCK45 if(.not. crystallite_todo(g,i,e)) cycle do stage = 1,5 - - plasticState(p)%RKCK45dotState(stage,:,c) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,c) + sizeDotState = plasticState(p)%sizeDotState + plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) + plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) do s = 1, phase_Nsources(p) - sourceState(p)%p(s)%RKCK45dotState(stage,:,c) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,c) + sizeDotState = sourceState(p)%p(s)%sizeDotState + source_RKdotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RKdotState(1:sizeDotState,1,s) enddo do n = 2, stage + sizeDotState = plasticState(p)%sizeDotState plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,c) + + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & - + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,c) + + A(n,stage) * source_RKdotState(1:sizeDotState,n,s) enddo enddo @@ -1606,12 +1611,12 @@ subroutine integrateStateRKCK45 sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%RKCK45dotState(6,:,c) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(B,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) + plastic_RKdotState(1:sizeDotState,6) = plasticState (p)%dotState(:,c) + plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:6),B) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - crystallite_todo(g,i,e) = converged(matmul(DB,plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,c)) & + crystallite_todo(g,i,e) = converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) & * crystallite_subdt(g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) @@ -1619,13 +1624,13 @@ subroutine integrateStateRKCK45 do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%RKCK45dotState(6,:,c) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = matmul(B,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) + source_RKdotState(1:sizeDotState,6,s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:6,s),B) sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & - converged(matmul(DB,sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,c)) & + converged(matmul(source_RKdotState(1:sizeDotState,1:6,s),DB) & * crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) diff --git a/src/material.f90 b/src/material.f90 index fc55bb17c..aefe44878 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -724,8 +724,6 @@ subroutine material_allocatePlasticState(phase,NipcMyPhase,& allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 5) & - allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) @@ -756,8 +754,6 @@ subroutine material_allocateSourceState(phase,of,NipcMyPhase,& allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - if (numerics_integrator == 5) & - allocate(sourceState(phase)%p(of)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) diff --git a/src/prec.f90 b/src/prec.f90 index 0afaeb536..dd39ddc05 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -47,8 +47,6 @@ module prec real(pReal), allocatable, dimension(:,:) :: & partionedState0, & subState0 - real(pReal), allocatable, dimension(:,:,:) :: & - RKCK45dotState end type type, extends(tState) :: tPlasticState From 5e9ff7947bfbd8ef1d122cc412b517fef7df9b75 Mon Sep 17 00:00:00 2001 From: Franz Roters Date: Fri, 3 Apr 2020 14:01:35 +0200 Subject: [PATCH 14/14] [skip ci] plastic_dotstate always before source_dotstate --- src/crystallite.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index fb12c439b..efa066696 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1359,10 +1359,10 @@ subroutine integrateStateRK4 logical :: & nonlocalBroken - real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState + real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RK4dotState,source_RK4dotState) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1526,8 +1526,8 @@ subroutine integrateStateRKCK45 sizeDotState logical :: & nonlocalBroken - real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState + real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState)