diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 781448b8c..e8a802b20 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -22,10 +22,10 @@ module crystallite use discretization use lattice use results - - implicit none + + implicit none private - + real(pReal), dimension(:,:,:), allocatable, public :: & crystallite_dt !< requested time increment of each grain real(pReal), dimension(:,:,:), allocatable :: & @@ -33,7 +33,7 @@ module crystallite crystallite_subFrac, & !< already calculated fraction of increment crystallite_subStep !< size of next integration step type(rotation), dimension(:,:,:), allocatable :: & - crystallite_orientation !< current orientation + crystallite_orientation !< current orientation real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_P !< 1st Piola-Kirchhoff stress per grain @@ -74,33 +74,33 @@ module crystallite crystallite_converged, & !< convergence flag crystallite_todo, & !< flag to indicate need for further computation crystallite_localPlasticity !< indicates this grain to have purely local constitutive law - + type :: tOutput !< new requested output (per phase) character(len=pStringLen), allocatable, dimension(:) :: & label end type tOutput type(tOutput), allocatable, dimension(:) :: output_constituent - + type :: tNumerics integer :: & iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp nState, & !< state loop limit - nStress !< stress loop limit + nStress !< stress loop limit real(pReal) :: & subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback subStepSizeCryst, & !< size of first substep when cutback subStepSizeLp, & !< size of first substep when cutback in Lp calculation subStepSizeLi, & !< size of first substep when cutback in Li calculation stepIncreaseCryst, & !< increase of next substep size when previous substep converged - rTol_crystalliteState, & !< relative tolerance in state loop + rTol_crystalliteState, & !< relative tolerance in state loop rTol_crystalliteStress, & !< relative tolerance in stress loop aTol_crystalliteStress !< absolute tolerance in stress loop end type tNumerics - + type(tNumerics) :: num ! numerics parameters. Better name? - + procedure(), pointer :: integrateState - + public :: & crystallite_init, & crystallite_stress, & @@ -116,7 +116,7 @@ contains !> @brief allocates and initialize per grain variables !-------------------------------------------------------------------------------------------------- subroutine crystallite_init - + logical, dimension(discretization_nIP,discretization_nElem) :: devNull integer :: & c, & !< counter in integration point component loop @@ -126,13 +126,13 @@ subroutine crystallite_init iMax, & !< maximum number of integration points eMax, & !< maximum number of elements myNcomponents !< number of components at current IP - + write(6,'(/,a)') ' <<<+- crystallite init -+>>>' - + cMax = homogenization_maxNgrains iMax = discretization_nIP eMax = discretization_nElem - + allocate(crystallite_S0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partionedS0(3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_S(3,3,cMax,iMax,eMax), source=0.0_pReal) @@ -172,23 +172,23 @@ subroutine crystallite_init allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) - + num%subStepMinCryst = config_numerics%getFloat('substepmincryst', defaultVal=1.0e-3_pReal) num%subStepSizeCryst = config_numerics%getFloat('substepsizecryst', defaultVal=0.25_pReal) num%stepIncreaseCryst = config_numerics%getFloat('stepincreasecryst', defaultVal=1.5_pReal) - + num%subStepSizeLp = config_numerics%getFloat('substepsizelp', defaultVal=0.5_pReal) num%subStepSizeLi = config_numerics%getFloat('substepsizeli', defaultVal=0.5_pReal) num%rTol_crystalliteState = config_numerics%getFloat('rtol_crystallitestate', defaultVal=1.0e-6_pReal) num%rTol_crystalliteStress = config_numerics%getFloat('rtol_crystallitestress',defaultVal=1.0e-6_pReal) num%aTol_crystalliteStress = config_numerics%getFloat('atol_crystallitestress',defaultVal=1.0e-8_pReal) - + num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1) - + num%nState = config_numerics%getInt ('nstate', defaultVal=20) num%nStress = config_numerics%getInt ('nstress', defaultVal=40) - + if(num%subStepMinCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinCryst') if(num%subStepSizeCryst <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeCryst') if(num%stepIncreaseCryst <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseCryst') @@ -199,12 +199,12 @@ subroutine crystallite_init if(num%rTol_crystalliteState <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteState') if(num%rTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rTol_crystalliteStress') if(num%aTol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='aTol_crystalliteStress') - + if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') - + if(num%nState < 1) call IO_error(301,ext_msg='nState') if(num%nStress< 1) call IO_error(301,ext_msg='nStress') - + select case(numerics_integrator) case(1) integrateState => integrateStateFPI @@ -217,11 +217,11 @@ subroutine crystallite_init case(5) integrateState => integrateStateRKCK45 end select - + allocate(output_constituent(size(config_phase))) do c = 1, size(config_phase) #if defined(__GFORTRAN__) - allocate(output_constituent(c)%label(1)) + allocate(output_constituent(c)%label(1)) output_constituent(c)%label(1)= 'GfortranBug86277' output_constituent(c)%label = config_phase(c)%getStrings('(output)',defaultVal=output_constituent(c)%label ) if (output_constituent(c)%label (1) == 'GfortranBug86277') output_constituent(c)%label = [character(len=pStringLen)::] @@ -239,6 +239,8 @@ subroutine crystallite_init myNcomponents = homogenization_Ngrains(material_homogenizationAt(e)) do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, myNcomponents crystallite_Fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! plastic def gradient reflects init orientation + crystallite_Fp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) & + / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e) crystallite_F0(1:3,1:3,c,i,e) = math_I3 crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phaseAt(c,e)) @@ -250,16 +252,16 @@ subroutine crystallite_init enddo; enddo enddo !$OMP END PARALLEL DO - + if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? - + crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFi0 = crystallite_Fi0 crystallite_partionedF0 = crystallite_F0 crystallite_partionedF = crystallite_F0 - + call crystallite_orientations() - + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -271,7 +273,7 @@ subroutine crystallite_init enddo enddo !$OMP END PARALLEL DO - + devNull = crystallite_stress() call crystallite_stressTangent @@ -283,7 +285,7 @@ subroutine crystallite_init write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) flush(6) endif - + call debug_info call debug_reset #endif @@ -295,7 +297,7 @@ end subroutine crystallite_init !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) - + logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress real(pReal), intent(in), optional :: & dummyArgumentToPreventInternalCompilerErrorWithGCC @@ -308,7 +310,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) e, & !< counter in element loop startIP, endIP, & s - + #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & .and. FEsolving_execElem(1) <= debug_e & @@ -460,7 +462,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) crystallite_stress = .false. elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) - crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) + crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) enddo enddo elementLooping5 @@ -605,12 +607,12 @@ end subroutine crystallite_stressTangent !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- subroutine crystallite_orientations - + integer & c, & !< counter in integration point component loop i, & !< counter in integration point loop e !< counter in element loop - + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -618,7 +620,7 @@ subroutine crystallite_orientations call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) enddo; enddo; enddo !$OMP END PARALLEL DO - + nonlocalPresent: if (any(plasticState%nonLocal)) then !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -636,7 +638,7 @@ end subroutine crystallite_orientations !> @brief Map 2nd order tensor to reference config !-------------------------------------------------------------------------------------------------- function crystallite_push33ToRef(ipc,ip,el, tensor33) - + real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3), intent(in) :: tensor33 real(pReal), dimension(3,3) :: T @@ -644,7 +646,7 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33) el, & ip, & ipc - + T = matmul(material_orientation0(ipc,ip,el)%asMatrix(), & ! ToDo: initial orientation correct? transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el)))) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) @@ -661,11 +663,11 @@ subroutine crystallite_results real(pReal), allocatable, dimension(:,:,:) :: selected_tensors type(rotation), allocatable, dimension(:) :: selected_rotations character(len=pStringLen) :: group,lattice_label - + do p=1,size(config_name_phase) group = trim('current/constituent')//'/'//trim(config_name_phase(p))//'/generic' - - call results_closeGroup(results_addGroup(group)) + + call results_closeGroup(results_addGroup(group)) do o = 1, size(output_constituent(p)%label) select case (output_constituent(p)%label(o)) @@ -722,19 +724,19 @@ subroutine crystallite_results end select enddo enddo - + contains !------------------------------------------------------------------------------------------------ !> @brief select tensors for output !------------------------------------------------------------------------------------------------ function select_tensors(dataset,instance) - + integer, intent(in) :: instance real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset real(pReal), allocatable, dimension(:,:,:) :: select_tensors integer :: e,i,c,j - + allocate(select_tensors(3,3,count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP)) j=0 @@ -748,20 +750,20 @@ subroutine crystallite_results enddo enddo enddo - + end function select_tensors - - + + !-------------------------------------------------------------------------------------------------- !> @brief select rotations for output -!-------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- function select_rotations(dataset,instance) - + integer, intent(in) :: instance type(rotation), dimension(:,:,:), intent(in) :: dataset type(rotation), allocatable, dimension(:) :: select_rotations integer :: e,i,c,j - + allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP)) j=0 @@ -775,7 +777,7 @@ subroutine crystallite_results enddo enddo enddo - + end function select_rotations end subroutine crystallite_results @@ -786,12 +788,12 @@ end subroutine crystallite_results !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- logical function integrateStress(ipc,ip,el,timeFraction) - + integer, intent(in):: el, & ! element index ip, & ! integration point index ipc ! grain index real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep - + real(pReal), dimension(3,3):: F, & ! deformation gradient at end of timestep Fp_new, & ! plastic deformation gradient at end of timestep Fe_new, & ! elastic deformation gradient at end of timestep @@ -848,7 +850,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) logical :: error external :: & dgesv - + integrateStress = .false. if (present(timeFraction)) then @@ -862,7 +864,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! take as first guess Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! take as first guess - + call math_invert33(invFp_current,devNull,error,crystallite_subFp0(1:3,1:3,ipc,ip,el)) if (error) return ! error call math_invert33(invFi_current,devNull,error,crystallite_subFi0(1:3,1:3,ipc,ip,el)) @@ -883,22 +885,22 @@ logical function integrateStress(ipc,ip,el,timeFraction) invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) Fi_new = math_inv33(invFi_new) detInvFi = math_det33(invFi_new) - + jacoCounterLp = 0 steplengthLp = 1.0_pReal residuumLp_old = 0.0_pReal Lpguess_old = Lpguess - + NiterationStressLp = 0 LpLoop: do NiterationStressLp = NiterationStressLp + 1 if (NiterationStressLp>num%nStress) return ! error - + B = math_I3 - dt*Lpguess Fe = matmul(matmul(A,B), invFi_new) call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi_new, ipc, ip, el) - + call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & S, Fi_new, ipc, ip, el) @@ -906,7 +908,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) aTolLp = max(num%rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error num%aTol_crystalliteStress) ! minimum lower cutoff residuumLp = Lpguess - Lp_constitutive - + if (any(IEEE_is_NaN(residuumLp))) then return ! error elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance @@ -964,11 +966,11 @@ logical function integrateStress(ipc,ip,el,timeFraction) + deltaLi * steplengthLi cycle LiLoop endif - + !* calculate Jacobian for correction term if (mod(jacoCounterLi, num%iJacoLpresiduum) == 0) then jacoCounterLi = jacoCounterLi + 1 - + temp_33 = matmul(matmul(A,B),invFi_current) do o=1,3; do p=1,3 dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) @@ -986,11 +988,11 @@ logical function integrateStress(ipc,ip,el,timeFraction) if (ierr /= 0) return ! error deltaLi = - math_9to33(temp_9) endif - + Liguess = Liguess & + deltaLi * steplengthLi enddo LiLoop - + invFp_new = matmul(invFp_current,B) call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error @@ -1051,7 +1053,7 @@ subroutine integrateStateFPI #endif ! store previousDotState and previousDotState2 - + !$OMP PARALLEL DO PRIVATE(p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -1078,7 +1080,7 @@ subroutine integrateStateFPI call update_dependentState call update_stress(1.0_pReal) call update_dotState(1.0_pReal) - + !$OMP PARALLEL !$OMP DO PRIVATE(sizeDotState,residuum_plastic,residuum_source,zeta,p,c) do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -1091,7 +1093,7 @@ subroutine integrateStateFPI zeta = damper(plasticState(p)%dotState (:,c), & plasticState(p)%previousDotState (:,c), & plasticState(p)%previousDotState2(:,c)) - + residuum_plastic(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - plasticState(p)%subState0(1:sizeDotState,c) & - ( plasticState(p)%dotState (:,c) * zeta & @@ -1099,18 +1101,18 @@ subroutine integrateStateFPI ) * crystallite_subdt(g,i,e) plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & - - residuum_plastic(1:sizeDotState) + - residuum_plastic(1:sizeDotState) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + plasticState(p)%previousDotState(:,c) * (1.0_pReal - zeta) - + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%aTolState(1:sizeDotState)) - + do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + zeta = damper(sourceState(p)%p(s)%dotState (:,c), & sourceState(p)%p(s)%previousDotState (:,c), & sourceState(p)%p(s)%previousDotState2(:,c)) @@ -1181,12 +1183,12 @@ subroutine integrateStateFPI !> @brief calculate the damping for correction of state and dot state !-------------------------------------------------------------------------------------------------- real(pReal) pure function damper(current,previous,previous2) - + real(pReal), dimension(:), intent(in) ::& current, previous, previous2 - + real(pReal) :: dot_prod12, dot_prod22 - + dot_prod12 = dot_product(current - previous, previous - previous2) dot_prod22 = dot_product(previous - previous2, previous - previous2) if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then @@ -1194,7 +1196,7 @@ subroutine integrateStateFPI else damper = 1.0_pReal endif - + end function damper end subroutine integrateStateFPI @@ -1229,7 +1231,7 @@ subroutine integrateStateAdaptiveEuler c, & s, & sizeDotState - + ! 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,discretization_nIP,discretization_nElem) :: & @@ -1250,14 +1252,14 @@ subroutine integrateStateAdaptiveEuler if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - + residuum_plastic(1:sizeDotState,g,i,e) = plasticState(p)%dotstate(1:sizeDotState,c) & * (- 0.5_pReal * crystallite_subdt(g,i,e)) 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, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + residuum_source(1:sizeDotState,s,g,i,e) = sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & * (- 0.5_pReal * crystallite_subdt(g,i,e)) sourceState(p)%p(s)%state(1:sizeDotState,c) = & @@ -1279,17 +1281,17 @@ subroutine integrateStateAdaptiveEuler if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) sizeDotState = plasticState(p)%sizeDotState - + residuum_plastic(1:sizeDotState,g,i,e) = residuum_plastic(1:sizeDotState,g,i,e) & + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e) - + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%aTolState(1:sizeDotState)) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + residuum_source(1:sizeDotState,s,g,i,e) = & residuum_source(1:sizeDotState,s,g,i,e) + 0.5_pReal * sourceState(p)%p(s)%dotState(:,c) * crystallite_subdt(g,i,e) @@ -1298,13 +1300,13 @@ subroutine integrateStateAdaptiveEuler sourceState(p)%p(s)%state(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 - + end subroutine integrateStateAdaptiveEuler @@ -1469,23 +1471,23 @@ subroutine integrateStateRKCK45 do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e) - + sizeDotState = plasticState(p)%sizeDotState - + 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) & ! 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) ! why transpose? Better to transpose constant B - + do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - + 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) @@ -1493,13 +1495,13 @@ subroutine integrateStateRKCK45 sourceState(p)%p(s)%dotState(:,cc) = & matmul(transpose(sourceState(p)%p(s)%RKCK45dotState(1:6,1:sizeDotState,cc)),B) enddo - + endif enddo; enddo; enddo !$OMP END PARALLEL DO call update_state(1.0_pReal) - + ! --- relative residui and state convergence --- !$OMP PARALLEL DO PRIVATE(sizeDotState,p,cc) @@ -1508,16 +1510,16 @@ subroutine integrateStateRKCK45 do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) if (crystallite_todo(g,i,e)) then p = material_phaseAt(g,e); cc = material_phaseMemberAt(g,i,e) - + sizeDotState = plasticState(p)%sizeDotState - + crystallite_todo(g,i,e) = converged(residuum_plastic(1:sizeDotState,g,i,e), & plasticState(p)%state(1:sizeDotState,cc), & plasticState(p)%aTolState(1:sizeDotState)) do s = 1, 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), & sourceState(p)%p(s)%state(1:sizeDotState,cc), & @@ -1532,7 +1534,7 @@ subroutine integrateStateRKCK45 call update_stress(1.0_pReal) call setConvergenceFlag if (any(plasticState(:)%nonlocal)) call nonlocalConvergenceCheck - + end subroutine integrateStateRKCK45 @@ -1541,7 +1543,7 @@ end subroutine integrateStateRKCK45 !> @detail one non-converged nonlocal sets all other nonlocals to non-converged to trigger cut back !-------------------------------------------------------------------------------------------------- subroutine nonlocalConvergenceCheck - + if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... where( .not. crystallite_localPlasticity) crystallite_converged = .false. @@ -1559,7 +1561,7 @@ subroutine setConvergenceFlag e, & !< element index in element loop i, & !< integration point index in ip loop g !< grain index in grain loop - + !OMP DO PARALLEL PRIVATE do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) @@ -1575,7 +1577,7 @@ end subroutine setConvergenceFlag !> @brief determines whether a point is converged !-------------------------------------------------------------------------------------------------- logical pure function converged(residuum,state,aTol) - + real(pReal), intent(in), dimension(:) ::& residuum, state, aTol real(pReal) :: & @@ -1695,11 +1697,11 @@ subroutine update_dotState(timeFraction) g, & !< grain index in grain loop p, & c, & - s + s logical :: & NaN, & nonlocalStop - + nonlocalStop = .false. !$OMP PARALLEL DO PRIVATE (p,c,NaN) @@ -1726,7 +1728,7 @@ subroutine update_dotState(timeFraction) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity + if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity end subroutine update_DotState @@ -1741,11 +1743,11 @@ subroutine update_deltaState mySize, & myOffset, & c, & - s + s logical :: & NaN, & nonlocalStop - + nonlocalStop = .false. !$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize,NaN) @@ -1762,23 +1764,23 @@ subroutine update_deltaState myOffset = plasticState(p)%offsetDeltaState mySize = plasticState(p)%sizeDeltaState NaN = any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c))) - + if (.not. NaN) then - + plasticState(p)%state(myOffset + 1: myOffset + mySize,c) = & plasticState(p)%state(myOffset + 1: myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c) do s = 1, 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(s)%state(myOffset + 1:myOffset + mySize,c) = & sourceState(p)%p(s)%state(myOffset + 1:myOffset + mySize,c) + sourceState(p)%p(s)%deltaState(1:mySize,c) endif enddo endif - + crystallite_todo(g,i,e) = .not. NaN if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken crystallite_converged(g,i,e) = .false. @@ -1788,7 +1790,7 @@ subroutine update_deltaState enddo; enddo; enddo !$OMP END PARALLEL DO if (nonlocalStop) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity - + end subroutine update_deltaState