From a458dc831be5b7148e05e426f61f8abc942e3bc5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 24 Jan 2019 06:56:43 +0100 Subject: [PATCH 1/3] prepare for consistent use of full tensor representation --- src/constitutive.f90 | 20 +++++++++----------- src/crystallite.f90 | 4 +++- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index ccaf01c33..ba6a554a3 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -854,7 +854,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac integer(pInt) :: & ho, & !< homogenization tme, & !< thermal member position - s, & !< counter in source loop + s, & !< counter in source loop instance, of ho = material_homogenizationAt(el) @@ -920,7 +920,7 @@ end subroutine constitutive_collectDotState !> @brief for constitutive models having an instantaneous change of state !> will return false if delta state is not needed/supported by the constitutive model !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el) +subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) use prec, only: & pReal, & pLongInt @@ -929,8 +929,7 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el) debug_constitutive, & debug_levelBasic use math, only: & - math_Mandel6to33, & - math_Mandel33to6, & + math_sym33to6, & math_mul33x33 use material, only: & phasememberAt, & @@ -954,18 +953,17 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el) ipc, & !< component-ID of integration point ip, & !< integration point el !< element - real(pReal), intent(in), dimension(6) :: & - S6 !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), intent(in), dimension(3,3) :: & + S, & !< 2nd Piola Kirchhoff stress Fe, & !< elastic deformation gradient Fi !< intermediate deformation gradient real(pReal), dimension(3,3) :: & Mp integer(pInt) :: & - s, & !< counter in source loop + i, & instance, of - Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) + Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) @@ -975,13 +973,13 @@ subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el) call plastic_kinehardening_deltaState(Mp,instance,of) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_deltaState(math_Mandel33to6(Mp),ip,el) + call plastic_nonlocal_deltaState(math_sym33to6(Mp),ip,el) end select plasticityType - sourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) + sourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) - sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) + sourceType: select case (phase_source(i,material_phase(ipc,ip,el))) case (SOURCE_damage_isoBrittle_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 1fcd92705..01c20a893 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2958,6 +2958,8 @@ logical function stateJump(ipc,ip,el) phaseAt, phasememberAt use constitutive, only: & constitutive_collectDeltaState + use math, only: & + math_6toSym33 implicit none integer(pInt), intent(in):: & @@ -2977,7 +2979,7 @@ logical function stateJump(ipc,ip,el) c = phasememberAt(ipc,ip,el) p = phaseAt(ipc,ip,el) - call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), & + call constitutive_collectDeltaState(math_6toSym33(crystallite_Tstar_v(1:6,ipc,ip,el)), & crystallite_Fe(1:3,1:3,ipc,ip,el), & crystallite_Fi(1:3,1:3,ipc,ip,el), & ipc,ip,el) From a8a5c8eec03500b1276bd91a418e5aef7e7c07cf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 24 Jan 2019 07:12:20 +0100 Subject: [PATCH 2/3] preparing function for deltaState essential a "stateJump" over all points --- src/crystallite.f90 | 76 +++++++++++++++++++++++++++++++-------------- 1 file changed, 52 insertions(+), 24 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 01c20a893..d1565c9aa 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2932,6 +2932,40 @@ subroutine update_dotState(timeFraction) end subroutine update_DotState +subroutine update_deltaState + + implicit none + integer(pInt) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g, & !< grain index in grain loop + p, & + c, & + s + + !$OMP PARALLEL DO + 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 FLUSH(crystallite_todo) + if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then ! converged and still alive... + crystallite_todo(g,i,e) = stateJump(g,i,e) + !$OMP FLUSH(crystallite_todo) + if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken + crystallite_converged(g,i,e) = .false. + if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... + !$OMP CRITICAL (checkTodo) + crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped + !$OMP END CRITICAL (checkTodo) + endif + endif + endif + enddo; enddo; enddo + !$OMP END PARALLEL DO + +end subroutine update_deltaState + + !-------------------------------------------------------------------------------------------------- !> @brief calculates a jump in the state according to the current state and the current stress !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state @@ -2971,10 +3005,8 @@ logical function stateJump(ipc,ip,el) c, & p, & mySource, & - myOffsetPlasticDeltaState, & - myOffsetSourceDeltaState, & - mySizePlasticDeltaState, & - mySizeSourceDeltaState + myOffset, & + mySize c = phasememberAt(ipc,ip,el) p = phaseAt(ipc,ip,el) @@ -2984,44 +3016,40 @@ logical function stateJump(ipc,ip,el) crystallite_Fi(1:3,1:3,ipc,ip,el), & ipc,ip,el) - myOffsetPlasticDeltaState = plasticState(p)%offsetDeltaState - mySizePlasticDeltaState = plasticState(p)%sizeDeltaState + myOffset = plasticState(p)%offsetDeltaState + mySize = plasticState(p)%sizeDeltaState - if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)))) then ! NaN occured in deltaState + if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then ! NaN occured in deltaState stateJump = .false. return endif - plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : & - myOffsetPlasticDeltaState + mySizePlasticDeltaState,c) = & - plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : & - myOffsetPlasticDeltaState + mySizePlasticDeltaState,c) + & - plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) + 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) - myOffsetSourceDeltaState = sourceState(p)%p(mySource)%offsetDeltaState - mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState - if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)))) then ! NaN occured in deltaState + myOffset = sourceState(p)%p(mySource)%offsetDeltaState + mySize = sourceState(p)%p(mySource)%sizeDeltaState + if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySize,c)))) then ! NaN occured in deltaState stateJump = .false. return endif - sourceState(p)%p(mySource)%state(myOffsetSourceDeltaState + 1_pInt : & - myOffsetSourceDeltaState + mySizeSourceDeltaState,c) = & - sourceState(p)%p(mySource)%state(myOffsetSourceDeltaState + 1_pInt : & - myOffsetSourceDeltaState + mySizeSourceDeltaState,c) + & - sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c) + 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) enddo #ifdef DEBUG - if (any(dNeq0(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c))) & + if (any(dNeq0(plasticState(p)%deltaState(1:mySize,c))) & .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then write(6,'(a,i8,1x,i2,1x,i3, /)') '<< CRYST >> update state at el ip ipc ',el,ip,ipc - write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) + write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> deltaState', plasticState(p)%deltaState(1:mySize,c) write(6,'(a,/,(12x,12(e12.5,1x)),/)') '<< CRYST >> new state', & - plasticState(p)%state(myOffsetPlasticDeltaState + 1_pInt : & - myOffsetPlasticDeltaState + mySizePlasticDeltaState,c) + plasticState(p)%state(myOffset + 1_pInt : & + myOffset + mySize,c) endif #endif From c3b48c348440a189b532954715652430ce38a63f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 24 Jan 2019 07:34:30 +0100 Subject: [PATCH 3/3] WIP: update_deltaState --- src/crystallite.f90 | 46 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 44 insertions(+), 2 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d1565c9aa..bfb5bf833 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -2933,22 +2933,64 @@ end subroutine update_DotState subroutine update_deltaState - + use, intrinsic :: & + IEEE_arithmetic + use prec, only: & + dNeq0 +#ifdef DEBUG + use debug, only: & + debug_e, & + debug_i, & + debug_g, & + debug_level, & + debug_crystallite, & + debug_levelExtensive, & + debug_levelSelective +#endif + use material, only: & + plasticState, & + sourceState, & + phase_Nsources, & + phaseAt, phasememberAt + use constitutive, only: & + constitutive_collectDeltaState + use math, only: & + math_6toSym33 implicit none integer(pInt) :: & e, & !< element index in element loop i, & !< integration point index in ip loop g, & !< grain index in grain loop p, & + mySize, & + myOffset, & c, & s + logical :: NaN - !$OMP PARALLEL DO + !$OMP PARALLEL DO PRIVATE(p,c,myOffset,mySize) 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 FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then ! converged and still alive... + 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), & + g,i,e) + p = phaseAt(g,i,e); c = phasememberAt(g,i,e) + 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_pInt: myOffset + mySize,c) = & + plasticState(p)%state(myOffset + 1_pInt: myOffset + mySize,c) + & + plasticState(p)%deltaState(1:mySize,c) + + endif + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken