diff --git a/code/constitutive.f90 b/code/constitutive.f90 index c48571197..5fe9529d0 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -94,6 +94,7 @@ subroutine constitutive_init mesh_NcpElems, & mesh_element, & mesh_ipNeighborhood, & + mesh_maxNipNeighbors, & FE_Nips, & FE_NipNeighbors, & FE_geomtype @@ -123,6 +124,7 @@ integer(pInt) g, & gMax, & ! maximum number of grains iMax, & ! maximum number of integration points eMax, & ! maximum number of elements + n, & p, & s, & myInstance,& @@ -763,7 +765,7 @@ end subroutine constitutive_hooke_TandItsTangent !* This subroutine contains the constitutive equation for * !* calculating the rate of change of microstructure * !********************************************************************* -subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, orientation, ipc, ip, el) +subroutine constitutive_collectDotState(Tstar_v, Fe, Fp, Temperature, subdt, subfrac, orientation, ipc, ip, el) use prec, only: pReal, pLongInt use debug, only: debug_cumDotStateCalls, & @@ -772,7 +774,8 @@ use debug, only: debug_cumDotStateCalls, & debug_constitutive, & debug_levelBasic use mesh, only: mesh_NcpElems, & - mesh_maxNips + mesh_maxNips, & + mesh_maxNipNeighbors use material, only: phase_plasticity, & material_phase, & homogenization_maxNgrains @@ -796,6 +799,8 @@ integer(pInt), intent(in) :: ipc, & ! component-ID of current integrat el ! current element real(pReal), intent(in) :: Temperature, & subdt ! timestep +real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + subfrac ! subfraction of timestep real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe, & ! elastic deformation gradient Fp ! plastic deformation gradient @@ -831,8 +836,8 @@ select case (phase_plasticity(material_phase(ipc,ip,el))) case (constitutive_nonlocal_label) constitutive_dotState(ipc,ip,el)%p = constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, constitutive_state, & - subdt, orientation, ipc, ip, el) - + constitutive_state0, subdt, subfrac, orientation, ipc, ip, el) + end select if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index c555fdc70..bc16992a9 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -2020,13 +2020,14 @@ endsubroutine !********************************************************************* !* rate of change of microstructure * !********************************************************************* -function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, state, timestep, orientation, g,ip,el) +function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, state, state0, timestep, subfrac, orientation, g,ip,el) use prec, only: pReal, & pInt, & p_vec, & DAMASK_NaN -use numerics, only: numerics_integrationMode +use numerics, only: numerics_integrationMode, & + numerics_timeSyncing use IO, only: IO_error use debug, only: debug_level, & debug_constitutive, & @@ -2048,6 +2049,7 @@ use math, only: math_norm3, & use mesh, only: mesh_NcpElems, & mesh_maxNips, & mesh_element, & + mesh_maxNipNeighbors, & mesh_ipNeighborhood, & mesh_ipVolume, & mesh_ipArea, & @@ -2072,13 +2074,17 @@ integer(pInt), intent(in) :: g, & ! current real(pReal), intent(in) :: Temperature, & ! temperature timestep ! substepped crystallite time increment real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation +real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + subfrac ! fraction of timestep at the beginning of the substepped crystallite time increment real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe, & ! elastic deformation gradient Fp ! plastic deformation gradient real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & orientation ! crystal lattice orientation type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state ! current microstructural state + state, & ! current microstructural state + state0 ! microstructural state at beginning of crystallite increment + !*** input/output variables !*** output variables @@ -2111,11 +2117,15 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance rhoDotAthermalAnnihilation, & ! density evolution by athermal annihilation rhoDotThermalAnnihilation ! density evolution by thermal annihilation real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & - neighboring_rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles) - rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles) + rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles) + rhoSgl0, & ! single dislocation densities at start of cryst inc (positive/negative screw and edge without dipoles) + rhoSglMe, & ! single dislocation densities of central ip (positive/negative screw and edge without dipoles) + neighboring_rhoSgl ! current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & - v, & ! dislocation glide velocity - neighboring_v, & ! dislocation glide velocity + v, & ! current dislocation glide velocity + v0, & ! dislocation glide velocity at start of cryst inc + vMe, & ! dislocation glide velocity of central ip + neighboring_v, & ! dislocation glide velocity of enighboring ip gdot ! shear rates real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & rhoForest, & ! forest dislocation density @@ -2188,7 +2198,6 @@ where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < constitutive_nonlocal rhoDip = 0.0_pReal - !*** sanity check for timestep if (timestep <= 0.0_pReal) then ! if illegal timestep... @@ -2305,6 +2314,17 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then constitutive_nonlocal_dotState = DAMASK_NaN ! -> return NaN and, hence, enforce cutback return endif + + + if (numerics_timeSyncing) then + forall (t = 1_pInt:4_pInt) & + v0(1_pInt:ns,t) = state0(g,ip,el)%p((12_pInt+t)*ns+1_pInt:(13_pInt+t)*ns) + forall (t = 1_pInt:8_pInt) & + rhoSgl0(1_pInt:ns,t) = state0(g,ip,el)%p((t-1_pInt)*ns+1_pInt:t*ns) + where (abs(rhoSgl0) * mesh_ipVolume(ip,el) ** 0.667_pReal < constitutive_nonlocal_significantN(myInstance) & + .or. abs(rhoSgl0) < constitutive_nonlocal_significantRho(myInstance)) & + rhoSgl0 = 0.0_pReal + endif !*** be aware of the definition of lattice_st = lattice_sd x lattice_sn !!! @@ -2352,12 +2372,18 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then endif if (considerEnteringFlux) then - forall (t = 1_pInt:4_pInt) & - neighboring_v(1_pInt:ns,t) = state(g,neighboring_ip,neighboring_el)%p((12_pInt+t)*ns+1_pInt:(13_pInt+t)*ns) - forall (t = 1_pInt:4_pInt) & - neighboring_rhoSgl(1_pInt:ns,t) = max(state(g,neighboring_ip,neighboring_el)%p((t-1_pInt)*ns+1_pInt:t*ns), 0.0_pReal) - forall (t = 5_pInt:8_pInt) & - neighboring_rhoSgl(1_pInt:ns,t) = state(g,neighboring_ip,neighboring_el)%p((t-1_pInt)*ns+1_pInt:t*ns) + if(numerics_timeSyncing .and. (subfrac(g,neighboring_ip,neighboring_el) == 0.0_pReal & + .or. subfrac(g,ip,el) == 0.0_pReal)) then + forall (t = 1_pInt:4_pInt) & + neighboring_v(1_pInt:ns,t) = state0(g,neighboring_ip,neighboring_el)%p((12_pInt+t)*ns+1_pInt:(13_pInt+t)*ns) + forall (t = 1_pInt:8_pInt) & + neighboring_rhoSgl(1_pInt:ns,t) = state0(g,neighboring_ip,neighboring_el)%p((t-1_pInt)*ns+1_pInt:t*ns) + else + forall (t = 1_pInt:4_pInt) & + neighboring_v(1_pInt:ns,t) = state(g,neighboring_ip,neighboring_el)%p((12_pInt+t)*ns+1_pInt:(13_pInt+t)*ns) + forall (t = 1_pInt:8_pInt) & + neighboring_rhoSgl(1_pInt:ns,t) = state(g,neighboring_ip,neighboring_el)%p((t-1_pInt)*ns+1_pInt:t*ns) + endif where (abs(neighboring_rhoSgl) * mesh_ipVolume(neighboring_ip,neighboring_el) ** 0.667_pReal & < constitutive_nonlocal_significantN(myInstance) & .or. abs(neighboring_rhoSgl) < constitutive_nonlocal_significantRho(myInstance)) & @@ -2404,6 +2430,14 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then endif if (considerLeavingFlux) then + if(numerics_timeSyncing .and. (subfrac(g,neighboring_ip,neighboring_el) == 0.0_pReal & + .or. subfrac(g,ip,el) == 0.0_pReal)) then + rhoSglMe = rhoSgl0 + vMe = v0 + else + rhoSglMe = rhoSgl + vMe = v + endif normal_me2neighbor_defConf = math_det33(Favg) * math_mul33x3(math_inv33(math_transpose33(Favg)), & mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!) normal_me2neighbor = math_mul33x3(math_transpose33(my_Fe), normal_me2neighbor_defConf) / math_det33(my_Fe) ! interface normal in my lattice configuration @@ -2411,18 +2445,18 @@ if (.not. phase_localPlasticity(material_phase(g,ip,el))) then normal_me2neighbor = normal_me2neighbor / math_norm3(normal_me2neighbor) ! normalize the surface normal to unit length do s = 1_pInt,ns do t = 1_pInt,4_pInt - c = (t + 1_pInt) / 2_pInt - if (v(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) - if (v(s,t) * neighboring_v(s,t) >= 0.0_pReal) then ! no sign change in flux density + c = (t + 1_pInt) / 2_pInt + if (vMe(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) + if (vMe(s,t) * neighboring_v(s,t) >= 0.0_pReal) then ! no sign change in flux density transmissivity = sum(constitutive_nonlocal_compatibility(c,1_pInt:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor transmissivity = 0.0_pReal endif - lineLength = rhoSgl(s,t) * v(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface + lineLength = rhoSglMe(s,t) * vMe(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current type rhoDotFlux(s,t+4_pInt) = rhoDotFlux(s,t+4_pInt) + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) & - * sign(1.0_pReal, v(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point - lineLength = rhoSgl(s,t+4_pInt) * v(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of deads that wants to leave through this interface + * sign(1.0_pReal, vMe(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point + lineLength = rhoSglMe(s,t+4_pInt) * vMe(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of deads that wants to leave through this interface rhoDotFlux(s,t+4_pInt) = rhoDotFlux(s,t+4_pInt) - lineLength / mesh_ipVolume(ip,el) * transmissivity ! dead dislocations leaving through this interface endif enddo @@ -3583,6 +3617,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) case('accumulatedshear') constitutive_nonlocal_accumulatedShear(1:ns,g,ip,el) = constitutive_nonlocal_accumulatedShear(1:ns,g,ip,el) + sum(gdot,2)*dt + !$OMP FLUSH(constitutive_nonlocal_accumulatedShear) constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = constitutive_nonlocal_accumulatedShear(1:ns,g,ip,el) cs = cs + ns diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 4d90f6c21..97e4c255b 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -99,6 +99,12 @@ module crystallite crystallite_requested, & !< flag to request crystallite calculation crystallite_todo, & !< flag to indicate need for further computation crystallite_converged !< convergence flag + logical, dimension (:,:), allocatable :: & + crystallite_clearToWindForward, & + crystallite_clearToCutback, & + crystallite_syncSubFrac, & + crystallite_syncSubFracCompleted, & + crystallite_neighborEnforcedCutback contains @@ -225,6 +231,11 @@ subroutine crystallite_init(Temperature) allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. allocate(crystallite_todo(gMax,iMax,eMax)); crystallite_todo = .false. allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true. + allocate(crystallite_clearToWindForward(iMax,eMax)); crystallite_clearToWindForward = .true. + allocate(crystallite_syncSubFrac(iMax,eMax)); crystallite_syncSubFrac = .false. + allocate(crystallite_syncSubFracCompleted(iMax,eMax)); crystallite_syncSubFracCompleted = .false. + allocate(crystallite_clearToCutback(iMax,eMax)); crystallite_clearToCutback = .true. + allocate(crystallite_neighborEnforcedCutback(iMax,eMax)); crystallite_neighborEnforcedCutback = .false. allocate(crystallite_output(maxval(crystallite_Noutput), & material_Ncrystallite)) ; crystallite_output = '' allocate(crystallite_sizePostResults(material_Ncrystallite)) ; crystallite_sizePostResults = 0_pInt @@ -449,6 +460,7 @@ use numerics, only: subStepMinCryst, & nCryst, & numerics_integrator, & numerics_integrationMode, & + numerics_timeSyncing, & relevantStrain, & analyticJaco use debug, only: debug_level, & @@ -474,7 +486,9 @@ use FEsolving, only: FEsolving_execElem, & FEsolving_execIP use mesh, only: mesh_element, & mesh_NcpElems, & - mesh_maxNips + mesh_maxNips, & + mesh_ipNeighborhood, & + FE_NipNeighbors use material, only: homogenization_Ngrains, & homogenization_maxNgrains use constitutive, only: constitutive_sizeState, & @@ -488,12 +502,14 @@ use constitutive, only: constitutive_sizeState, & constitutive_dotState_backup, & constitutive_TandItsTangent + implicit none !*** input variables ***! logical, intent(in) :: updateJaco, rate_sensitivity ! flag indicating wehther we want to update the Jacobian (stiffness) or not !*** local variables ***! real(pReal) myPert, & ! perturbation with correct sign - formerSubStep + formerSubStep, & + subFracIntermediate real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient Tstar ! 2nd Piola-Kirchhoff stress tensor @@ -517,6 +533,9 @@ integer(pInt) NiterationCrystallite, & g, & ! grain index k, & l, & + n, & + neighboring_e, & + neighboring_i, & o, & p, & perturbation , & ! loop counter for forward,backward perturbation mode @@ -583,32 +602,207 @@ NiterationCrystallite = 0_pInt numerics_integrationMode = 1_pInt do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2)))) ! cutback loop for crystallites -!$OMP PARALLEL DO PRIVATE(myNgrains,formerSubStep) - do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed - myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do g = 1,myNgrains - - ! --- wind forward --- - - if (crystallite_converged(g,i,e)) then + if (any(.not. crystallite_localPlasticity) .and. numerics_timeSyncing) then + + ! Time synchronization can only be used for nonlocal calculations, and only there it makes sense. + ! The idea is that in nonlocal calculations often the vast amjority of the ips + ! converges in one iteration whereas a small fraction of ips has to do a lot of cutbacks. + ! Hence, we try to minimize the computational effort by just doing a lot of cutbacks + ! in the vicinity of the "bad" ips and leave the easily converged volume more or less as it is. + ! However, some synchronization of the time step has to be done at the border between "bad" ips + ! and the ones that immediately converged. + + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i6)') '<< CRYST >> crystallite iteration ',NiterationCrystallite + !$OMP END CRITICAL (write2out) + endif + + if (any(crystallite_syncSubFrac)) then + + ! Just did a time synchronization. + ! Dont do anything else than winding the synchronizers forward. + + crystallite_clearToWindForward = crystallite_localPlasticity(1,:,:) .or. crystallite_syncSubFrac + crystallite_clearToCutback = crystallite_localPlasticity(1,:,:) + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i6)') '<< CRYST >> time synchronization: wind forward' + !$OMP END CRITICAL (write2out) + endif + + elseif (any(crystallite_syncSubFracCompleted)) then + + ! Just completed a time synchronization. + ! Make sure that the ips that synchronized their time step start non-converged + + where(crystallite_syncSubFracCompleted) & + crystallite_converged(1,:,:) = .false. + crystallite_syncSubFracCompleted = .false. + crystallite_clearToWindForward = crystallite_localPlasticity(1,:,:) + crystallite_clearToCutback = crystallite_localPlasticity(1,:,:) .or. .not. crystallite_converged(1,:,:) + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i6)') '<< CRYST >> time synchronization: done, proceed with cutback' + !$OMP END CRITICAL (write2out) + endif + + else + + ! Normal calculation. + ! If all converged and are at the end of the time increment, then just do a final wind forward. + ! If all converged, but not all reached the end of the time increment, then we only wind + ! those forward that are still on their way, all others have to wait. + ! If some did not converge and all are still at the start of the time increment, + ! then all non-convergers force their converged neighbors to also do a cutback. + ! In case that some ips have already wound forward to an intermediate time (subfrac), + ! then all those ips that converged in the first iteration, but now have a non-converged neighbor + ! have to synchronize their time step to the same intermediate time. If such a synchronization + ! takes place, all other ips have to wait and only the synchronizers do a cutback. In the next + ! iteration those will do a wind forward while all others still wait. + + crystallite_clearToWindForward = crystallite_localPlasticity(1,:,:) + crystallite_clearToCutback = crystallite_localPlasticity(1,:,:) + if (all(crystallite_localPlasticity .or. crystallite_converged)) then + if (all(crystallite_localPlasticity .or. crystallite_subStep + crystallite_subFrac >= 1.0_pReal)) then + crystallite_clearToWindForward = .true. ! final wind forward + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i6)') '<< CRYST >> final wind forward' + !$OMP END CRITICAL (write2out) + endif + else + crystallite_clearToWindForward = crystallite_localPlasticity(1,:,:) .or. crystallite_subStep(1,:,:) < 1.0_pReal + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i6)') '<< CRYST >> wind forward' + !$OMP END CRITICAL (write2out) + endif + endif + else + subFracIntermediate = maxval(crystallite_subFrac, mask=.not.crystallite_localPlasticity) + if (subFracIntermediate == 0.0_pReal) then + crystallite_neighborEnforcedCutback = .false. ! look for ips that require a cutback because of a nonconverged neighbor + !$OMP PARALLEL DO PRIVATE(neighboring_e,neighboring_i) + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + if (.not. crystallite_localPlasticity(1,i,e) .and. crystallite_converged(1,i,e)) then + do n = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) + neighboring_e = mesh_ipNeighborhood(1,n,i,e) + neighboring_i = mesh_ipNeighborhood(2,n,i,e) + if (neighboring_e > 0_pInt .and. neighboring_i > 0_pInt) then + if (.not. crystallite_localPlasticity(1,neighboring_i,neighboring_e) & + .and. .not. crystallite_converged(1,neighboring_i,neighboring_e)) then + crystallite_neighborEnforcedCutback(i,e) = .true. +#ifndef _OPENMP + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & + write(6,'(a12,i5,1x,i2,a,i5,1x,i2)') '<< CRYST >> ', neighboring_e,neighboring_i, & + ' enforced cutback at ',e,i +#endif + endif + endif + enddo + endif + enddo + enddo + !$OMP END PARALLEL DO + where(crystallite_neighborEnforcedCutback) & + crystallite_converged(1,:,:) = .false. + else + crystallite_syncSubFrac = .false. ! look for ips that have to do a time synchronization because of a nonconverged neighbor + !$OMP PARALLEL DO PRIVATE(neighboring_e,neighboring_i) + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + if (.not. crystallite_localPlasticity(1,i,e) .and. crystallite_subFrac(1,i,e) == 0.0_pReal) then + do n = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) + neighboring_e = mesh_ipNeighborhood(1,n,i,e) + neighboring_i = mesh_ipNeighborhood(2,n,i,e) + if (neighboring_e > 0_pInt .and. neighboring_i > 0_pInt) then + if (.not. crystallite_localPlasticity(1,neighboring_i,neighboring_e) & + .and. .not. crystallite_converged(1,neighboring_i,neighboring_e)) then + crystallite_syncSubFrac(i,e) = .true. +#ifndef _OPENMP + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & + write(6,'(a12,i5,1x,i2,a,i5,1x,i2)') '<< CRYST >> ',neighboring_e,neighboring_i, & + ' enforced time synchronization at ',e,i +#endif + endif + endif + enddo + endif + enddo + enddo + !$OMP END PARALLEL DO + where(crystallite_syncSubFrac) & + crystallite_converged(1,:,:) = .false. + endif + where(.not. crystallite_localPlasticity .and. crystallite_subStep < 1.0_pReal) & + crystallite_converged = .false. + if (any(crystallite_syncSubFrac)) then ! have to do syncing now, so all wait except for the synchronizers which do a cutback + crystallite_clearToWindForward = crystallite_localPlasticity(1,:,:) + crystallite_clearToCutback = crystallite_localPlasticity(1,:,:) .or. crystallite_syncSubFrac + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i6)') '<< CRYST >> time synchronization: cutback' + !$OMP END CRITICAL (write2out) + endif + else + where(.not. crystallite_converged(1,:,:)) & + crystallite_clearToCutback = .true. + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,'(a,i6)') '<< CRYST >> cutback' + !$OMP END CRITICAL (write2out) + endif + endif + endif + endif + + ! Make sure that all cutbackers start with the same substep + + where(.not. crystallite_localPlasticity .and. .not. crystallite_converged) & + crystallite_subStep = minval(crystallite_subStep, mask=.not. crystallite_localPlasticity & + .and. .not. crystallite_converged) + + ! Those that do neither wind forward nor cutback are not to do + + where(.not. crystallite_clearToWindForward .and. .not. crystallite_clearToCutback) & + crystallite_todo(1,:,:) = .false. + + endif + + !$OMP PARALLEL DO PRIVATE(myNgrains,formerSubStep) + do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed + myNgrains = homogenization_Ngrains(mesh_element(3,e)) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed + do g = 1,myNgrains + + ! --- wind forward --- + + if (crystallite_converged(g,i,e) .and. crystallite_clearToWindForward(i,e)) then formerSubStep = crystallite_subStep(g,i,e) crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) -!$OMP FLUSH(crystallite_subFrac) - crystallite_subStep(g,i,e) = min( 1.0_pReal - crystallite_subFrac(g,i,e), & - stepIncreaseCryst * crystallite_subStep(g,i,e) ) -!$OMP FLUSH(crystallite_subStep) - if (crystallite_subStep(g,i,e) > 0.0_pReal) then - crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... - crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad -!$OMP FLUSH(crystallite_subF0) - crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad - crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation - crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient - constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure - crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress - crystallite_todo(g,i,e) = .true. -!$OMP FLUSH(crystallite_todo) + !$OMP FLUSH(crystallite_subFrac) + crystallite_subStep(g,i,e) = min(1.0_pReal - crystallite_subFrac(g,i,e), & + stepIncreaseCryst * crystallite_subStep(g,i,e)) + !$OMP FLUSH(crystallite_subStep) + if (crystallite_subStep(g,i,e) > 0.0_pReal) then + crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... + crystallite_subF0(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ...def grad + !$OMP FLUSH(crystallite_subF0) + crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad + crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation + crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient + constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure + crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress + if (crystallite_syncSubFrac(i,e)) then ! if we just did a synchronization of states, then we wind forward without any further time integration + crystallite_syncSubFracCompleted(i,e) = .true. + crystallite_syncSubFrac(i,e) = .false. + crystallite_todo(g,i,e) = .false. + else + crystallite_todo(g,i,e) = .true. + endif + !$OMP FLUSH(crystallite_todo) #ifndef _OPENMP if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & @@ -619,33 +813,37 @@ do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) write(6,*) endif #endif - elseif (formerSubStep > 0.0_pReal) then ! this crystallite just converged - crystallite_todo(g,i,e) = .false. ! so done here -!$OMP FLUSH(crystallite_todo) - if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt) then - !$OMP CRITICAL (distributionCrystallite) - debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) = & - debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) + 1_pInt - !$OMP END CRITICAL (distributionCrystallite) + elseif (formerSubStep > 0.0_pReal) then ! this crystallite just converged for the entire timestep + crystallite_todo(g,i,e) = .false. ! so done here + !$OMP FLUSH(crystallite_todo) + if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt) then + !$OMP CRITICAL (distributionCrystallite) + debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) = & + debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) + 1_pInt + !$OMP END CRITICAL (distributionCrystallite) + endif endif - endif - - ! --- cutback --- - - else - crystallite_subStep(g,i,e) = subStepSizeCryst * crystallite_subStep(g,i,e) ! cut step in half and restore... -!$OMP FLUSH(crystallite_subStep) - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature - crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) ! ...plastic def grad -!$OMP FLUSH(crystallite_Fp) + + ! --- cutback --- + + elseif (.not. crystallite_converged(g,i,e) .and. crystallite_clearToCutback(i,e)) then + if (crystallite_syncSubFrac(i,e)) then ! synchronize time + crystallite_subStep(g,i,e) = subFracIntermediate + else + crystallite_subStep(g,i,e) = subStepSizeCryst * crystallite_subStep(g,i,e) ! cut step in half and restore... + endif + !$OMP FLUSH(crystallite_subStep) + crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature + crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) ! ...plastic def grad + !$OMP FLUSH(crystallite_Fp) crystallite_invFp(1:3,1:3,g,i,e) = math_inv33(crystallite_Fp(1:3,1:3,g,i,e)) -!$OMP FLUSH(crystallite_invFp) + !$OMP FLUSH(crystallite_invFp) crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e) ! ...2nd PK stress ! cant restore dotState here, since not yet calculated in first cutback after initialization crystallite_todo(g,i,e) = crystallite_subStep(g,i,e) > subStepMinCryst ! still on track or already done (beyond repair) -!$OMP FLUSH(crystallite_todo) + !$OMP FLUSH(crystallite_todo) #ifndef _OPENMP if (crystallite_todo(g,i,e) & .and. iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt & @@ -659,15 +857,15 @@ do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) endif ! --- prepare for integration --- - - if (crystallite_todo(g,i,e)) then - crystallite_subF(1:3,1:3,g,i,e) = crystallite_subF0(1:3,1:3,g,i,e) & - + crystallite_subStep(g,i,e) & - * (crystallite_partionedF(1:3,1:3,g,i,e) - crystallite_partionedF0(1:3,1:3,g,i,e)) -!$OMP FLUSH(crystallite_subF) - crystallite_Fe(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) - crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) - crystallite_converged(g,i,e) = .false. ! start out non-converged + + if (crystallite_todo(g,i,e) .and. (crystallite_clearToWindForward(i,e) .or. crystallite_clearToCutback(i,e))) then + crystallite_subF(1:3,1:3,g,i,e) = crystallite_subF0(1:3,1:3,g,i,e) & + + crystallite_subStep(g,i,e) & + * (crystallite_partionedF(1:3,1:3,g,i,e) - crystallite_partionedF0(1:3,1:3,g,i,e)) + !$OMP FLUSH(crystallite_subF) + crystallite_Fe(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) + crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) + crystallite_converged(g,i,e) = .false. ! start out non-converged endif enddo ! grains @@ -675,8 +873,29 @@ do while (any(crystallite_todo(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) enddo ! elements !$OMP END PARALLEL DO + if(numerics_timeSyncing) then + if (any(.not. crystallite_localPlasticity .and. .not. crystallite_todo .and. .not. crystallite_converged & + .and. crystallite_subStep <= subStepMinCryst)) then ! no way of rescuing a nonlocal ip that violated the lower time step limit, ... + where(.not. crystallite_localPlasticity) + crystallite_todo = .false. ! ... so let all nonlocal ips die peacefully + crystallite_subStep = 0.0_pReal + endwhere + endif + endif + + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then + !$OMP CRITICAL (write2out) + write(6,*) + write(6,'(a,e12.5)') '<< CRYST >> min(subStep) ',minval(crystallite_subStep) + write(6,'(a,e12.5)') '<< CRYST >> max(subStep) ',maxval(crystallite_subStep) + write(6,'(a,e12.5)') '<< CRYST >> min(subFrac) ',minval(crystallite_subFrac) + write(6,'(a,e12.5)') '<< CRYST >> max(subFrac) ',maxval(crystallite_subFrac) + write(6,*) + !$OMP END CRITICAL (write2out) + endif + ! --- integrate --- requires fully defined state array (basic + dependent state) - + if (any(crystallite_todo)) then select case(numerics_integrator(numerics_integrationMode)) case(1_pInt) @@ -1068,7 +1287,8 @@ RK4dotTemperature = 0.0_pReal constitutive_RK4dotState(g,i,e)%p = 0.0_pReal ! initialize Runge-Kutta dotState if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) + crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, & + crystallite_orientation, g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif @@ -1193,7 +1413,7 @@ do n = 1_pInt,4_pInt if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep - crystallite_orientation, g,i,e) + crystallite_subFrac, crystallite_orientation, g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif @@ -1390,7 +1610,8 @@ endif do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) + crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, & + crystallite_orientation, g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif @@ -1539,7 +1760,7 @@ do n = 1_pInt,5_pInt if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & crystallite_Temperature(g,i,e), c(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep - crystallite_orientation, g,i,e) + crystallite_subFrac, crystallite_orientation, g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif @@ -1644,9 +1865,9 @@ relTemperatureResiduum = 0.0_pReal mySizeDotState = constitutive_sizeDotState(g,i,e) forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s) - if (crystallite_Temperature(g,i,e) > 0) & + if (crystallite_Temperature(g,i,e) > 0) then relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e) - + endif !$OMP FLUSH(relStateResiduum,relTemperatureResiduum) crystallite_todo(g,i,e) = & @@ -1848,7 +2069,8 @@ if (numerics_integrationMode == 1_pInt) then do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) + crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, & + crystallite_orientation, g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif @@ -1942,7 +2164,8 @@ if (numerics_integrationMode == 1_pInt) then do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) + crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, & + crystallite_orientation, g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif @@ -1988,8 +2211,9 @@ if (numerics_integrationMode == 1_pInt) then forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) & relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s) - if (crystallite_Temperature(g,i,e) > 0_pInt) & + if (crystallite_Temperature(g,i,e) > 0_pInt) then relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e) + endif !$OMP FLUSH(relStateResiduum,relTemperatureResiduum) #ifndef _OPENMP @@ -2076,7 +2300,8 @@ end subroutine crystallite_integrateStateAdaptiveEuler subroutine crystallite_integrateStateEuler(gg,ii,ee) !*** variables and functions from other modules ***! -use numerics, only: numerics_integrationMode +use numerics, only: numerics_integrationMode, & + numerics_timeSyncing use debug, only: debug_level, & debug_crystallite, & debug_levelBasic, & @@ -2138,9 +2363,10 @@ if (numerics_integrationMode == 1_pInt) then !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - if (crystallite_todo(g,i,e)) then + if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) + crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, & + crystallite_orientation, g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif @@ -2148,10 +2374,10 @@ if (numerics_integrationMode == 1_pInt) then !$OMP ENDDO !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - if (crystallite_todo(g,i,e)) then + if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature - if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... + if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped !$OMP END CRITICAL (checkTodo) @@ -2168,7 +2394,7 @@ if (numerics_integrationMode == 1_pInt) then !$OMP DO PRIVATE(mySizeDotState) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - if (crystallite_todo(g,i,e)) then + if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then mySizeDotState = constitutive_sizeDotState(g,i,e) constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) & + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) @@ -2190,14 +2416,15 @@ if (numerics_integrationMode == 1_pInt) then enddo; enddo; enddo !$OMP ENDDO - - ! --- STATE JUMP --- + ! --- STATE JUMP --- + !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - if (crystallite_todo(g,i,e)) then + if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) - if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... + if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local... + .and. .not. numerics_timeSyncing) then !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped !$OMP END CRITICAL (checkTodo) @@ -2211,7 +2438,7 @@ if (numerics_integrationMode == 1_pInt) then !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - if (crystallite_todo(g,i,e)) then + if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states endif @@ -2225,9 +2452,10 @@ endif !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - if (crystallite_todo(g,i,e)) then + if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) - if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... + if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local... + .and. .not. numerics_timeSyncing) then !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped !$OMP END CRITICAL (checkTodo) @@ -2241,7 +2469,7 @@ endif !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains - if (crystallite_todo(g,i,e)) then + if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then crystallite_converged(g,i,e) = .true. ! if still "to do" then converged per definitionem if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then !$OMP CRITICAL (distributionState) @@ -2259,7 +2487,8 @@ endif ! --- CHECK NON-LOCAL CONVERGENCE --- if (.not. singleRun) then ! if not requesting Integration of just a single IP - if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) then ! any non-local not yet converged (or broken)... + if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity) & ! any non-local not yet converged (or broken)... + .and. .not. numerics_timeSyncing) then crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged endif endif @@ -2365,7 +2594,8 @@ endif do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) + crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, & + crystallite_orientation, g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif @@ -2457,7 +2687,8 @@ do while (any(crystallite_todo .and. .not. crystallite_converged) .and. Niterati do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, & - crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_orientation, g,i,e) + crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, & + crystallite_orientation, g,i,e) crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), & crystallite_Temperature(g,i,e),g,i,e) endif diff --git a/code/numerics.f90 b/code/numerics.f90 index 70972b108..bdd5bd7ab 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -72,7 +72,8 @@ real(pReal), protected, public :: & volDiscrMod_RGC = 1.0e+12_pReal, & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint) volDiscrPow_RGC = 5.0_pReal !< powerlaw penalty for volume discrepancy logical, protected, public :: & - analyticJaco = .false. !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations + analyticJaco = .false., & !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations + numerics_timeSyncing = .false. !< flag indicating if time synchronization in crystallite is used for nonlocal plasticity !* Random seeding parameters integer(pInt), protected, public :: & fixedSeed = 0_pInt !< fixed seeding for pseudo-random number generator, Default 0: use random seed @@ -222,6 +223,8 @@ subroutine numerics_init numerics_integrator(2) = IO_intValue(line,positions,2_pInt) case ('analyticjaco') analyticJaco = IO_intValue(line,positions,2_pInt) > 0_pInt + case ('timesyncing') + numerics_timeSyncing = IO_intValue(line,positions,2_pInt) > 0_pInt case ('unitlength') numerics_unitlength = IO_floatValue(line,positions,2_pInt) @@ -337,6 +340,9 @@ subroutine numerics_init write(6,'(a)') ' Initializing PETSc' CHKERRQ(ierr) #endif + +numerics_timeSyncing = numerics_timeSyncing .and. all(numerics_integrator==2_pInt) ! timeSyncing only allowed for explicit Euler integrator + !* writing parameters to output file write(6,'(a24,1x,es8.1)') ' relevantStrain: ',relevantStrain @@ -356,6 +362,7 @@ subroutine numerics_init write(6,'(a24,1x,es8.1)') ' rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,1x,es8.1)') ' aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,2(1x,i8))') ' integrator: ',numerics_integrator + write(6,'(a24,1x,L8)') ' timeSyncing: ',numerics_timeSyncing write(6,'(a24,1x,L8)') ' analytic Jacobian: ',analyticJaco write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength