From 73f39136c48b5b33f7f515b545dd36ac54c26d32 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Jan 2019 22:19:38 +0100 Subject: [PATCH 01/31] taking over from old branch --- src/crystallite.f90 | 37 ++++++++++++++----------------------- 1 file changed, 14 insertions(+), 23 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 74eef259e..b720c4101 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1551,7 +1551,6 @@ subroutine integrateStateFPI() homogenization_Ngrains use constitutive, only: & constitutive_collectDotState, & - constitutive_microstructure, & constitutive_plasticity_maxSizeDotState, & constitutive_source_maxSizeDotState @@ -1569,9 +1568,9 @@ subroutine integrateStateFPI() real(pReal) :: & stateDamper real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & - plasticStateResiduum + residuum_plastic ! residuum for plastic state real(pReal), dimension(constitutive_source_maxSizeDotState, maxval(phase_Nsources)) :: & - sourceStateResiduum + residuum_source ! residuum for source state logical :: & converged, & doneWithIntegration @@ -1616,7 +1615,7 @@ subroutine integrateStateFPI() !$OMP PARALLEL !$OMP DO PRIVATE(sizeDotState, & - !$OMP& plasticStateResiduum,sourceStateResiduum, & + !$OMP& residuum_plastic,residuum_source, & !$OMP& stateDamper, converged,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) @@ -1629,21 +1628,20 @@ subroutine integrateStateFPI() plasticState(p)%previousDotState2(:,c)) sizeDotState = plasticState(p)%sizeDotState - plasticStateResiduum(1:sizeDotState) = plasticState(p)%state (1:sizeDotState,c) & + residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - plasticState(p)%subState0(1:sizeDotState,c) & - ( plasticState(p)%dotState (:,c) * stateDamper & + plasticState(p)%previousDotState(:,c) * (1.0_pReal-stateDamper) & ) * crystallite_subdt(g,i,e) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & - - plasticStateResiduum(1:sizeDotState) - + - residuum_plastic(1:sizeDotState) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * stateDamper & + plasticState(p)%previousDotState(:,c) * (1.0_pReal - stateDamper) - converged = all( abs(plasticStateResiduum(1:sizeDotState)) < & + converged = all( abs(residuum_plastic(1:sizeDotState)) < & plasticState(p)%aTolState(1:sizeDotState) & - .or. abs(plasticStateResiduum(1:sizeDotState)) < & + .or. abs(residuum_plastic(1:sizeDotState)) < & rTol_crystalliteState * abs( plasticState(p)%state(1:sizeDotState,c))) @@ -1652,26 +1650,21 @@ subroutine integrateStateFPI() sourceState(p)%p(s)%previousDotState (:,c), & sourceState(p)%p(s)%previousDotState2(:,c)) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceStateResiduum(1:sizeDotState,s) = sourceState(p)%p(s)%state (1:sizeDotState,c) & + residuum_source(1:sizeDotState,s) = sourceState(p)%p(s)%state (1:sizeDotState,c) & - sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - ( sourceState(p)%p(s)%dotState (:,c) * stateDamper & + sourceState(p)%p(s)%previousDotState(:,c) * (1.0_pReal - stateDamper) & ) * crystallite_subdt(g,i,e) - ! --- correct state with residuum --- - sourceState(p)%p(s)%state(1:sizeDotState,c) = & - sourceState(p)%p(s)%state(1:sizeDotState,c) & - - sourceStateResiduum(1:sizeDotState,s) ! need to copy to local variable, since we cant flush a pointer in openmp - - sourceState(p)%p(s)%dotState(:,c) = & - sourceState(p)%p(s)%dotState(:,c) * stateDamper & - + sourceState(p)%p(s)%previousDotState(:,c) & - * (1.0_pReal - stateDamper) + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) & + - residuum_source(1:sizeDotState,s) + sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * stateDamper & + + sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - stateDamper) converged = converged .and. & - all( abs(sourceStateResiduum(1:sizeDotState,s)) < & + all( abs(residuum_source(1:sizeDotState,s)) < & sourceState(p)%p(s)%aTolState(1:sizeDotState) & - .or. abs(sourceStateResiduum(1:sizeDotState,s)) < & + .or. abs(residuum_source(1:sizeDotState,s)) < & rTol_crystalliteState * abs(sourceState(p)%p(s)%state(1:sizeDotState,c))) enddo if (converged) crystallite_converged(g,i,e) = .true. ! ... converged per definition @@ -1771,8 +1764,6 @@ end subroutine integrateStateEuler !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- subroutine integrateStateAdaptiveEuler() - use, intrinsic :: & - IEEE_arithmetic use numerics, only: & rTol_crystalliteState use mesh, only: & From b4afc303be3b2cdbd98bbc629d413cd96f3c17c3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Jan 2019 22:29:36 +0100 Subject: [PATCH 02/31] clearer logic --- src/crystallite.f90 | 96 ++++++++++++++++++++++----------------------- 1 file changed, 46 insertions(+), 50 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index b720c4101..be14f801a 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1566,13 +1566,12 @@ subroutine integrateStateFPI() s, & sizeDotState real(pReal) :: & - stateDamper + zeta real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & residuum_plastic ! residuum for plastic state - real(pReal), dimension(constitutive_source_maxSizeDotState, maxval(phase_Nsources)) :: & + real(pReal), dimension(constitutive_source_maxSizeDotState) :: & residuum_source ! residuum for source state logical :: & - converged, & doneWithIntegration ! --+>> PREGUESS FOR STATE <<+-- @@ -1614,65 +1613,59 @@ subroutine integrateStateFPI() call update_dotState(1.0_pReal) !$OMP PARALLEL - !$OMP DO PRIVATE(sizeDotState, & - !$OMP& residuum_plastic,residuum_source, & - !$OMP& stateDamper, converged,p,c) + !$OMP DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then p = phaseAt(g,i,e); c = phasememberAt(g,i,e) - StateDamper = damper(plasticState(p)%dotState (:,c), & - plasticState(p)%previousDotState (:,c), & - plasticState(p)%previousDotState2(:,c)) + zeta = damper(plasticState(p)%dotState (:,c), & + plasticState(p)%previousDotState (:,c), & + plasticState(p)%previousDotState2(:,c)) sizeDotState = plasticState(p)%sizeDotState residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - plasticState(p)%subState0(1:sizeDotState,c) & - - ( plasticState(p)%dotState (:,c) * stateDamper & - + plasticState(p)%previousDotState(:,c) * (1.0_pReal-stateDamper) & + - ( plasticState(p)%dotState (:,c) * zeta & + + plasticState(p)%previousDotState(:,c) * (1.0_pReal-zeta) & ) * crystallite_subdt(g,i,e) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & - residuum_plastic(1:sizeDotState) - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * stateDamper & - + plasticState(p)%previousDotState(:,c) * (1.0_pReal - stateDamper) - - converged = all( abs(residuum_plastic(1:sizeDotState)) < & - plasticState(p)%aTolState(1:sizeDotState) & - .or. abs(residuum_plastic(1:sizeDotState)) < & - rTol_crystalliteState * abs( plasticState(p)%state(1:sizeDotState,c))) + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) + + crystallite_converged(g,i,e) = all(abs(residuum_plastic(1:sizeDotState)) & + < min(plasticState(p)%aTolState(1:sizeDotState), & + abs(plasticState(p)%state(1:sizeDotState,c)*rTol_crystalliteState))) do s = 1_pInt, phase_Nsources(p) - stateDamper = 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), & + sourceState(p)%p(s)%previousDotState (:,c), & + sourceState(p)%p(s)%previousDotState2(:,c)) sizeDotState = sourceState(p)%p(s)%sizeDotState - residuum_source(1:sizeDotState,s) = sourceState(p)%p(s)%state (1:sizeDotState,c) & - - sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - - ( sourceState(p)%p(s)%dotState (:,c) * stateDamper & - + sourceState(p)%p(s)%previousDotState(:,c) * (1.0_pReal - stateDamper) & - ) * crystallite_subdt(g,i,e) + + 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 (:,c) * zeta & + + sourceState(p)%p(s)%previousDotState(:,c) * (1.0_pReal - zeta) & + ) * 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,s) - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * stateDamper & - + sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - stateDamper) + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) & + - residuum_source(1:sizeDotState) + sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & + + sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - zeta) - converged = converged .and. & - all( abs(residuum_source(1:sizeDotState,s)) < & - sourceState(p)%p(s)%aTolState(1:sizeDotState) & - .or. abs(residuum_source(1:sizeDotState,s)) < & - rTol_crystalliteState * abs(sourceState(p)%p(s)%state(1:sizeDotState,c))) - enddo - if (converged) crystallite_converged(g,i,e) = .true. ! ... converged per definition - - endif - enddo; enddo; enddo + crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and. & + all(abs(residuum_source(1:sizeDotState)) & + < min(sourceState(p)%p(s)%aTolState(1:sizeDotState), & + abs(sourceState(p)%p(s)%state(1:sizeDotState,c)*rTol_crystalliteState))) + enddo + endif + enddo; enddo; enddo !$OMP ENDDO - ! --- STATE JUMP --- !$OMP DO do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -1870,6 +1863,17 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & plasticStateResiduum(1:mySizePlasticDotState,g,i,e) & + 0.5_pReal * plasticState(p)%dotState(:,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state + + converged = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & + rTol_crystalliteState .or. & + abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & + plasticState(p)%aTolState(1:mySizePlasticDotState)) + + forall (s = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%dotState(s,c)) > 0.0_pReal) & + relPlasticStateResiduum(s,g,i,e) = & + plasticStateResiduum(s,g,i,e) / plasticState(p)%dotState(s,c) + + do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e) = & @@ -1878,10 +1882,7 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state enddo - ! --- relative residui --- - forall (s = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%dotState(s,c)) > 0.0_pReal) & - relPlasticStateResiduum(s,g,i,e) = & - plasticStateResiduum(s,g,i,e) / plasticState(p)%dotState(s,c) + do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState forall (s = 1_pInt:mySizeSourceDotState,abs(sourceState(p)%p(mySource)%dotState(s,c)) > 0.0_pReal) & @@ -1889,11 +1890,6 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & sourceStateResiduum(s,mySource,g,i,e) / sourceState(p)%p(mySource)%dotState(s,c) enddo - ! --- converged ? --- - converged = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & - rTol_crystalliteState .or. & - abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & - plasticState(p)%aTolState(1:mySizePlasticDotState)) do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState converged = converged .and. & From 0be05b3ee1c4a894b2e442ff3156945bb3a5efe5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Jan 2019 22:46:21 +0100 Subject: [PATCH 03/31] one variable is enough --- src/crystallite.f90 | 63 ++++++++++++++++++++++----------------------- 1 file changed, 31 insertions(+), 32 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index be14f801a..100bd1aa4 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1785,8 +1785,7 @@ subroutine integrateStateAdaptiveEuler() p, & c, & mySource, & - mySizePlasticDotState, & ! size of dot states - mySizeSourceDotState + sizeDotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & plasticStateResiduum, & ! residuum from evolution in micrstructure @@ -1810,31 +1809,31 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & ! contribution to state and relative residui and from Euler integration call update_dotState(1.0_pReal) - !$OMP PARALLEL DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then p = phaseAt(g,i,e); c = phasememberAt(g,i,e) - mySizePlasticDotState = plasticState(p)%sizeDotState - plasticStateResiduum(1:mySizePlasticDotState,g,i,e) = & + sizeDotState = plasticState(p)%sizeDotState + plasticStateResiduum(1:sizeDotState,g,i,e) = & - 0.5_pReal & - * plasticState(p)%dotstate(1:mySizePlasticDotState,c) & + * plasticState(p)%dotstate(1:sizeDotState,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - plasticState(p)%state (1:mySizePlasticDotState,c) = & - plasticState(p)%state (1:mySizePlasticDotState,c) & - + plasticState(p)%dotstate(1:mySizePlasticDotState,c) & + plasticState(p)%state (1:sizeDotState,c) = & + plasticState(p)%state (1:sizeDotState,c) & + + plasticState(p)%dotstate(1:sizeDotState,c) & * crystallite_subdt(g,i,e) do mySource = 1_pInt, phase_Nsources(p) - mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState - sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e) = & + sizeDotState = sourceState(p)%p(mySource)%sizeDotState + sourceStateResiduum(1:sizeDotState,mySource,g,i,e) = & - 0.5_pReal & - * sourceState(p)%p(mySource)%dotstate(1:mySizeSourceDotState,c) & + * sourceState(p)%p(mySource)%dotstate(1:sizeDotState,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - sourceState(p)%p(mySource)%state (1:mySizeSourceDotState,c) = & - sourceState(p)%p(mySource)%state (1:mySizeSourceDotState,c) & - + sourceState(p)%p(mySource)%dotstate(1:mySizeSourceDotState,c) & + sourceState(p)%p(mySource)%state (1:sizeDotState,c) = & + sourceState(p)%p(mySource)%state (1:sizeDotState,c) & + + sourceState(p)%p(mySource)%dotstate(1:sizeDotState,c) & * crystallite_subdt(g,i,e) enddo endif @@ -1850,7 +1849,7 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & relSourceStateResiduum = 0.0_pReal - !$OMP PARALLEL DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,converged,p,c,s) + !$OMP PARALLEL DO PRIVATE(sizeDotState,converged,p,c,s) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) @@ -1858,45 +1857,45 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & p = phaseAt(g,i,e); c = phasememberAt(g,i,e) ! --- contribution of heun step to absolute residui --- - mySizePlasticDotState = plasticState(p)%sizeDotState - plasticStateResiduum(1:mySizePlasticDotState,g,i,e) = & - plasticStateResiduum(1:mySizePlasticDotState,g,i,e) & + sizeDotState = plasticState(p)%sizeDotState + plasticStateResiduum(1:sizeDotState,g,i,e) = & + plasticStateResiduum(1:sizeDotState,g,i,e) & + 0.5_pReal * plasticState(p)%dotState(:,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - converged = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & + converged = all(abs(relPlasticStateResiduum(1:sizeDotState,g,i,e)) < & rTol_crystalliteState .or. & - abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & - plasticState(p)%aTolState(1:mySizePlasticDotState)) + abs(plasticStateResiduum(1:sizeDotState,g,i,e)) < & + plasticState(p)%aTolState(1:sizeDotState)) - forall (s = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%dotState(s,c)) > 0.0_pReal) & + forall (s = 1_pInt:sizeDotState, abs(plasticState(p)%dotState(s,c)) > 0.0_pReal) & relPlasticStateResiduum(s,g,i,e) = & plasticStateResiduum(s,g,i,e) / plasticState(p)%dotState(s,c) do mySource = 1_pInt, phase_Nsources(p) - mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState - sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e) = & - sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e) & + sizeDotState = sourceState(p)%p(mySource)%sizeDotState + sourceStateResiduum(1:sizeDotState,mySource,g,i,e) = & + sourceStateResiduum(1:sizeDotState,mySource,g,i,e) & + 0.5_pReal * sourceState(p)%p(mySource)%dotState(:,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state enddo do mySource = 1_pInt, phase_Nsources(p) - mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState - forall (s = 1_pInt:mySizeSourceDotState,abs(sourceState(p)%p(mySource)%dotState(s,c)) > 0.0_pReal) & + sizeDotState = sourceState(p)%p(mySource)%sizeDotState + forall (s = 1_pInt:sizeDotState,abs(sourceState(p)%p(mySource)%dotState(s,c)) > 0.0_pReal) & relSourceStateResiduum(s,mySource,g,i,e) = & sourceStateResiduum(s,mySource,g,i,e) / sourceState(p)%p(mySource)%dotState(s,c) enddo do mySource = 1_pInt, phase_Nsources(p) - mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState + sizeDotState = sourceState(p)%p(mySource)%sizeDotState converged = converged .and. & - all(abs(relSourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e)) < & + all(abs(relSourceStateResiduum(1:sizeDotState,mySource,g,i,e)) < & rTol_crystalliteState .or. & - abs(sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e)) < & - sourceState(p)%p(mySource)%aTolState(1:mySizeSourceDotState)) + abs(sourceStateResiduum(1:sizeDotState,mySource,g,i,e)) < & + sourceState(p)%p(mySource)%aTolState(1:sizeDotState)) enddo if (converged) crystallite_converged(g,i,e) = .true. ! ... converged per definition endif From 1408d66c0caec280768039732cb09ad53b579475 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Jan 2019 23:02:59 +0100 Subject: [PATCH 04/31] s is used for source --- src/crystallite.f90 | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 100bd1aa4..b416573de 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1781,7 +1781,7 @@ subroutine integrateStateAdaptiveEuler() e, & ! element index in element loop i, & ! integration point index in ip loop g, & ! grain index in grain loop - s, & ! state index + u, & ! state index p, & c, & mySource, & @@ -1849,7 +1849,7 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & relSourceStateResiduum = 0.0_pReal - !$OMP PARALLEL DO PRIVATE(sizeDotState,converged,p,c,s) + !$OMP PARALLEL DO PRIVATE(sizeDotState,converged,p,c,u) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) @@ -1868,9 +1868,9 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & abs(plasticStateResiduum(1:sizeDotState,g,i,e)) < & plasticState(p)%aTolState(1:sizeDotState)) - forall (s = 1_pInt:sizeDotState, abs(plasticState(p)%dotState(s,c)) > 0.0_pReal) & - relPlasticStateResiduum(s,g,i,e) = & - plasticStateResiduum(s,g,i,e) / plasticState(p)%dotState(s,c) + forall (u = 1_pInt:sizeDotState, abs(plasticState(p)%dotState(u,c)) > 0.0_pReal) & + relPlasticStateResiduum(u,g,i,e) = & + plasticStateResiduum(u,g,i,e) / plasticState(p)%dotState(u,c) do mySource = 1_pInt, phase_Nsources(p) @@ -1879,17 +1879,11 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & sourceStateResiduum(1:sizeDotState,mySource,g,i,e) & + 0.5_pReal * sourceState(p)%p(mySource)%dotState(:,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - enddo + forall (u = 1_pInt:sizeDotState,abs(sourceState(p)%p(mySource)%dotState(u,c)) > 0.0_pReal) & + relSourceStateResiduum(u,mySource,g,i,e) = & + sourceStateResiduum(u,mySource,g,i,e) / sourceState(p)%p(mySource)%dotState(u,c) - do mySource = 1_pInt, phase_Nsources(p) - sizeDotState = sourceState(p)%p(mySource)%sizeDotState - forall (s = 1_pInt:sizeDotState,abs(sourceState(p)%p(mySource)%dotState(s,c)) > 0.0_pReal) & - relSourceStateResiduum(s,mySource,g,i,e) = & - sourceStateResiduum(s,mySource,g,i,e) / sourceState(p)%p(mySource)%dotState(s,c) - enddo - - do mySource = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(mySource)%sizeDotState converged = converged .and. & all(abs(relSourceStateResiduum(1:sizeDotState,mySource,g,i,e)) < & From eade54a68f49c31e274577292970271efef2d915 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Jan 2019 23:04:50 +0100 Subject: [PATCH 05/31] consistent variable names --- src/crystallite.f90 | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index b416573de..054bb9d22 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1784,7 +1784,7 @@ subroutine integrateStateAdaptiveEuler() u, & ! state index p, & c, & - mySource, & + s, & sizeDotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & @@ -1825,15 +1825,15 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & plasticState(p)%state (1:sizeDotState,c) & + plasticState(p)%dotstate(1:sizeDotState,c) & * crystallite_subdt(g,i,e) - do mySource = 1_pInt, phase_Nsources(p) - sizeDotState = sourceState(p)%p(mySource)%sizeDotState - sourceStateResiduum(1:sizeDotState,mySource,g,i,e) = & + do s = 1_pInt, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceStateResiduum(1:sizeDotState,s,g,i,e) = & - 0.5_pReal & - * sourceState(p)%p(mySource)%dotstate(1:sizeDotState,c) & + * sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - sourceState(p)%p(mySource)%state (1:sizeDotState,c) = & - sourceState(p)%p(mySource)%state (1:sizeDotState,c) & - + sourceState(p)%p(mySource)%dotstate(1:sizeDotState,c) & + sourceState(p)%p(s)%state (1:sizeDotState,c) = & + sourceState(p)%p(s)%state (1:sizeDotState,c) & + + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & * crystallite_subdt(g,i,e) enddo endif @@ -1873,23 +1873,23 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & plasticStateResiduum(u,g,i,e) / plasticState(p)%dotState(u,c) - do mySource = 1_pInt, phase_Nsources(p) - sizeDotState = sourceState(p)%p(mySource)%sizeDotState - sourceStateResiduum(1:sizeDotState,mySource,g,i,e) = & - sourceStateResiduum(1:sizeDotState,mySource,g,i,e) & - + 0.5_pReal * sourceState(p)%p(mySource)%dotState(:,c) & + do s = 1_pInt, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceStateResiduum(1:sizeDotState,s,g,i,e) = & + sourceStateResiduum(1:sizeDotState,s,g,i,e) & + + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - forall (u = 1_pInt:sizeDotState,abs(sourceState(p)%p(mySource)%dotState(u,c)) > 0.0_pReal) & - relSourceStateResiduum(u,mySource,g,i,e) = & - sourceStateResiduum(u,mySource,g,i,e) / sourceState(p)%p(mySource)%dotState(u,c) + forall (u = 1_pInt:sizeDotState,abs(sourceState(p)%p(s)%dotState(u,c)) > 0.0_pReal) & + relSourceStateResiduum(u,s,g,i,e) = & + sourceStateResiduum(u,s,g,i,e) / sourceState(p)%p(s)%dotState(u,c) - sizeDotState = sourceState(p)%p(mySource)%sizeDotState + sizeDotState = sourceState(p)%p(s)%sizeDotState converged = converged .and. & - all(abs(relSourceStateResiduum(1:sizeDotState,mySource,g,i,e)) < & + all(abs(relSourceStateResiduum(1:sizeDotState,s,g,i,e)) < & rTol_crystalliteState .or. & - abs(sourceStateResiduum(1:sizeDotState,mySource,g,i,e)) < & - sourceState(p)%p(mySource)%aTolState(1:sizeDotState)) + abs(sourceStateResiduum(1:sizeDotState,s,g,i,e)) < & + sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo if (converged) crystallite_converged(g,i,e) = .true. ! ... converged per definition endif From bdd193fbd73eb9a1e2214684b90d4425edd94519 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Jan 2019 23:31:26 +0100 Subject: [PATCH 06/31] now readable (kind of) --- src/crystallite.f90 | 105 ++++++++++++++++++-------------------------- 1 file changed, 43 insertions(+), 62 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 054bb9d22..2e68ae756 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1619,12 +1619,12 @@ subroutine integrateStateFPI() do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then p = phaseAt(g,i,e); c = phasememberAt(g,i,e) + sizeDotState = plasticState(p)%sizeDotState zeta = damper(plasticState(p)%dotState (:,c), & plasticState(p)%previousDotState (:,c), & plasticState(p)%previousDotState2(:,c)) - sizeDotState = plasticState(p)%sizeDotState - + residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - plasticState(p)%subState0(1:sizeDotState,c) & - ( plasticState(p)%dotState (:,c) * zeta & @@ -1642,11 +1642,12 @@ subroutine integrateStateFPI() do s = 1_pInt, 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)) - sizeDotState = sourceState(p)%p(s)%sizeDotState - + 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 (:,c) * zeta & @@ -1771,8 +1772,6 @@ subroutine integrateStateAdaptiveEuler() phase_Nsources, & homogenization_maxNgrains use constitutive, only: & - constitutive_collectDotState, & - constitutive_microstructure, & constitutive_plasticity_maxSizeDotState, & constitutive_source_maxSizeDotState @@ -1786,6 +1785,8 @@ subroutine integrateStateAdaptiveEuler() c, & s, & sizeDotState + + ! ToDo: MD: once all constitutives use allocate state, attach these arrays to the state in case of adaptive Euler real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & plasticStateResiduum, & ! residuum from evolution in micrstructure @@ -1796,45 +1797,29 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & sourceStateResiduum, & ! residuum from evolution in micrstructure relSourceStateResiduum ! relative residuum from evolution in microstructure - logical :: & - converged - - - plasticStateResiduum = 0.0_pReal - relPlasticStateResiduum = 0.0_pReal - sourceStateResiduum = 0.0_pReal - relSourceStateResiduum = 0.0_pReal - !-------------------------------------------------------------------------------------------------- ! contribution to state and relative residui and from Euler integration call update_dotState(1.0_pReal) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,homogenization_Ngrains(mesh_element(3,e)) - if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - p = phaseAt(g,i,e); c = phasememberAt(g,i,e) - + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + do g = 1,homogenization_Ngrains(mesh_element(3,e)) + if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then + p = phaseAt(g,i,e); c = phasememberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - plasticStateResiduum(1:sizeDotState,g,i,e) = & - - 0.5_pReal & - * plasticState(p)%dotstate(1:sizeDotState,c) & - * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - plasticState(p)%state (1:sizeDotState,c) = & - plasticState(p)%state (1:sizeDotState,c) & - + plasticState(p)%dotstate(1:sizeDotState,c) & - * crystallite_subdt(g,i,e) + + plasticStateResiduum(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) & + * (- 0.5_pReal * crystallite_subdt(g,i,e)) + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & + + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) !ToDo: state, partitioned state? do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceStateResiduum(1:sizeDotState,s,g,i,e) = & - - 0.5_pReal & - * sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - sourceState(p)%p(s)%state (1:sizeDotState,c) = & - sourceState(p)%p(s)%state (1:sizeDotState,c) & - + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * crystallite_subdt(g,i,e) + + sourceStateResiduum(1:sizeDotState,s,g,i,e) = 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)%state(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) !ToDo: state, partitioned state? enddo endif enddo; enddo; enddo @@ -1845,55 +1830,51 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & call update_stress(1.0_pReal) call update_dotState(1.0_pReal) - relPlasticStateResiduum = 0.0_pReal - relSourceStateResiduum = 0.0_pReal + relPlasticStateResiduum = 0.0_pReal + relSourceStateResiduum = 0.0_pReal - !$OMP PARALLEL DO PRIVATE(sizeDotState,converged,p,c,u) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,homogenization_Ngrains(mesh_element(3,e)) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,u) + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then p = phaseAt(g,i,e); c = phasememberAt(g,i,e) + sizeDotState = plasticState(p)%sizeDotState ! --- contribution of heun step to absolute residui --- - sizeDotState = plasticState(p)%sizeDotState - plasticStateResiduum(1:sizeDotState,g,i,e) = & - plasticStateResiduum(1:sizeDotState,g,i,e) & - + 0.5_pReal * plasticState(p)%dotState(:,c) & - * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state + + plasticStateResiduum(1:sizeDotState,g,i,e) = plasticStateResiduum(1:sizeDotState,g,i,e) & + + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e) - converged = all(abs(relPlasticStateResiduum(1:sizeDotState,g,i,e)) < & + crystallite_converged(g,i,e) = all(abs(relPlasticStateResiduum(1:sizeDotState,g,i,e)) < & rTol_crystalliteState .or. & abs(plasticStateResiduum(1:sizeDotState,g,i,e)) < & plasticState(p)%aTolState(1:sizeDotState)) forall (u = 1_pInt:sizeDotState, abs(plasticState(p)%dotState(u,c)) > 0.0_pReal) & - relPlasticStateResiduum(u,g,i,e) = & - plasticStateResiduum(u,g,i,e) / plasticState(p)%dotState(u,c) + relPlasticStateResiduum(u,g,i,e) = plasticStateResiduum(u,g,i,e) & + / plasticState(p)%dotState(u,c) do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceStateResiduum(1:sizeDotState,s,g,i,e) = & - sourceStateResiduum(1:sizeDotState,s,g,i,e) & - + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) & - * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state - + + sourceStateResiduum(1:sizeDotState,s,g,i,e) = sourceStateResiduum(1:sizeDotState,s,g,i,e) & + + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e) + forall (u = 1_pInt:sizeDotState,abs(sourceState(p)%p(s)%dotState(u,c)) > 0.0_pReal) & - relSourceStateResiduum(u,s,g,i,e) = & - sourceStateResiduum(u,s,g,i,e) / sourceState(p)%p(s)%dotState(u,c) + relSourceStateResiduum(u,s,g,i,e) = sourceStateResiduum(u,s,g,i,e) & + / sourceState(p)%p(s)%dotState(u,c) - sizeDotState = sourceState(p)%p(s)%sizeDotState - converged = converged .and. & + crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and. & all(abs(relSourceStateResiduum(1:sizeDotState,s,g,i,e)) < & rTol_crystalliteState .or. & abs(sourceStateResiduum(1:sizeDotState,s,g,i,e)) < & sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo - if (converged) crystallite_converged(g,i,e) = .true. ! ... converged per definition endif - enddo; enddo; enddo + enddo; enddo; enddo !$OMP END PARALLEL DO if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck From 6a3dac1df2ba3fe6fd0cd1bb457cdcc3b175ee72 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Jan 2019 23:45:41 +0100 Subject: [PATCH 07/31] still improving readability --- src/crystallite.f90 | 72 ++++++++++++++++++++++++--------------------- 1 file changed, 38 insertions(+), 34 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 2e68ae756..7fb3aefe6 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1758,6 +1758,8 @@ end subroutine integrateStateEuler !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- subroutine integrateStateAdaptiveEuler() + use prec, only: & + dNeq0 use numerics, only: & rTol_crystalliteState use mesh, only: & @@ -1780,22 +1782,23 @@ subroutine integrateStateAdaptiveEuler() e, & ! element index in element loop i, & ! integration point index in ip loop g, & ! grain index in grain loop - u, & ! state index p, & c, & s, & sizeDotState - ! ToDo: MD: once all constitutives use allocate state, attach these arrays to the state in case of adaptive Euler + ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler + ! ToDo: MD: rel residuu don't have to be pointwise + real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - plasticStateResiduum, & ! residuum from evolution in micrstructure - relPlasticStateResiduum ! relative residuum from evolution in microstructure + residuum_plastic, & + residuum_plastic_rel real(pReal), dimension(constitutive_source_maxSizeDotState,& maxval(phase_Nsources), & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - sourceStateResiduum, & ! residuum from evolution in micrstructure - relSourceStateResiduum ! relative residuum from evolution in microstructure + residuum_source_rel, & + residuum_source !-------------------------------------------------------------------------------------------------- ! contribution to state and relative residui and from Euler integration @@ -1809,15 +1812,15 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & p = phaseAt(g,i,e); c = phasememberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - plasticStateResiduum(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) & - * (- 0.5_pReal * crystallite_subdt(g,i,e)) + residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) & + * (- 0.5_pReal * crystallite_subdt(g,i,e)) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) !ToDo: state, partitioned state? do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceStateResiduum(1:sizeDotState,s,g,i,e) = sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * (- 0.5_pReal * crystallite_subdt(g,i,e)) + residuum_source(1:sizeDotState,s,g,i,e) = 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)%state(1:sizeDotState,c) & + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) !ToDo: state, partitioned state? enddo @@ -1829,12 +1832,8 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & call update_dependentState call update_stress(1.0_pReal) call update_dotState(1.0_pReal) - - relPlasticStateResiduum = 0.0_pReal - relSourceStateResiduum = 0.0_pReal - - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,u) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) @@ -1844,33 +1843,38 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & ! --- contribution of heun step to absolute residui --- - plasticStateResiduum(1:sizeDotState,g,i,e) = plasticStateResiduum(1:sizeDotState,g,i,e) & - + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e) - - crystallite_converged(g,i,e) = all(abs(relPlasticStateResiduum(1:sizeDotState,g,i,e)) < & - rTol_crystalliteState .or. & - abs(plasticStateResiduum(1:sizeDotState,g,i,e)) < & - plasticState(p)%aTolState(1:sizeDotState)) - - forall (u = 1_pInt:sizeDotState, abs(plasticState(p)%dotState(u,c)) > 0.0_pReal) & - relPlasticStateResiduum(u,g,i,e) = plasticStateResiduum(u,g,i,e) & - / plasticState(p)%dotState(u,c) + 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) - + where(dNeq0(plasticState(p)%dotState(1:sizeDotState,c))) + residuum_plastic_rel(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) & + / plasticState(p)%dotState(1:sizeDotState,c) + else where + residuum_plastic_rel(1:sizeDotState,g,i,e) = 0.0_pReal + end where + + crystallite_converged(g,i,e) = all(abs(residuum_plastic_rel(1:sizeDotState,g,i,e)) < & + rTol_crystalliteState .or. & + abs(residuum_plastic(1:sizeDotState,g,i,e)) < & + plasticState(p)%aTolState(1:sizeDotState)) + do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceStateResiduum(1:sizeDotState,s,g,i,e) = sourceStateResiduum(1:sizeDotState,s,g,i,e) & - + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e) + 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) - forall (u = 1_pInt:sizeDotState,abs(sourceState(p)%p(s)%dotState(u,c)) > 0.0_pReal) & - relSourceStateResiduum(u,s,g,i,e) = sourceStateResiduum(u,s,g,i,e) & - / sourceState(p)%p(s)%dotState(u,c) + where(dNeq0(sourceState(p)%p(s)%dotState(1:sizeDotState,c))) + residuum_source_rel(1:sizeDotState,s,g,i,e) = residuum_source(1:sizeDotState,s,g,i,e) & + / sourceState(p)%p(s)%dotState(1:sizeDotState,c) + else where + residuum_source_rel(1:SizeDotState,s,g,i,e) = 0.0_pReal + end where crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and. & - all(abs(relSourceStateResiduum(1:sizeDotState,s,g,i,e)) < & + all(abs(residuum_source_rel(1:sizeDotState,s,g,i,e)) < & rTol_crystalliteState .or. & - abs(sourceStateResiduum(1:sizeDotState,s,g,i,e)) < & + abs(residuum_source(1:sizeDotState,s,g,i,e)) < & sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo endif From 1a66f976b7e5d7f2edf6493562e92e10ea8f10d1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 00:01:40 +0100 Subject: [PATCH 08/31] common variable name --- src/crystallite.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 7fb3aefe6..b47f3334f 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1919,7 +1919,7 @@ subroutine integrateStateRK4() p, & ! phase loop c, & n, & - mySource + s integer(pInt), dimension(2) :: eIter ! bounds for element iteration integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration gIter ! bounds for grain iteration @@ -1938,8 +1938,8 @@ subroutine integrateStateRK4() if (.not. singleRun) then do p = 1_pInt, material_Nphase plasticState(p)%RK4dotState = 0.0_pReal - do mySource = 1_pInt, phase_Nsources(p) - sourceState(p)%p(mySource)%RK4dotState = 0.0_pReal + do s = 1_pInt, phase_Nsources(p) + sourceState(p)%p(s)%RK4dotState = 0.0_pReal enddo enddo else @@ -1947,8 +1947,8 @@ subroutine integrateStateRK4() i = iIter(1,e) do g = gIter(1,e), gIter(2,e) plasticState(phaseAt(g,i,e))%RK4dotState(:,phasememberAt(g,i,e)) = 0.0_pReal - do mySource = 1_pInt, phase_Nsources(phaseAt(g,i,e)) - sourceState(phaseAt(g,i,e))%p(mySource)%RK4dotState(:,phasememberAt(g,i,e)) = 0.0_pReal + do s = 1_pInt, phase_Nsources(phaseAt(g,i,e)) + sourceState(phaseAt(g,i,e))%p(s)%RK4dotState(:,phasememberAt(g,i,e)) = 0.0_pReal enddo enddo endif @@ -1967,13 +1967,13 @@ subroutine integrateStateRK4() do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e)) then - p = phaseAt(g,i,e) - c = phasememberAt(g,i,e) + p = phaseAt(g,i,e); c = phasememberAt(g,i,e) + plasticState(p)%RK4dotState(:,c) = plasticState(p)%RK4dotState(:,c) & + weight(n)*plasticState(p)%dotState(:,c) - do mySource = 1_pInt, phase_Nsources(p) - sourceState(p)%p(mySource)%RK4dotState(:,c) = sourceState(p)%p(mySource)%RK4dotState(:,c) & - + weight(n)*sourceState(p)%p(mySource)%dotState(:,c) + do s = 1_pInt, phase_Nsources(p) + sourceState(p)%p(s)%RK4dotState(:,c) = sourceState(p)%p(s)%RK4dotState(:,c) & + + weight(n)*sourceState(p)%p(s)%dotState(:,c) enddo endif enddo; enddo; enddo From a09036ff4824531e4faebf787752bca6b60fdfbb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 00:11:10 +0100 Subject: [PATCH 09/31] on-the-fly initialization --- src/crystallite.f90 | 79 ++++++++++----------------------------------- 1 file changed, 17 insertions(+), 62 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index b47f3334f..7767eb6f3 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1893,19 +1893,13 @@ subroutine integrateStateRK4() use, intrinsic :: & IEEE_arithmetic use mesh, only: & - mesh_element, & - mesh_NcpElems + mesh_element use material, only: & homogenization_Ngrains, & plasticState, & sourceState, & phase_Nsources, & phaseAt, phasememberAt - use config, only: & - material_Nphase - use constitutive, only: & - constitutive_collectDotState, & - constitutive_microstructure implicit none real(pReal), dimension(4), parameter :: & @@ -1920,65 +1914,28 @@ subroutine integrateStateRK4() c, & n, & s - integer(pInt), dimension(2) :: eIter ! bounds for element iteration - integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration - gIter ! bounds for grain iteration - logical :: singleRun ! flag indicating computation for single (g,i,e) triple - - eIter = FEsolving_execElem(1:2) - do e = eIter(1),eIter(2) - iIter(1:2,e) = FEsolving_execIP(1:2,e) - gIter(1:2,e) = [ 1_pInt,homogenization_Ngrains(mesh_element(3,e))] - enddo - - singleRun = (eIter(1) == eIter(2) .and. iIter(1,eIter(1)) == iIter(2,eIter(2))) - -!-------------------------------------------------------------------------------------------------- -! initialize dotState - if (.not. singleRun) then - do p = 1_pInt, material_Nphase - plasticState(p)%RK4dotState = 0.0_pReal - do s = 1_pInt, phase_Nsources(p) - sourceState(p)%p(s)%RK4dotState = 0.0_pReal - enddo - enddo - else - e = eIter(1) - i = iIter(1,e) - do g = gIter(1,e), gIter(2,e) - plasticState(phaseAt(g,i,e))%RK4dotState(:,phasememberAt(g,i,e)) = 0.0_pReal - do s = 1_pInt, phase_Nsources(phaseAt(g,i,e)) - sourceState(phaseAt(g,i,e))%p(s)%RK4dotState(:,phasememberAt(g,i,e)) = 0.0_pReal - enddo - enddo - endif call update_dotState(1.0_pReal) -!-------------------------------------------------------------------------------------------------- -! --- SECOND TO FOURTH RUNGE KUTTA STEP PLUS FINAL INTEGRATION --- do n = 1_pInt,4_pInt - ! --- state update --- - !$OMP PARALLEL - !$OMP DO PRIVATE(p,c) + !$OMP PARALLEL DO PRIVATE(p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e); c = phasememberAt(g,i,e) - plasticState(p)%RK4dotState(:,c) = plasticState(p)%RK4dotState(:,c) & - + weight(n)*plasticState(p)%dotState(:,c) + plasticState(p)%RK4dotState(:,c) = WEIGHT(n)*plasticState(p)%dotState(:,c) & + + merge(plasticState(p)%RK4dotState(:,c),0.0_pReal,n>1_pInt) do s = 1_pInt, phase_Nsources(p) - sourceState(p)%p(s)%RK4dotState(:,c) = sourceState(p)%p(s)%RK4dotState(:,c) & - + weight(n)*sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%RK4dotState(:,c) = WEIGHT(n)*sourceState(p)%p(s)%dotState(:,c) & + + merge(sourceState(p)%p(s)%RK4dotState(:,c),0.0_pReal,n>1_pInt) enddo endif enddo; enddo; enddo - !$OMP ENDDO - !$OMP END PARALLEL + !$OMP END PARALLEL DO call update_state(TIMESTEPFRACTION(n)) call update_deltaState @@ -1988,9 +1945,8 @@ subroutine integrateStateRK4() ! --- dot state and RK dot state--- first3steps: if (n < 4) then - call update_dotState(timeStepFraction(n)) + call update_dotState(TIMESTEPFRACTION(n)) endif first3steps - enddo @@ -2458,9 +2414,8 @@ subroutine update_deltaState i, & !< integration point index in ip loop g, & !< grain index in grain loop p, & - mySize, & + mySize, & myOffset, & - mySource, & c, & s logical :: & @@ -2469,7 +2424,7 @@ subroutine update_deltaState nonlocalStop = .false. - !$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,mySource,NaN) + !$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,NaN) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) @@ -2489,15 +2444,15 @@ subroutine update_deltaState plasticState(p)%state(myOffset + 1_pInt: myOffset + mySize,c) = & plasticState(p)%state(myOffset + 1_pInt: myOffset + mySize,c) + & plasticState(p)%deltaState(1:mySize,c) - do mySource = 1_pInt, phase_Nsources(p) - myOffset = sourceState(p)%p(mySource)%offsetDeltaState - mySize = sourceState(p)%p(mySource)%sizeDeltaState - NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c))) + do s = 1_pInt, phase_Nsources(p) + myOffset = sourceState(p)%p(s)%offsetDeltaState + mySize = sourceState(p)%p(s)%sizeDeltaState + NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(s)%deltaState(1:mySize,c))) if (.not. NaN) then - sourceState(p)%p(mySource)%state(myOffset + 1_pInt:myOffset +mySize,c) = & - sourceState(p)%p(mySource)%state(myOffset + 1_pInt:myOffset +mySize,c) + & - sourceState(p)%p(mySource)%deltaState(1:mySize,c) + sourceState(p)%p(s)%state(myOffset + 1_pInt:myOffset +mySize,c) = & + sourceState(p)%p(s)%state(myOffset + 1_pInt:myOffset +mySize,c) + & + sourceState(p)%p(s)%deltaState(1:mySize,c) endif enddo endif From 77f1f45c231d4bcfd4dd3b6844d1cd7cbf4e1c32 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 00:17:04 +0100 Subject: [PATCH 10/31] just figured out that RK4 integrator is totally broken readable code helps ;) --- src/crystallite.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 7767eb6f3..a5f7592d6 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1888,6 +1888,7 @@ end subroutine integrateStateAdaptiveEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 4th order explicit Runge Kutta method +! ToDo: This is totally BROKEN: RK4dotState is never used!!! !-------------------------------------------------------------------------------------------------- subroutine integrateStateRK4() use, intrinsic :: & @@ -1941,7 +1942,6 @@ subroutine integrateStateRK4() call update_deltaState call update_dependentState call update_stress(TIMESTEPFRACTION(n)) - ! --- dot state and RK dot state--- first3steps: if (n < 4) then From 5908e3fd3486e1584081908ab56cfc5d1ad3a022 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 06:44:26 +0100 Subject: [PATCH 11/31] wrong tolerance selection --- src/crystallite.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index a5f7592d6..b29ede160 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1637,7 +1637,7 @@ subroutine integrateStateFPI() + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) crystallite_converged(g,i,e) = all(abs(residuum_plastic(1:sizeDotState)) & - < min(plasticState(p)%aTolState(1:sizeDotState), & + < max(plasticState(p)%aTolState(1:sizeDotState), & abs(plasticState(p)%state(1:sizeDotState,c)*rTol_crystalliteState))) @@ -1661,7 +1661,7 @@ subroutine integrateStateFPI() crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and. & all(abs(residuum_source(1:sizeDotState)) & - < min(sourceState(p)%p(s)%aTolState(1:sizeDotState), & + < max(sourceState(p)%p(s)%aTolState(1:sizeDotState), & abs(sourceState(p)%p(s)%state(1:sizeDotState,c)*rTol_crystalliteState))) enddo endif From 462b1b7c189e8370c9930736b15ffc8ed22306f9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 06:47:36 +0100 Subject: [PATCH 12/31] sorted according to importance --- src/crystallite.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index b29ede160..0150d68b0 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1723,7 +1723,7 @@ subroutine integrateStateFPI() dot_prod12 = dot_product(current - previous, previous - previous2) dot_prod22 = dot_product(previous - previous2, previous - previous2) - if (dot_prod22 > 0.0_pReal .and. (dot_prod12 < 0.0_pReal .or. dot_product(current,previous) < 0.0_pReal)) then + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) else damper = 1.0_pReal From ca7c105f363c80d49bce0fc5b9fd9add335961cf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 08:56:16 +0100 Subject: [PATCH 13/31] only one loop needed --- src/crystallite.f90 | 35 +++++++++++------------------------ 1 file changed, 11 insertions(+), 24 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0150d68b0..b0f1c1f94 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2038,46 +2038,33 @@ subroutine integrateStateRKCK45() ! --- state update --- - !$OMP PARALLEL - !$OMP DO PRIVATE(p,cc) + !$OMP PARALLEL DO PRIVATE(p,cc) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e) cc = phasememberAt(g,i,e) - plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc) ! store Runge-Kutta dotState + plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc) + plasticState(p)%dotState(:,cc) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,cc) + do mySource = 1_pInt, phase_Nsources(p) sourceState(p)%p(mySource)%RKCK45dotState(stage,:,cc) = sourceState(p)%p(mySource)%dotState(:,cc) - enddo - endif - enddo; enddo; enddo - !$OMP ENDDO - - !$OMP DO PRIVATE(p,cc,n) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,homogenization_Ngrains(mesh_element(3,e)) - if (crystallite_todo(g,i,e)) then - p = phaseAt(g,i,e) - cc = phasememberAt(g,i,e) - - plasticState(p)%dotState(:,cc) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,cc) - do mySource = 1_pInt, phase_Nsources(p) sourceState(p)%p(mySource)%dotState(:,cc) = A(1,stage) * sourceState(p)%p(mySource)%RKCK45dotState(1,:,cc) enddo + do n = 2_pInt, stage - plasticState(p)%dotState(:,cc) = & - plasticState(p)%dotState(:,cc) + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,cc) + plasticState(p)%dotState(:,cc) = plasticState(p)%dotState(:,cc) & + + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,cc) do mySource = 1_pInt, phase_Nsources(p) - sourceState(p)%p(mySource)%dotState(:,cc) = & - sourceState(p)%p(mySource)%dotState(:,cc) + A(n,stage) * sourceState(p)%p(mySource)%RKCK45dotState(n,:,cc) + sourceState(p)%p(mySource)%dotState(:,cc) = sourceState(p)%p(mySource)%dotState(:,cc) & + + A(n,stage) * sourceState(p)%p(mySource)%RKCK45dotState(n,:,cc) enddo enddo + endif enddo; enddo; enddo - !$OMP ENDDO - !$OMP END PARALLEL + !$OMP END PARALLEL DO call update_state(1.0_pReal) !MD: 1.0 correct? call update_deltaState From df6ec59f76cdfa25e69e8e80486de4d47b56787b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 09:11:12 +0100 Subject: [PATCH 14/31] use "s" for source --- src/crystallite.f90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index b0f1c1f94..de535e8c2 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2008,7 +2008,7 @@ subroutine integrateStateRKCK45() i, & ! integration point index in ip loop g, & ! grain index in grain loop stage, & ! stage index in integration stage loop - s, & ! state index + u, & ! state index n, & p, & cc, & @@ -2043,8 +2043,8 @@ subroutine integrateStateRKCK45() do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e)) then - p = phaseAt(g,i,e) - cc = phasememberAt(g,i,e) + p = phaseAt(g,i,e); cc = phasememberAt(g,i,e) + plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc) plasticState(p)%dotState(:,cc) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,cc) @@ -2134,7 +2134,7 @@ subroutine integrateStateRKCK45() !$OMP PARALLEL ! --- relative residui and state convergence --- - !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc,s) + !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc,u) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) @@ -2142,15 +2142,15 @@ subroutine integrateStateRKCK45() p = phaseAt(g,i,e) cc = phasememberAt(g,i,e) mySizePlasticDotState = plasticState(p)%sizeDotState - forall (s = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%state(s,cc)) > 0.0_pReal) & - relPlasticStateResiduum(s,g,i,e) = & - plasticStateResiduum(s,g,i,e) / plasticState(p)%state(s,cc) + forall (u = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%state(u,cc)) > 0.0_pReal) & + relPlasticStateResiduum(u,g,i,e) = & + plasticStateResiduum(u,g,i,e) / plasticState(p)%state(u,cc) do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState - forall (s = 1_pInt:mySizeSourceDotState,abs(sourceState(p)%p(mySource)%state(s,cc)) > 0.0_pReal) & - relSourceStateResiduum(s,mySource,g,i,e) = & - sourceStateResiduum(s,mySource,g,i,e) / sourceState(p)%p(mySource)%state(s,cc) + forall (u = 1_pInt:mySizeSourceDotState,abs(sourceState(p)%p(mySource)%state(u,cc)) > 0.0_pReal) & + relSourceStateResiduum(u,mySource,g,i,e) = & + sourceStateResiduum(u,mySource,g,i,e) / sourceState(p)%p(mySource)%state(u,cc) enddo crystallite_todo(g,i,e) = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & rTol_crystalliteState .or. & From 31906e3ebd70ea033e9f7fb652cd6280278f1fbd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 09:21:33 +0100 Subject: [PATCH 15/31] no need for 2 loops --- src/crystallite.f90 | 47 +++++++++++++++++---------------------------- 1 file changed, 18 insertions(+), 29 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index de535e8c2..3239f12c4 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2016,6 +2016,8 @@ subroutine integrateStateRKCK45() mySizePlasticDotState, & ! size of dot States mySizeSourceDotState + ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler + ! ToDo: MD: rel residuu don't have to be pointwise real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & @@ -2080,54 +2082,41 @@ subroutine integrateStateRKCK45() relPlasticStateResiduum = 0.0_pReal relSourceStateResiduum = 0.0_pReal - !$OMP PARALLEL - !$OMP DO PRIVATE(p,cc) + !$OMP PARALLEL DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e)) then - p = phaseAt(g,i,e) - cc = phasememberAt(g,i,e) - plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) ! store Runge-Kutta dotState - do mySource = 1_pInt, phase_Nsources(p) - sourceState(p)%p(mySource)%RKCK45dotState(6,:,cc) = sourceState(p)%p(mySource)%dotState(:,cc) ! store Runge-Kutta dotState - enddo - endif - enddo; enddo; enddo - !$OMP ENDDO - - !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,homogenization_Ngrains(mesh_element(3,e)) - if (crystallite_todo(g,i,e)) then - p = phaseAt(g,i,e) - cc = phasememberAt(g,i,e) - - ! --- absolute residuum in state --- - mySizePlasticDotState = plasticState(p)%sizeDotState + p = phaseAt(g,i,e); cc = phasememberAt(g,i,e) + + mySizePlasticDotState = plasticState(p)%sizeDotState + + plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) + plasticStateResiduum(1:mySizePlasticDotState,g,i,e) = & matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:mySizePlasticDotState,cc)),DB) & * crystallite_subdt(g,i,e) + + plasticState(p)%dotState(:,cc) = & + matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:mySizePlasticDotState,cc)), B) + do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState + + sourceState(p)%p(mySource)%RKCK45dotState(6,:,cc) = sourceState(p)%p(mySource)%dotState(:,cc) ! store Runge-Kutta dotState + sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e) = & matmul(transpose(sourceState(p)%p(mySource)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),DB) & * crystallite_subdt(g,i,e) - enddo - ! --- dot state --- - plasticState(p)%dotState(:,cc) = & - matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:mySizePlasticDotState,cc)), B) - do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState sourceState(p)%p(mySource)%dotState(:,cc) = & matmul(transpose(sourceState(p)%p(mySource)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),B) enddo + endif enddo; enddo; enddo - !$OMP ENDDO - !$OMP END PARALLEL + !$OMP END PARALLEL DO call update_state(1.0_pReal) From 46be595ea803004d09a157c440a2848ad33e7f9e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 10:28:47 +0100 Subject: [PATCH 16/31] no need to store relative residual for all points --- src/crystallite.f90 | 43 +++++++++++++++++++++---------------------- 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 3239f12c4..d24b16dbf 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1787,18 +1787,19 @@ subroutine integrateStateAdaptiveEuler() s, & sizeDotState - ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler - ! ToDo: MD: rel residuu don't have to be pointwise - -real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & + ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler + real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - residuum_plastic, & - residuum_plastic_rel + residuum_plastic real(pReal), dimension(constitutive_source_maxSizeDotState,& maxval(phase_Nsources), & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - residuum_source_rel, & residuum_source + + real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & + residuum_plastic_rel + real(pReal), dimension(constitutive_source_maxSizeDotState) :: & + residuum_source_rel !-------------------------------------------------------------------------------------------------- ! contribution to state and relative residui and from Euler integration @@ -1828,10 +1829,10 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & enddo; enddo; enddo !$OMP END PARALLEL DO - call update_deltaState - call update_dependentState - call update_stress(1.0_pReal) - call update_dotState(1.0_pReal) + call update_deltaState + call update_dependentState + call update_stress(1.0_pReal) + call update_dotState(1.0_pReal) !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -1840,20 +1841,18 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then p = phaseAt(g,i,e); c = phasememberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - - ! --- contribution of heun step to absolute residui --- 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) where(dNeq0(plasticState(p)%dotState(1:sizeDotState,c))) - residuum_plastic_rel(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) & - / plasticState(p)%dotState(1:sizeDotState,c) + residuum_plastic_rel(1:sizeDotState) = residuum_plastic(1:sizeDotState,g,i,e) & + / plasticState(p)%dotState(1:sizeDotState,c) else where - residuum_plastic_rel(1:sizeDotState,g,i,e) = 0.0_pReal + residuum_plastic_rel(1:sizeDotState) = 0.0_pReal end where - crystallite_converged(g,i,e) = all(abs(residuum_plastic_rel(1:sizeDotState,g,i,e)) < & + crystallite_converged(g,i,e) = all(abs(residuum_plastic_rel(1:sizeDotState)) < & rTol_crystalliteState .or. & abs(residuum_plastic(1:sizeDotState,g,i,e)) < & plasticState(p)%aTolState(1:sizeDotState)) @@ -1865,14 +1864,14 @@ real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e) where(dNeq0(sourceState(p)%p(s)%dotState(1:sizeDotState,c))) - residuum_source_rel(1:sizeDotState,s,g,i,e) = residuum_source(1:sizeDotState,s,g,i,e) & - / sourceState(p)%p(s)%dotState(1:sizeDotState,c) + residuum_source_rel(1:sizeDotState) = residuum_source(1:sizeDotState,s,g,i,e) & + / sourceState(p)%p(s)%dotState(1:sizeDotState,c) else where - residuum_source_rel(1:SizeDotState,s,g,i,e) = 0.0_pReal + residuum_source_rel(1:SizeDotState) = 0.0_pReal end where - crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and. & - all(abs(residuum_source_rel(1:sizeDotState,s,g,i,e)) < & + crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and. & + all(abs(residuum_source_rel(1:sizeDotState)) < & rTol_crystalliteState .or. & abs(residuum_source(1:sizeDotState,s,g,i,e)) < & sourceState(p)%p(s)%aTolState(1:sizeDotState)) From 0745d7ebc20ab6803869d88ca88cb56a8afb0aca Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 10:33:57 +0100 Subject: [PATCH 17/31] convergence flag is set only later --- src/crystallite.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d24b16dbf..053aa35eb 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1809,7 +1809,7 @@ subroutine integrateStateAdaptiveEuler() do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) - if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then + if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e); c = phasememberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState @@ -1838,7 +1838,7 @@ subroutine integrateStateAdaptiveEuler() do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) - if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then + if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e); c = phasememberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState From 72c4f2b25fae73f0c9ec665fbb1132571f80aa51 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 10:37:18 +0100 Subject: [PATCH 18/31] same names everywhere if possible --- src/crystallite.f90 | 54 ++++++++++++++++++++++----------------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 053aa35eb..80e1a7ed6 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2011,7 +2011,7 @@ subroutine integrateStateRKCK45() n, & p, & cc, & - mySource, & + s, & mySizePlasticDotState, & ! size of dot States mySizeSourceDotState @@ -2049,17 +2049,17 @@ subroutine integrateStateRKCK45() plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc) plasticState(p)%dotState(:,cc) = A(1,stage) * plasticState(p)%RKCK45dotState(1,:,cc) - do mySource = 1_pInt, phase_Nsources(p) - sourceState(p)%p(mySource)%RKCK45dotState(stage,:,cc) = sourceState(p)%p(mySource)%dotState(:,cc) - sourceState(p)%p(mySource)%dotState(:,cc) = A(1,stage) * sourceState(p)%p(mySource)%RKCK45dotState(1,:,cc) + do s = 1_pInt, phase_Nsources(p) + sourceState(p)%p(s)%RKCK45dotState(stage,:,cc) = sourceState(p)%p(s)%dotState(:,cc) + sourceState(p)%p(s)%dotState(:,cc) = A(1,stage) * sourceState(p)%p(s)%RKCK45dotState(1,:,cc) enddo do n = 2_pInt, stage plasticState(p)%dotState(:,cc) = plasticState(p)%dotState(:,cc) & + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,cc) - do mySource = 1_pInt, phase_Nsources(p) - sourceState(p)%p(mySource)%dotState(:,cc) = sourceState(p)%p(mySource)%dotState(:,cc) & - + A(n,stage) * sourceState(p)%p(mySource)%RKCK45dotState(n,:,cc) + do s = 1_pInt, phase_Nsources(p) + sourceState(p)%p(s)%dotState(:,cc) = sourceState(p)%p(s)%dotState(:,cc) & + + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,cc) enddo enddo @@ -2099,18 +2099,18 @@ subroutine integrateStateRKCK45() plasticState(p)%dotState(:,cc) = & matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:mySizePlasticDotState,cc)), B) - do mySource = 1_pInt, phase_Nsources(p) - mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState + do s = 1_pInt, phase_Nsources(p) + mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(mySource)%RKCK45dotState(6,:,cc) = sourceState(p)%p(mySource)%dotState(:,cc) ! store Runge-Kutta dotState + sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc) ! store Runge-Kutta dotState - sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e) = & - matmul(transpose(sourceState(p)%p(mySource)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),DB) & + sourceStateResiduum(1:mySizeSourceDotState,s,g,i,e) = & + matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),DB) & * crystallite_subdt(g,i,e) - mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState - sourceState(p)%p(mySource)%dotState(:,cc) = & - matmul(transpose(sourceState(p)%p(mySource)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),B) + mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%dotState(:,cc) = & + matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),B) enddo endif @@ -2127,30 +2127,30 @@ subroutine integrateStateRKCK45() do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e)) then - p = phaseAt(g,i,e) - cc = phasememberAt(g,i,e) + p = phaseAt(g,i,e); cc = phasememberAt(g,i,e) + mySizePlasticDotState = plasticState(p)%sizeDotState forall (u = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%state(u,cc)) > 0.0_pReal) & relPlasticStateResiduum(u,g,i,e) = & plasticStateResiduum(u,g,i,e) / plasticState(p)%state(u,cc) - do mySource = 1_pInt, phase_Nsources(p) - mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState - forall (u = 1_pInt:mySizeSourceDotState,abs(sourceState(p)%p(mySource)%state(u,cc)) > 0.0_pReal) & - relSourceStateResiduum(u,mySource,g,i,e) = & - sourceStateResiduum(u,mySource,g,i,e) / sourceState(p)%p(mySource)%state(u,cc) + do s = 1_pInt, phase_Nsources(p) + mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState + forall (u = 1_pInt:mySizeSourceDotState,abs(sourceState(p)%p(s)%state(u,cc)) > 0.0_pReal) & + relSourceStateResiduum(u,s,g,i,e) = & + sourceStateResiduum(u,s,g,i,e) / sourceState(p)%p(s)%state(u,cc) enddo crystallite_todo(g,i,e) = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & rTol_crystalliteState .or. & abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & plasticState(p)%aTolState(1:mySizePlasticDotState)) - do mySource = 1_pInt, phase_Nsources(p) - mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState + do s = 1_pInt, phase_Nsources(p) + mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & - all(abs(relSourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e)) < & + all(abs(relSourceStateResiduum(1:mySizeSourceDotState,s,g,i,e)) < & rTol_crystalliteState .or. & - abs(sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e)) < & - sourceState(p)%p(mySource)%aTolState(1:mySizeSourceDotState)) + abs(sourceStateResiduum(1:mySizeSourceDotState,s,g,i,e)) < & + sourceState(p)%p(s)%aTolState(1:mySizeSourceDotState)) enddo endif enddo; enddo; enddo From 0876787e3c56bd201214d599a74c1ce1c11ef9ce Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 10:46:53 +0100 Subject: [PATCH 19/31] avoid loops --- src/crystallite.f90 | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 80e1a7ed6..cc261726f 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1976,10 +1976,8 @@ subroutine integrateStateRKCK45() phaseAt, phasememberAt, & homogenization_maxNgrains use constitutive, only: & - constitutive_collectDotState, & constitutive_plasticity_maxSizeDotState, & - constitutive_source_maxSizeDotState, & - constitutive_microstructure + constitutive_source_maxSizeDotState implicit none real(pReal), dimension(5,5), parameter :: & @@ -2059,7 +2057,7 @@ subroutine integrateStateRKCK45() + A(n,stage) * plasticState(p)%RKCK45dotState(n,:,cc) do s = 1_pInt, phase_Nsources(p) sourceState(p)%p(s)%dotState(:,cc) = sourceState(p)%p(s)%dotState(:,cc) & - + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,cc) + + A(n,stage) * sourceState(p)%p(s)%RKCK45dotState(n,:,cc) enddo enddo @@ -2088,7 +2086,7 @@ subroutine integrateStateRKCK45() if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e); cc = phasememberAt(g,i,e) - mySizePlasticDotState = plasticState(p)%sizeDotState + mySizePlasticDotState = plasticState(p)%sizeDotState plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) @@ -2133,18 +2131,18 @@ subroutine integrateStateRKCK45() forall (u = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%state(u,cc)) > 0.0_pReal) & relPlasticStateResiduum(u,g,i,e) = & plasticStateResiduum(u,g,i,e) / plasticState(p)%state(u,cc) + + crystallite_todo(g,i,e) = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & + rTol_crystalliteState .or. & + abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & + plasticState(p)%aTolState(1:mySizePlasticDotState)) do s = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState forall (u = 1_pInt:mySizeSourceDotState,abs(sourceState(p)%p(s)%state(u,cc)) > 0.0_pReal) & relSourceStateResiduum(u,s,g,i,e) = & sourceStateResiduum(u,s,g,i,e) / sourceState(p)%p(s)%state(u,cc) - enddo - crystallite_todo(g,i,e) = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & - rTol_crystalliteState .or. & - abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & - plasticState(p)%aTolState(1:mySizePlasticDotState)) - do s = 1_pInt, phase_Nsources(p) + mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & all(abs(relSourceStateResiduum(1:mySizeSourceDotState,s,g,i,e)) < & From 4ec0fd70a2574e2caa84497fea165d53fcffb608 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 10:48:59 +0100 Subject: [PATCH 20/31] only one variable needed --- src/crystallite.f90 | 47 ++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index cc261726f..81b730aad 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2010,8 +2010,7 @@ subroutine integrateStateRKCK45() p, & cc, & s, & - mySizePlasticDotState, & ! size of dot States - mySizeSourceDotState + sizeDotState ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler ! ToDo: MD: rel residuu don't have to be pointwise @@ -2079,36 +2078,36 @@ subroutine integrateStateRKCK45() relPlasticStateResiduum = 0.0_pReal relSourceStateResiduum = 0.0_pReal - !$OMP PARALLEL DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e); cc = phasememberAt(g,i,e) - mySizePlasticDotState = plasticState(p)%sizeDotState + sizeDotState = plasticState(p)%sizeDotState plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) - plasticStateResiduum(1:mySizePlasticDotState,g,i,e) = & - matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:mySizePlasticDotState,cc)),DB) & + plasticStateResiduum(1:sizeDotState,g,i,e) = & + matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & * crystallite_subdt(g,i,e) plasticState(p)%dotState(:,cc) = & - matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:mySizePlasticDotState,cc)), B) + matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)), B) do s = 1_pInt, phase_Nsources(p) - mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState + sizeDotState = sourceState(p)%p(s)%sizeDotState sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc) ! store Runge-Kutta dotState - sourceStateResiduum(1:mySizeSourceDotState,s,g,i,e) = & - matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),DB) & + sourceStateResiduum(1:sizeDotState,s,g,i,e) = & + matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & * crystallite_subdt(g,i,e) - mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState + sizeDotState = sourceState(p)%p(s)%sizeDotState sourceState(p)%p(s)%dotState(:,cc) = & - matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),B) + matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),B) enddo endif @@ -2120,35 +2119,35 @@ subroutine integrateStateRKCK45() !$OMP PARALLEL ! --- relative residui and state convergence --- - !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc,u) + !$OMP DO PRIVATE(sizeDotState,p,cc,u) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) if (crystallite_todo(g,i,e)) then p = phaseAt(g,i,e); cc = phasememberAt(g,i,e) - mySizePlasticDotState = plasticState(p)%sizeDotState - forall (u = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%state(u,cc)) > 0.0_pReal) & + sizeDotState = plasticState(p)%sizeDotState + forall (u = 1_pInt:sizeDotState, abs(plasticState(p)%state(u,cc)) > 0.0_pReal) & relPlasticStateResiduum(u,g,i,e) = & plasticStateResiduum(u,g,i,e) / plasticState(p)%state(u,cc) - crystallite_todo(g,i,e) = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & + crystallite_todo(g,i,e) = all(abs(relPlasticStateResiduum(1:sizeDotState,g,i,e)) < & rTol_crystalliteState .or. & - abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & - plasticState(p)%aTolState(1:mySizePlasticDotState)) + abs(plasticStateResiduum(1:sizeDotState,g,i,e)) < & + plasticState(p)%aTolState(1:sizeDotState)) do s = 1_pInt, phase_Nsources(p) - mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState - forall (u = 1_pInt:mySizeSourceDotState,abs(sourceState(p)%p(s)%state(u,cc)) > 0.0_pReal) & + sizeDotState = sourceState(p)%p(s)%sizeDotState + forall (u = 1_pInt:sizeDotState,abs(sourceState(p)%p(s)%state(u,cc)) > 0.0_pReal) & relSourceStateResiduum(u,s,g,i,e) = & sourceStateResiduum(u,s,g,i,e) / sourceState(p)%p(s)%state(u,cc) - mySizeSourceDotState = sourceState(p)%p(s)%sizeDotState + sizeDotState = sourceState(p)%p(s)%sizeDotState crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & - all(abs(relSourceStateResiduum(1:mySizeSourceDotState,s,g,i,e)) < & + all(abs(relSourceStateResiduum(1:sizeDotState,s,g,i,e)) < & rTol_crystalliteState .or. & - abs(sourceStateResiduum(1:mySizeSourceDotState,s,g,i,e)) < & - sourceState(p)%p(s)%aTolState(1:mySizeSourceDotState)) + abs(sourceStateResiduum(1:sizeDotState,s,g,i,e)) < & + sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo endif enddo; enddo; enddo From fd069a96cdc68941a44de7aca9083a2b03b5d5d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 10:51:24 +0100 Subject: [PATCH 21/31] unifying name --- src/crystallite.f90 | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 81b730aad..f9f469a5d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2017,13 +2017,13 @@ subroutine integrateStateRKCK45() real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - plasticStateResiduum, & ! residuum from evolution in microstructure - relPlasticStateResiduum ! relative residuum from evolution in microstructure + residuum_plastic, & ! residuum from evolution in microstructure + residuum_plastic_rel ! relative residuum from evolution in microstructure real(pReal), dimension(constitutive_source_maxSizeDotState, & maxval(phase_Nsources), & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - sourceStateResiduum, & ! residuum from evolution in microstructure - relSourceStateResiduum ! relative residuum from evolution in microstructure + residuum_source, & ! residuum from evolution in microstructure + residuum_source_rel ! relative residuum from evolution in microstructure @@ -2076,8 +2076,8 @@ subroutine integrateStateRKCK45() !-------------------------------------------------------------------------------------------------- ! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE --- - relPlasticStateResiduum = 0.0_pReal - relSourceStateResiduum = 0.0_pReal + residuum_plastic_rel = 0.0_pReal + residuum_source_rel = 0.0_pReal !$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) @@ -2089,7 +2089,7 @@ subroutine integrateStateRKCK45() plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) - plasticStateResiduum(1:sizeDotState,g,i,e) = & + residuum_plastic(1:sizeDotState,g,i,e) = & matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & * crystallite_subdt(g,i,e) @@ -2101,7 +2101,7 @@ subroutine integrateStateRKCK45() sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc) ! store Runge-Kutta dotState - sourceStateResiduum(1:sizeDotState,s,g,i,e) = & + residuum_source(1:sizeDotState,s,g,i,e) = & matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & * crystallite_subdt(g,i,e) @@ -2128,25 +2128,25 @@ subroutine integrateStateRKCK45() sizeDotState = plasticState(p)%sizeDotState forall (u = 1_pInt:sizeDotState, abs(plasticState(p)%state(u,cc)) > 0.0_pReal) & - relPlasticStateResiduum(u,g,i,e) = & - plasticStateResiduum(u,g,i,e) / plasticState(p)%state(u,cc) + residuum_plastic_rel(u,g,i,e) = & + residuum_plastic(u,g,i,e) / plasticState(p)%state(u,cc) - crystallite_todo(g,i,e) = all(abs(relPlasticStateResiduum(1:sizeDotState,g,i,e)) < & + crystallite_todo(g,i,e) = all(abs(residuum_plastic_rel(1:sizeDotState,g,i,e)) < & rTol_crystalliteState .or. & - abs(plasticStateResiduum(1:sizeDotState,g,i,e)) < & + abs(residuum_plastic(1:sizeDotState,g,i,e)) < & plasticState(p)%aTolState(1:sizeDotState)) do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState forall (u = 1_pInt:sizeDotState,abs(sourceState(p)%p(s)%state(u,cc)) > 0.0_pReal) & - relSourceStateResiduum(u,s,g,i,e) = & - sourceStateResiduum(u,s,g,i,e) / sourceState(p)%p(s)%state(u,cc) + residuum_source_rel(u,s,g,i,e) = & + residuum_source(u,s,g,i,e) / sourceState(p)%p(s)%state(u,cc) sizeDotState = sourceState(p)%p(s)%sizeDotState crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & - all(abs(relSourceStateResiduum(1:sizeDotState,s,g,i,e)) < & + all(abs(residuum_source_rel(1:sizeDotState,s,g,i,e)) < & rTol_crystalliteState .or. & - abs(sourceStateResiduum(1:sizeDotState,s,g,i,e)) < & + abs(residuum_source(1:sizeDotState,s,g,i,e)) < & sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo endif From 3dd21177a0464faf46aca285cb7d3f8ca7325743 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 11:04:49 +0100 Subject: [PATCH 22/31] no need to store relative residual pointwise --- src/crystallite.f90 | 76 ++++++++++++++++++++++----------------------- 1 file changed, 38 insertions(+), 38 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index f9f469a5d..0a190e364 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1890,8 +1890,6 @@ end subroutine integrateStateAdaptiveEuler ! ToDo: This is totally BROKEN: RK4dotState is never used!!! !-------------------------------------------------------------------------------------------------- subroutine integrateStateRK4() - use, intrinsic :: & - IEEE_arithmetic use mesh, only: & mesh_element use material, only: & @@ -1960,8 +1958,8 @@ end subroutine integrateStateRK4 !> adaptive step size (use 5th order solution to advance = "local extrapolation") !-------------------------------------------------------------------------------------------------- subroutine integrateStateRKCK45() - use, intrinsic :: & - IEEE_arithmetic + use prec, only: & + dNeq0 use numerics, only: & rTol_crystalliteState use mesh, only: & @@ -2005,26 +2003,25 @@ subroutine integrateStateRKCK45() i, & ! integration point index in ip loop g, & ! grain index in grain loop stage, & ! stage index in integration stage loop - u, & ! state index n, & p, & cc, & s, & sizeDotState - ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of adaptive Euler - ! ToDo: MD: rel residuu don't have to be pointwise + ! ToDo: MD: once all constitutives use allocate state, attach residuum arrays to the state in case of RKCK45 real(pReal), dimension(constitutive_plasticity_maxSizeDotState, & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - residuum_plastic, & ! residuum from evolution in microstructure - residuum_plastic_rel ! relative residuum from evolution in microstructure + residuum_plastic ! relative residuum from evolution in microstructure real(pReal), dimension(constitutive_source_maxSizeDotState, & maxval(phase_Nsources), & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & - residuum_source, & ! residuum from evolution in microstructure - residuum_source_rel ! relative residuum from evolution in microstructure - + residuum_source ! relative residuum from evolution in microstructure + real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & + residuum_plastic_rel + real(pReal), dimension(constitutive_source_maxSizeDotState) :: & + residuum_source_rel call update_dotState(1.0_pReal) @@ -2076,8 +2073,6 @@ subroutine integrateStateRKCK45() !-------------------------------------------------------------------------------------------------- ! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE --- - residuum_plastic_rel = 0.0_pReal - residuum_source_rel = 0.0_pReal !$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) @@ -2116,43 +2111,48 @@ subroutine integrateStateRKCK45() call update_state(1.0_pReal) -!$OMP PARALLEL ! --- relative residui and state convergence --- - !$OMP DO PRIVATE(sizeDotState,p,cc,u) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do g = 1,homogenization_Ngrains(mesh_element(3,e)) - if (crystallite_todo(g,i,e)) then - p = phaseAt(g,i,e); cc = phasememberAt(g,i,e) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc) + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + do g = 1,homogenization_Ngrains(mesh_element(3,e)) + if (crystallite_todo(g,i,e)) then + p = phaseAt(g,i,e); cc = phasememberAt(g,i,e) - sizeDotState = plasticState(p)%sizeDotState - forall (u = 1_pInt:sizeDotState, abs(plasticState(p)%state(u,cc)) > 0.0_pReal) & - residuum_plastic_rel(u,g,i,e) = & - residuum_plastic(u,g,i,e) / plasticState(p)%state(u,cc) + sizeDotState = plasticState(p)%sizeDotState + where(dNeq0(plasticState(p)%dotState(1:sizeDotState,cc))) + residuum_plastic_rel(1:sizeDotState) = residuum_plastic(1:sizeDotState,g,i,e) & + / plasticState(p)%state(1:sizeDotState,cc) + else where + residuum_plastic_rel(1:sizeDotState) = 0.0_pReal + end where + - crystallite_todo(g,i,e) = all(abs(residuum_plastic_rel(1:sizeDotState,g,i,e)) < & - rTol_crystalliteState .or. & - abs(residuum_plastic(1:sizeDotState,g,i,e)) < & - plasticState(p)%aTolState(1:sizeDotState)) + crystallite_todo(g,i,e) = all(abs(residuum_plastic_rel(1:sizeDotState)) < & + rTol_crystalliteState .or. & + abs(residuum_plastic(1:sizeDotState,g,i,e)) < & + plasticState(p)%aTolState(1:sizeDotState)) - do s = 1_pInt, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - forall (u = 1_pInt:sizeDotState,abs(sourceState(p)%p(s)%state(u,cc)) > 0.0_pReal) & - residuum_source_rel(u,s,g,i,e) = & - residuum_source(u,s,g,i,e) / sourceState(p)%p(s)%state(u,cc) + do s = 1_pInt, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + + where(dNeq0(sourceState(p)%p(s)%dotState(1:sizeDotState,cc))) + residuum_source_rel(1:sizeDotState) = residuum_source(1:sizeDotState,s,g,i,e) & + / sourceState(p)%p(s)%state(1:sizeDotState,cc) + else where + residuum_source_rel(1:SizeDotState) = 0.0_pReal + end where - sizeDotState = sourceState(p)%p(s)%sizeDotState crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & - all(abs(residuum_source_rel(1:sizeDotState,s,g,i,e)) < & + all(abs(residuum_source_rel(1:sizeDotState)) < & rTol_crystalliteState .or. & abs(residuum_source(1:sizeDotState,s,g,i,e)) < & sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo endif enddo; enddo; enddo - !$OMP ENDDO -!$OMP END PARALLEL + !$OMP END PARALLEL DO call update_deltaState call update_dependentState From 39e766bba006e52742a952fbf523e413fa02750d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 12:36:02 +0100 Subject: [PATCH 23/31] improved readability --- src/crystallite.f90 | 139 +++++++++++++++++++++++--------------------- 1 file changed, 73 insertions(+), 66 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0a190e364..210bf8198 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1758,10 +1758,6 @@ end subroutine integrateStateEuler !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- subroutine integrateStateAdaptiveEuler() - use prec, only: & - dNeq0 - use numerics, only: & - rTol_crystalliteState use mesh, only: & mesh_element, & mesh_NcpElems, & @@ -1795,11 +1791,6 @@ subroutine integrateStateAdaptiveEuler() maxval(phase_Nsources), & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & residuum_source - - real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & - residuum_plastic_rel - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & - residuum_source_rel !-------------------------------------------------------------------------------------------------- ! contribution to state and relative residui and from Euler integration @@ -1845,42 +1836,55 @@ subroutine integrateStateAdaptiveEuler() 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) - where(dNeq0(plasticState(p)%dotState(1:sizeDotState,c))) - residuum_plastic_rel(1:sizeDotState) = residuum_plastic(1:sizeDotState,g,i,e) & - / plasticState(p)%dotState(1:sizeDotState,c) - else where - residuum_plastic_rel(1:sizeDotState) = 0.0_pReal - end where - - crystallite_converged(g,i,e) = all(abs(residuum_plastic_rel(1:sizeDotState)) < & - rTol_crystalliteState .or. & - abs(residuum_plastic(1:sizeDotState,g,i,e)) < & - plasticState(p)%aTolState(1:sizeDotState)) + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & + plasticState(p)%dotState(1:sizeDotState,c), & + plasticState(p)%aTolState(1:sizeDotState)) do s = 1_pInt, 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) - - where(dNeq0(sourceState(p)%p(s)%dotState(1:sizeDotState,c))) - residuum_source_rel(1:sizeDotState) = residuum_source(1:sizeDotState,s,g,i,e) & - / sourceState(p)%p(s)%dotState(1:sizeDotState,c) - else where - residuum_source_rel(1:SizeDotState) = 0.0_pReal - end where - crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and. & - all(abs(residuum_source_rel(1:sizeDotState)) < & - rTol_crystalliteState .or. & - abs(residuum_source(1:sizeDotState,s,g,i,e)) < & - sourceState(p)%p(s)%aTolState(1:sizeDotState)) - enddo + crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and.& + converged(residuum_source(1:sizeDotState,s,g,i,e), & + sourceState(p)%p(s)%dotState(1:sizeDotState,c), & + sourceState(p)%p(s)%aTolState(1:sizeDotState)) + enddo + endif enddo; enddo; enddo !$OMP END PARALLEL DO if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief determines whether a point is converged + !-------------------------------------------------------------------------------------------------- + logical pure function converged(residuum,dotState,absoluteTolerance) + use prec, only: & + dNeq0 + use numerics, only: & + rTol_crystalliteState + + implicit none + real(pReal), dimension(:), intent(in) ::& + residuum, dotState, absoluteTolerance + real(pReal), dimension(size(residuum,1)) ::& + residuum_rel + + where(dNeq0(dotState)) + residuum_rel = residuum/dotState + else where + residuum_rel = 0.0_pReal + end where + + converged = all(abs(residuum_rel) < rTol_crystalliteState .or. & + abs(residuum) < absoluteTolerance) + + end function converged end subroutine integrateStateAdaptiveEuler @@ -1958,10 +1962,6 @@ end subroutine integrateStateRK4 !> adaptive step size (use 5th order solution to advance = "local extrapolation") !-------------------------------------------------------------------------------------------------- subroutine integrateStateRKCK45() - use prec, only: & - dNeq0 - use numerics, only: & - rTol_crystalliteState use mesh, only: & mesh_element, & mesh_NcpElems, & @@ -2018,15 +2018,10 @@ subroutine integrateStateRKCK45() maxval(phase_Nsources), & homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & residuum_source ! relative residuum from evolution in microstructure - real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: & - residuum_plastic_rel - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & - residuum_source_rel call update_dotState(1.0_pReal) - ! --- SECOND TO SIXTH RUNGE KUTTA STEP --- do stage = 1_pInt,5_pInt @@ -2121,34 +2116,18 @@ subroutine integrateStateRKCK45() p = phaseAt(g,i,e); cc = phasememberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - where(dNeq0(plasticState(p)%dotState(1:sizeDotState,cc))) - residuum_plastic_rel(1:sizeDotState) = residuum_plastic(1:sizeDotState,g,i,e) & - / plasticState(p)%state(1:sizeDotState,cc) - else where - residuum_plastic_rel(1:sizeDotState) = 0.0_pReal - end where - - - crystallite_todo(g,i,e) = all(abs(residuum_plastic_rel(1:sizeDotState)) < & - rTol_crystalliteState .or. & - abs(residuum_plastic(1:sizeDotState,g,i,e)) < & - plasticState(p)%aTolState(1:sizeDotState)) + + crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & + plasticState(p)%dotState(1:sizeDotState,cc), & + plasticState(p)%aTolState(1:sizeDotState)) do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - where(dNeq0(sourceState(p)%p(s)%dotState(1:sizeDotState,cc))) - residuum_source_rel(1:sizeDotState) = residuum_source(1:sizeDotState,s,g,i,e) & - / sourceState(p)%p(s)%state(1:sizeDotState,cc) - else where - residuum_source_rel(1:SizeDotState) = 0.0_pReal - end where - - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & - all(abs(residuum_source_rel(1:sizeDotState)) < & - rTol_crystalliteState .or. & - abs(residuum_source(1:sizeDotState,s,g,i,e)) < & - sourceState(p)%p(s)%aTolState(1:sizeDotState)) + crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and.& + converged(residuum_source(1:sizeDotState,s,g,i,e), & + sourceState(p)%p(s)%dotState(1:sizeDotState,cc), & + sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo endif enddo; enddo; enddo @@ -2159,6 +2138,34 @@ subroutine integrateStateRKCK45() call update_stress(1.0_pReal) call setConvergenceFlag if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief determines whether a point is converged + !-------------------------------------------------------------------------------------------------- + logical pure function converged(residuum,dotState,absoluteTolerance) + use prec, only: & + dNeq0 + use numerics, only: & + rTol_crystalliteState + + implicit none + real(pReal), dimension(:), intent(in) ::& + residuum, dotState, absoluteTolerance + real(pReal), dimension(size(residuum,1)) ::& + residuum_rel + + where(dNeq0(dotState)) + residuum_rel = residuum/dotState + else where + residuum_rel = 0.0_pReal + end where + + converged = all(abs(residuum_rel) < rTol_crystalliteState .or. & + abs(residuum) < absoluteTolerance) + + end function converged end subroutine integrateStateRKCK45 From 64b89484d2b75693b15be359427bb22244b336aa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 12:56:02 +0100 Subject: [PATCH 24/31] logic better visible --- src/crystallite.f90 | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 210bf8198..dc3e5b154 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1872,17 +1872,16 @@ subroutine integrateStateAdaptiveEuler() implicit none real(pReal), dimension(:), intent(in) ::& residuum, dotState, absoluteTolerance - real(pReal), dimension(size(residuum,1)) ::& - residuum_rel - + logical, dimension(size(residuum,1)) ::& + converged_array + where(dNeq0(dotState)) - residuum_rel = residuum/dotState + converged_array = abs(residuum) < absoluteTolerance .or. (abs(residuum/dotState) < rTol_crystalliteState) else where - residuum_rel = 0.0_pReal + converged_array = .true. end where - converged = all(abs(residuum_rel) < rTol_crystalliteState .or. & - abs(residuum) < absoluteTolerance) + converged = all(converged_array) end function converged @@ -2153,17 +2152,16 @@ subroutine integrateStateRKCK45() implicit none real(pReal), dimension(:), intent(in) ::& residuum, dotState, absoluteTolerance - real(pReal), dimension(size(residuum,1)) ::& - residuum_rel - + logical, dimension(size(residuum,1)) ::& + converged_array + where(dNeq0(dotState)) - residuum_rel = residuum/dotState + converged_array = abs(residuum) < absoluteTolerance .or. (abs(residuum/dotState) < rTol_crystalliteState) else where - residuum_rel = 0.0_pReal + converged_array = .true. end where - converged = all(abs(residuum_rel) < rTol_crystalliteState .or. & - abs(residuum) < absoluteTolerance) + converged = all(converged_array) end function converged From 1d88057ce42c7069d5a0b5d6c7ff8c1f13d29589 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 13:24:35 +0100 Subject: [PATCH 25/31] avoid superflous variables --- src/crystallite.f90 | 60 ++++++++++++++++++++------------------------- 1 file changed, 27 insertions(+), 33 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index dc3e5b154..4adae2a19 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1836,9 +1836,9 @@ subroutine integrateStateAdaptiveEuler() 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) = all(converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%dotState(1:sizeDotState,c), & - plasticState(p)%aTolState(1:sizeDotState)) + plasticState(p)%aTolState(1:sizeDotState))) do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState @@ -1847,9 +1847,9 @@ subroutine integrateStateAdaptiveEuler() + 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), & + all(converged(residuum_source(1:sizeDotState,s,g,i,e), & sourceState(p)%p(s)%dotState(1:sizeDotState,c), & - sourceState(p)%p(s)%aTolState(1:sizeDotState)) + sourceState(p)%p(s)%aTolState(1:sizeDotState))) enddo endif @@ -1863,25 +1863,22 @@ subroutine integrateStateAdaptiveEuler() !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- - logical pure function converged(residuum,dotState,absoluteTolerance) + logical pure elemental function converged(residuum,dotState,absoluteTolerance) use prec, only: & - dNeq0 + dEq0 use numerics, only: & rTol_crystalliteState implicit none - real(pReal), dimension(:), intent(in) ::& + real(pReal), intent(in) ::& residuum, dotState, absoluteTolerance - logical, dimension(size(residuum,1)) ::& - converged_array - where(dNeq0(dotState)) - converged_array = abs(residuum) < absoluteTolerance .or. (abs(residuum/dotState) < rTol_crystalliteState) - else where - converged_array = .true. - end where - - converged = all(converged_array) + if (dEq0(dotState)) then + converged = .true. + else + converged = abs(residuum) < absoluteTolerance & + .or. abs(residuum/dotState) < rTol_crystalliteState + endif end function converged @@ -2116,17 +2113,17 @@ subroutine integrateStateRKCK45() sizeDotState = plasticState(p)%sizeDotState - crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & + crystallite_todo(g,i,e) = all(converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%dotState(1:sizeDotState,cc), & - plasticState(p)%aTolState(1:sizeDotState)) + plasticState(p)%aTolState(1:sizeDotState))) do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and.& - converged(residuum_source(1:sizeDotState,s,g,i,e), & + all(converged(residuum_source(1:sizeDotState,s,g,i,e), & sourceState(p)%p(s)%dotState(1:sizeDotState,cc), & - sourceState(p)%p(s)%aTolState(1:sizeDotState)) + sourceState(p)%p(s)%aTolState(1:sizeDotState))) enddo endif enddo; enddo; enddo @@ -2138,30 +2135,27 @@ subroutine integrateStateRKCK45() call setConvergenceFlag if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck - contains + contains !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- - logical pure function converged(residuum,dotState,absoluteTolerance) + logical pure elemental function converged(residuum,dotState,absoluteTolerance) use prec, only: & - dNeq0 + dEq0 use numerics, only: & rTol_crystalliteState implicit none - real(pReal), dimension(:), intent(in) ::& + real(pReal), intent(in) ::& residuum, dotState, absoluteTolerance - logical, dimension(size(residuum,1)) ::& - converged_array - where(dNeq0(dotState)) - converged_array = abs(residuum) < absoluteTolerance .or. (abs(residuum/dotState) < rTol_crystalliteState) - else where - converged_array = .true. - end where - - converged = all(converged_array) + if (dEq0(dotState)) then + converged = .true. + else + converged = abs(residuum) < absoluteTolerance & + .or. abs(residuum/dotState) < rTol_crystalliteState + endif end function converged From fe88e5bf9cda3e38e2c7f46ce058052316f9b465 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 14:52:12 +0100 Subject: [PATCH 26/31] [skip ci] cleaning --- src/crystallite.f90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 4adae2a19..b089e2f77 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2085,13 +2085,12 @@ subroutine integrateStateRKCK45() do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc) ! store Runge-Kutta dotState - + sourceState(p)%p(s)%RKCK45dotState(6,:,cc) = sourceState(p)%p(s)%dotState(:,cc) + residuum_source(1:sizeDotState,s,g,i,e) = & matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & * crystallite_subdt(g,i,e) - sizeDotState = sourceState(p)%p(s)%sizeDotState sourceState(p)%p(s)%dotState(:,cc) = & matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),B) enddo @@ -2124,7 +2123,7 @@ subroutine integrateStateRKCK45() all(converged(residuum_source(1:sizeDotState,s,g,i,e), & sourceState(p)%p(s)%dotState(1:sizeDotState,cc), & sourceState(p)%p(s)%aTolState(1:sizeDotState))) - enddo + enddo endif enddo; enddo; enddo !$OMP END PARALLEL DO From e1c2747393392543bfee7cdcfe25a243d35116d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 16:06:14 +0100 Subject: [PATCH 27/31] logic error for nonlocal --- src/crystallite.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index b089e2f77..3ad592147 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2340,7 +2340,7 @@ subroutine update_dotState(timeFraction) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) !$OMP FLUSH(nonlocalStop) - if (nonlocalStop .or. (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e))) then + if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Fe, & crystallite_Fi(1:3,1:3,g,i,e), & @@ -2399,7 +2399,7 @@ subroutine update_deltaState do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1,homogenization_Ngrains(mesh_element(3,e)) !$OMP FLUSH(nonlocalStop) - if (nonlocalStop .or. (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e))) then + if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then call constitutive_collectDeltaState(math_6toSym33(crystallite_Tstar_v(1:6,g,i,e)), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), & From 3b13a1af6376314ef50f2f240bd2def7f7570c5e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Jan 2019 17:04:58 +0100 Subject: [PATCH 28/31] calculated convergence criteria wrongly --- src/crystallite.f90 | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 3ad592147..749f202e4 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1837,7 +1837,7 @@ subroutine integrateStateAdaptiveEuler() + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e) crystallite_converged(g,i,e) = all(converged(residuum_plastic(1:sizeDotState,g,i,e), & - plasticState(p)%dotState(1:sizeDotState,c), & + plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%aTolState(1:sizeDotState))) do s = 1_pInt, phase_Nsources(p) @@ -1848,7 +1848,7 @@ subroutine integrateStateAdaptiveEuler() crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and.& all(converged(residuum_source(1:sizeDotState,s,g,i,e), & - sourceState(p)%p(s)%dotState(1:sizeDotState,c), & + sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%aTolState(1:sizeDotState))) enddo @@ -1863,21 +1863,21 @@ subroutine integrateStateAdaptiveEuler() !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- - logical pure elemental function converged(residuum,dotState,absoluteTolerance) + logical pure elemental function converged(residuum,state,aTol) use prec, only: & dEq0 use numerics, only: & - rTol_crystalliteState + rTol => rTol_crystalliteState implicit none real(pReal), intent(in) ::& - residuum, dotState, absoluteTolerance + residuum, state, aTol - if (dEq0(dotState)) then - converged = .true. + if (dEq0(state)) then + converged = .true. ! ToDo: intended behavior? Not rely on absoluteTolerance else - converged = abs(residuum) < absoluteTolerance & - .or. abs(residuum/dotState) < rTol_crystalliteState + converged = abs(residuum) < aTol & + .or. abs(residuum/state) < rTol endif end function converged @@ -2113,7 +2113,7 @@ subroutine integrateStateRKCK45() sizeDotState = plasticState(p)%sizeDotState crystallite_todo(g,i,e) = all(converged(residuum_plastic(1:sizeDotState,g,i,e), & - plasticState(p)%dotState(1:sizeDotState,cc), & + plasticState(p)%state(1:sizeDotState,cc), & plasticState(p)%aTolState(1:sizeDotState))) do s = 1_pInt, phase_Nsources(p) @@ -2121,7 +2121,7 @@ subroutine integrateStateRKCK45() crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and.& all(converged(residuum_source(1:sizeDotState,s,g,i,e), & - sourceState(p)%p(s)%dotState(1:sizeDotState,cc), & + sourceState(p)%p(s)%state(1:sizeDotState,cc), & sourceState(p)%p(s)%aTolState(1:sizeDotState))) enddo endif @@ -2139,21 +2139,21 @@ subroutine integrateStateRKCK45() !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- - logical pure elemental function converged(residuum,dotState,absoluteTolerance) + logical pure elemental function converged(residuum,state,aTol) use prec, only: & dEq0 use numerics, only: & - rTol_crystalliteState + rTol => rTol_crystalliteState implicit none real(pReal), intent(in) ::& - residuum, dotState, absoluteTolerance + residuum, state, aTol - if (dEq0(dotState)) then - converged = .true. + if (dEq0(state)) then + converged = .true. ! ToDo: intended behavior? Not rely on absoluteTolerance else - converged = abs(residuum) < absoluteTolerance & - .or. abs(residuum/dotState) < rTol_crystalliteState + converged = abs(residuum) < aTol & + .or. abs(residuum/state) < rTol endif end function converged From 5eaeb37ea48d2d8b23721d981f24cc8a9a25eda7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 31 Jan 2019 05:17:46 +0100 Subject: [PATCH 29/31] just polishing --- src/crystallite.f90 | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 749f202e4..7c99c4d7a 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -819,8 +819,8 @@ subroutine crystallite_stressTangent() crystallite_invFi(1:3,1:3,c,i,e)) & + math_mul33x33(temp_33_3,dLidS(1:3,1:3,p,o)) end forall - lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + & - math_mul3333xx3333(dSdFi,dFidS) + lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) & + + math_mul3333xx3333(dSdFi,dFidS) call math_invert2(temp_99,error,math_identity2nd(9_pInt)+math_3333to99(lhs_3333)) if (error) then @@ -1350,11 +1350,10 @@ logical function integrateStress(& !* calculate Jacobian for correction term if (mod(jacoCounterLp, iJacoLpresiduum) == 0_pInt) then - forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) & - 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 = - dt * dFe_dLp - dRLp_dLp = math_identity2nd(9_pInt) & - - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) + forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) 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 = - dt * dFe_dLp + dRLp_dLp = math_identity2nd(9_pInt) & + - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & @@ -2076,11 +2075,11 @@ subroutine integrateStateRKCK45() plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) residuum_plastic(1:sizeDotState,g,i,e) = & - matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & + matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)),DB) & ! why transpose? Better to transpose constant DB * crystallite_subdt(g,i,e) plasticState(p)%dotState(:,cc) = & - matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)), B) + matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:sizeDotState,cc)), B) ! why transpose? Better to transpose constant B do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState From cbeb3dcff0133022622f1b16a2bf1375f463d4bf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 31 Jan 2019 09:12:44 +0100 Subject: [PATCH 30/31] use the same formulation for convergence every where --- src/crystallite.f90 | 73 +++++++++++++++++++++++---------------------- 1 file changed, 38 insertions(+), 35 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 7c99c4d7a..f9ceab03c 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1535,11 +1535,8 @@ end function integrateStress !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- subroutine integrateStateFPI() - use, intrinsic :: & - IEEE_arithmetic use numerics, only: & - nState, & - rTol_crystalliteState + nState use mesh, only: & mesh_element use material, only: & @@ -1549,7 +1546,6 @@ subroutine integrateStateFPI() phase_Nsources, & homogenization_Ngrains use constitutive, only: & - constitutive_collectDotState, & constitutive_plasticity_maxSizeDotState, & constitutive_source_maxSizeDotState @@ -1635,9 +1631,9 @@ subroutine integrateStateFPI() plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) - crystallite_converged(g,i,e) = all(abs(residuum_plastic(1:sizeDotState)) & - < max(plasticState(p)%aTolState(1:sizeDotState), & - abs(plasticState(p)%state(1:sizeDotState,c)*rTol_crystalliteState))) + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%aTolState(1:sizeDotState)) do s = 1_pInt, phase_Nsources(p) @@ -1659,9 +1655,9 @@ subroutine integrateStateFPI() + sourceState(p)%p(s)%previousDotState(:,c)* (1.0_pReal - zeta) crystallite_converged(g,i,e) = crystallite_converged(g,i,e) .and. & - all(abs(residuum_source(1:sizeDotState)) & - < max(sourceState(p)%p(s)%aTolState(1:sizeDotState), & - abs(sourceState(p)%p(s)%state(1:sizeDotState,c)*rTol_crystalliteState))) + converged(residuum_source(1:sizeDotState), & + sourceState(p)%p(s)%state(1:sizeDotState,c), & + sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo endif enddo; enddo; enddo @@ -1729,6 +1725,23 @@ subroutine integrateStateFPI() endif end function damper + + !-------------------------------------------------------------------------------------------------- + !> @brief determines whether a point is converged + !-------------------------------------------------------------------------------------------------- + logical pure function converged(residuum,state,aTol) + use prec, only: & + dEq0 + use numerics, only: & + rTol => rTol_crystalliteState + + implicit none + real(pReal), intent(in), dimension(:) ::& + residuum, state, aTol + + converged = all(abs(residuum) <= max(aTol, rTol*abs(state))) + + end function converged end subroutine integrateStateFPI @@ -1835,9 +1848,9 @@ subroutine integrateStateAdaptiveEuler() 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) = all(converged(residuum_plastic(1:sizeDotState,g,i,e), & + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%aTolState(1:sizeDotState))) + plasticState(p)%aTolState(1:sizeDotState)) do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState @@ -1846,9 +1859,9 @@ subroutine integrateStateAdaptiveEuler() + 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.& - all(converged(residuum_source(1:sizeDotState,s,g,i,e), & + converged(residuum_source(1:sizeDotState,s,g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%aTolState(1:sizeDotState))) + sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo endif @@ -1862,22 +1875,17 @@ subroutine integrateStateAdaptiveEuler() !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- - logical pure elemental function converged(residuum,state,aTol) + logical pure function converged(residuum,state,aTol) use prec, only: & dEq0 use numerics, only: & rTol => rTol_crystalliteState implicit none - real(pReal), intent(in) ::& + real(pReal), intent(in), dimension(:) ::& residuum, state, aTol - if (dEq0(state)) then - converged = .true. ! ToDo: intended behavior? Not rely on absoluteTolerance - else - converged = abs(residuum) < aTol & - .or. abs(residuum/state) < rTol - endif + converged = all(abs(residuum) <= max(aTol, rTol*abs(state))) end function converged @@ -2111,17 +2119,17 @@ subroutine integrateStateRKCK45() sizeDotState = plasticState(p)%sizeDotState - crystallite_todo(g,i,e) = all(converged(residuum_plastic(1:sizeDotState,g,i,e), & + crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%state(1:sizeDotState,cc), & - plasticState(p)%aTolState(1:sizeDotState))) + plasticState(p)%aTolState(1:sizeDotState)) do s = 1_pInt, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and.& - all(converged(residuum_source(1:sizeDotState,s,g,i,e), & + converged(residuum_source(1:sizeDotState,s,g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,cc), & - sourceState(p)%p(s)%aTolState(1:sizeDotState))) + sourceState(p)%p(s)%aTolState(1:sizeDotState)) enddo endif enddo; enddo; enddo @@ -2138,22 +2146,17 @@ subroutine integrateStateRKCK45() !-------------------------------------------------------------------------------------------------- !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- - logical pure elemental function converged(residuum,state,aTol) + logical pure function converged(residuum,state,aTol) use prec, only: & dEq0 use numerics, only: & rTol => rTol_crystalliteState implicit none - real(pReal), intent(in) ::& + real(pReal), intent(in), dimension(:) ::& residuum, state, aTol - if (dEq0(state)) then - converged = .true. ! ToDo: intended behavior? Not rely on absoluteTolerance - else - converged = abs(residuum) < aTol & - .or. abs(residuum/state) < rTol - endif + converged = all(abs(residuum) <= max(aTol, rTol*abs(state))) end function converged From aabd98bee9fe4d0a8eaec49bc545ebfe6f073b91 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 31 Jan 2019 09:14:02 +0100 Subject: [PATCH 31/31] no need to repeat the same code --- src/crystallite.f90 | 73 +++++++++++---------------------------------- 1 file changed, 18 insertions(+), 55 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index f9ceab03c..45aca46d1 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1725,23 +1725,6 @@ subroutine integrateStateFPI() endif end function damper - - !-------------------------------------------------------------------------------------------------- - !> @brief determines whether a point is converged - !-------------------------------------------------------------------------------------------------- - logical pure function converged(residuum,state,aTol) - use prec, only: & - dEq0 - use numerics, only: & - rTol => rTol_crystalliteState - - implicit none - real(pReal), intent(in), dimension(:) ::& - residuum, state, aTol - - converged = all(abs(residuum) <= max(aTol, rTol*abs(state))) - - end function converged end subroutine integrateStateFPI @@ -1870,25 +1853,6 @@ subroutine integrateStateAdaptiveEuler() if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief determines whether a point is converged - !-------------------------------------------------------------------------------------------------- - logical pure function converged(residuum,state,aTol) - use prec, only: & - dEq0 - use numerics, only: & - rTol => rTol_crystalliteState - - implicit none - real(pReal), intent(in), dimension(:) ::& - residuum, state, aTol - - converged = all(abs(residuum) <= max(aTol, rTol*abs(state))) - - end function converged - end subroutine integrateStateAdaptiveEuler @@ -2141,25 +2105,6 @@ subroutine integrateStateRKCK45() call setConvergenceFlag if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief determines whether a point is converged - !-------------------------------------------------------------------------------------------------- - logical pure function converged(residuum,state,aTol) - use prec, only: & - dEq0 - use numerics, only: & - rTol => rTol_crystalliteState - - implicit none - real(pReal), intent(in), dimension(:) ::& - residuum, state, aTol - - converged = all(abs(residuum) <= max(aTol, rTol*abs(state))) - - end function converged - end subroutine integrateStateRKCK45 @@ -2201,6 +2146,24 @@ subroutine setConvergenceFlag() end subroutine setConvergenceFlag + !-------------------------------------------------------------------------------------------------- + !> @brief determines whether a point is converged + !-------------------------------------------------------------------------------------------------- + logical pure function converged(residuum,state,aTol) + use prec, only: & + dEq0 + use numerics, only: & + rTol => rTol_crystalliteState + + implicit none + real(pReal), intent(in), dimension(:) ::& + residuum, state, aTol + + converged = all(abs(residuum) <= max(aTol, rTol*abs(state))) + + end function converged + + !-------------------------------------------------------------------------------------------------- !> @brief Standard forwarding of state as state = state0 + dotState * (delta t) !--------------------------------------------------------------------------------------------------