From 570c0421000461cb193f3032c31f86a01d78a7f5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 07:18:37 +0200 Subject: [PATCH 01/47] update of dependent state always in conjunction with state integration --- src/crystallite.f90 | 35 +++-------------------------------- 1 file changed, 3 insertions(+), 32 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index fb12c439b..d657b5b55 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -849,6 +849,9 @@ logical function integrateStress(ipc,ip,el,timeFraction) F = crystallite_subF(1:3,1:3,ipc,ip,el) endif + call constitutive_dependentState(crystallite_partionedF(1:3,1:3,ipc,ip,el), & + crystallite_Fp(1:3,1:3,ipc,ip,el),ipc,ip,el) + 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 @@ -1060,10 +1063,6 @@ subroutine integrateStateFPI source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) enddo - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. crystallite_todo(g,i,e)) exit iteration @@ -1202,10 +1201,6 @@ subroutine integrateStateEuler nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1281,10 +1276,6 @@ subroutine integrateStateAdaptiveEuler nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1415,10 +1406,6 @@ subroutine integrateStateRK4 * crystallite_subdt(g,i,e) enddo - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1466,14 +1453,6 @@ subroutine integrateStateRK4 nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & - nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1583,10 +1562,6 @@ subroutine integrateStateRKCK45 * crystallite_subdt(g,i,e) enddo - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. @@ -1644,10 +1619,6 @@ subroutine integrateStateRKCK45 nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_dependentState(crystallite_partionedF(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) - crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. From 9e926f1545d2f7218d63b5a0472e5d2221141fbe Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 07:41:55 +0200 Subject: [PATCH 02/47] centralize test for error --- src/constitutive.f90 | 20 +++++++++----- src/crystallite.f90 | 62 ++++++++++++-------------------------------- 2 files changed, 30 insertions(+), 52 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 94e590ab3..4402820e8 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -709,7 +709,7 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) +function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -727,19 +727,22 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, real(pReal), dimension(3,3) :: & Mp integer :: & + phase, & ho, & !< homogenization tme, & !< thermal member position i, & !< counter in source loop instance, of + logical :: broken ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) of = material_phasememberAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + phase = material_phaseAt(ipc,el) + instance = phase_plasticityInstance(phase) Mp = matmul(matmul(transpose(Fi),Fi),S) - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + plasticityType: select case (phase_plasticity(phase)) case (PLASTICITY_ISOTROPIC_ID) plasticityType call plastic_isotropic_dotState (Mp,instance,of) @@ -760,10 +763,11 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, call plastic_nonlocal_dotState (Mp,FArray,FpArray,temperature(ho)%p(tme),subdt, & instance,of,ip,el) end select plasticityType + broken = any(IEEE_is_NaN(plasticState(phase)%dotState(:,of))) - SourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el)) + SourceLoop: do i = 1, phase_Nsources(phase) - sourceType: select case (phase_source(i,material_phaseAt(ipc,el))) + sourceType: select case (phase_source(i,phase)) case (SOURCE_damage_anisoBrittle_ID) sourceType call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? @@ -775,13 +779,15 @@ subroutine constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, call source_damage_anisoDuctile_dotState ( ipc, ip, el) case (SOURCE_thermal_externalheat_ID) sourceType - call source_thermal_externalheat_dotState(material_phaseAt(ipc,el),of) + call source_thermal_externalheat_dotState(phase,of) end select sourceType + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of))) + enddo SourceLoop -end subroutine constitutive_collectDotState +end function constitutive_collectDotState !-------------------------------------------------------------------------------------------------- diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d657b5b55..03de4dea3 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1027,15 +1027,12 @@ subroutine integrateStateFPI p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1066,15 +1063,12 @@ subroutine integrateStateFPI crystallite_todo(g,i,e) = integrateStress(g,i,e) if(.not. crystallite_todo(g,i,e)) exit iteration - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. crystallite_todo(g,i,e)) exit iteration sizeDotState = plasticState(p)%sizeDotState @@ -1172,15 +1166,12 @@ subroutine integrateStateEuler p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1244,15 +1235,12 @@ subroutine integrateStateAdaptiveEuler p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1281,15 +1269,12 @@ subroutine integrateStateAdaptiveEuler nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1361,15 +1346,12 @@ subroutine integrateStateRK4 p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1411,15 +1393,12 @@ subroutine integrateStateRK4 nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e)*CC(stage), g,i,e) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit @@ -1517,15 +1496,12 @@ subroutine integrateStateRKCK45 p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo + if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) cycle @@ -1567,15 +1543,11 @@ subroutine integrateStateRKCK45 nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit - call constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e)*CC(stage), g,i,e) - crystallite_todo(g,i,e) = all(.not. IEEE_is_NaN(plasticState(p)%dotState(:,c))) - do s = 1, phase_Nsources(p) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. all(.not. IEEE_is_NaN(sourceState(p)%p(s)%dotState(:,c))) - enddo if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & nonlocalBroken = .true. if(.not. crystallite_todo(g,i,e)) exit From 6b11d438424ccc74241719434de678d1d156758d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 08:02:08 +0200 Subject: [PATCH 03/47] handle error checking centrally --- src/constitutive.f90 | 22 +++++++++++++++------- src/crystallite.f90 | 13 ++----------- 2 files changed, 17 insertions(+), 18 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 4402820e8..582222dae 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -327,7 +327,7 @@ module constitutive constitutive_initialFi, & constitutive_SandItsTangents, & constitutive_collectDotState, & - constitutive_collectDeltaState, & + constitutive_deltaState, & constitutive_results contains @@ -794,7 +794,7 @@ end function 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(S, Fe, Fi, ipc, ip, el) +function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -808,13 +808,18 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) Mp integer :: & i, & - instance, of + instance, of, & + phase + logical :: & + broken + Mp = matmul(matmul(transpose(Fi),Fi),S) of = material_phasememberAt(ipc,ip,el) + phase = material_phaseAt(ipc,el) instance = phase_plasticityInstance(material_phaseAt(ipc,el)) - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + plasticityType: select case (phase_plasticity(phase)) case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_deltaState(Mp,instance,of) @@ -823,10 +828,11 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) end select plasticityType + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) - sourceLoop: do i = 1, phase_Nsources(material_phaseAt(ipc,el)) + sourceLoop: do i = 1, phase_Nsources(phase) - sourceType: select case (phase_source(i,material_phaseAt(ipc,el))) + sourceType: select case (phase_source(i,phase)) case (SOURCE_damage_isoBrittle_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & @@ -834,9 +840,11 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el) end select sourceType + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) + enddo SourceLoop -end subroutine constitutive_collectDeltaState +end function constitutive_deltaState !-------------------------------------------------------------------------------------------------- diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 03de4dea3..81119f2e2 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1654,18 +1654,15 @@ logical function stateJump(ipc,ip,el) c = material_phaseMemberAt(ipc,ip,el) p = material_phaseAt(ipc,el) - call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), & + stateJump = .not. constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & crystallite_Fe(1:3,1:3,ipc,ip,el), & crystallite_Fi(1:3,1:3,ipc,ip,el), & ipc,ip,el) + if(.not. stateJump) return myOffset = plasticState(p)%offsetDeltaState mySize = plasticState(p)%sizeDeltaState - if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then - stateJump = .false. - return - endif plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = & plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c) @@ -1673,16 +1670,10 @@ logical function stateJump(ipc,ip,el) do mySource = 1, phase_Nsources(p) 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 - stateJump = .false. - return - endif sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = & sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c) enddo - stateJump = .true. - end function stateJump From ec53e4c318328f5d928f45650e09842bda7851ae Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 08:39:49 +0200 Subject: [PATCH 04/47] avoid writing to public variable crystallite_todo --- src/crystallite.f90 | 152 ++++++++++++++++++++++---------------------- 1 file changed, 76 insertions(+), 76 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 81119f2e2..caedf7ef5 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1016,10 +1016,10 @@ subroutine integrateStateFPI real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & - nonlocalBroken + nonlocalBroken, broken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState) + !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1027,15 +1027,15 @@ subroutine integrateStateFPI p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & @@ -1060,16 +1060,16 @@ subroutine integrateStateFPI source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) enddo - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. crystallite_todo(g,i,e)) exit iteration + broken = .not. integrateStress(g,i,e) + if(broken) exit iteration - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(.not. crystallite_todo(g,i,e)) exit iteration + if(broken) exit iteration sizeDotState = plasticState(p)%sizeDotState zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2) @@ -1102,12 +1102,12 @@ subroutine integrateStateFPI enddo if(crystallite_converged(g,i,e)) then - crystallite_todo(g,i,e) = stateJump(g,i,e) + broken = .not. stateJump(g,i,e) exit iteration endif enddo iteration - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. endif @@ -1155,10 +1155,10 @@ subroutine integrateStateEuler s, & sizeDotState logical :: & - nonlocalBroken + nonlocalBroken, broken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c) + !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1166,15 +1166,15 @@ subroutine integrateStateEuler p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & @@ -1187,16 +1187,16 @@ subroutine integrateStateEuler * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. stateJump(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - crystallite_converged(g,i,e) = crystallite_todo(g,i,e) + crystallite_converged(g,i,e) = .not. broken endif enddo; enddo; enddo @@ -1221,13 +1221,13 @@ subroutine integrateStateAdaptiveEuler s, & sizeDotState logical :: & - nonlocalBroken + nonlocalBroken, broken real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1235,15 +1235,15 @@ subroutine integrateStateAdaptiveEuler p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1259,25 +1259,25 @@ subroutine integrateStateAdaptiveEuler + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. stateJump(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1333,12 +1333,12 @@ subroutine integrateStateRK4 s, & sizeDotState logical :: & - nonlocalBroken + nonlocalBroken, broken real(pReal), dimension(constitutive_source_maxSizeDotState,4,maxval(phase_Nsources)) :: source_RK4dotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,4) :: plastic_RK4dotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,source_RK4dotState,plastic_RK4dotState,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1346,15 +1346,15 @@ subroutine integrateStateRK4 p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle do stage = 1,3 sizeDotState = plasticState(p)%sizeDotState @@ -1388,24 +1388,24 @@ subroutine integrateStateRK4 * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e,CC(stage)) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) exit + if(broken) exit - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e)*CC(stage), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) exit + if(broken) exit enddo - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1427,15 +1427,15 @@ subroutine integrateStateRK4 * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. stateJump(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken + crystallite_converged(g,i,e) = .not. broken endif enddo; enddo; enddo @@ -1483,12 +1483,12 @@ subroutine integrateStateRKCK45 s, & sizeDotState logical :: & - nonlocalBroken + nonlocalBroken, broken real(pReal), dimension(constitutive_source_maxSizeDotState,6,maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,6) :: plastic_RKdotState nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState) + !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1496,15 +1496,15 @@ subroutine integrateStateRKCK45 p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle do stage = 1,5 sizeDotState = plasticState(p)%sizeDotState @@ -1538,23 +1538,23 @@ subroutine integrateStateRKCK45 * crystallite_subdt(g,i,e) enddo - crystallite_todo(g,i,e) = integrateStress(g,i,e,CC(stage)) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e,CC(stage)) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) exit + if(broken) exit - crystallite_todo(g,i,e) = .not. constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e)*CC(stage), g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) exit + if(broken) exit enddo - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1563,7 +1563,7 @@ subroutine integrateStateRKCK45 plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + plasticState(p)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - crystallite_todo(g,i,e) = converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) & + broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:6),DB) & * crystallite_subdt(g,i,e), & plasticState(p)%state(1:sizeDotState,c), & plasticState(p)%atol(1:sizeDotState)) @@ -1576,25 +1576,25 @@ subroutine integrateStateRKCK45 sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & + broken = broken .and. .not. & converged(matmul(source_RKdotState(1:sizeDotState,1:6,s),DB) & * crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = stateJump(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. stateJump(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - if(.not. crystallite_todo(g,i,e)) cycle + if(broken) cycle - crystallite_todo(g,i,e) = integrateStress(g,i,e) - if(.not. (crystallite_todo(g,i,e) .or. crystallite_localPlasticity(g,i,e))) & + broken = .not. integrateStress(g,i,e) + if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. - crystallite_converged(g,i,e) = crystallite_todo(g,i,e) ! consider converged if not broken + crystallite_converged(g,i,e) = .not. broken endif enddo; enddo; enddo From d50d55cef3ac57d0e86386c8103be832f889b8c2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 08:54:20 +0200 Subject: [PATCH 05/47] avoid public variables --- src/crystallite.f90 | 50 ++++++++++++++++++++++++++------------------- 1 file changed, 29 insertions(+), 21 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index caedf7ef5..908d57ef8 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -71,7 +71,6 @@ module crystallite crystallite_requested !< used by upper level (homogenization) to request crystallite calculation logical, dimension(:,:,:), allocatable :: & 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) @@ -98,7 +97,7 @@ module crystallite type(tNumerics) :: num ! numerics parameters. Better name? - procedure(), pointer :: integrateState + procedure(integrateStateFPI), pointer :: integrateState public :: & crystallite_init, & @@ -161,7 +160,6 @@ subroutine crystallite_init allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) 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) @@ -301,6 +299,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) e, & !< counter in element loop startIP, endIP, & s + logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false in hase of different Ngrains #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & @@ -344,7 +343,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst - crystallite_todo(c,i,e) = .true. + todo(c,i,e) = .true. crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst endif homogenizationRequestsCalculation enddo; enddo @@ -361,7 +360,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) endif singleRun NiterationCrystallite = 0 - cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) + cutbackLooping: do while (any(todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) NiterationCrystallite = NiterationCrystallite + 1 #ifdef DEBUG @@ -380,8 +379,8 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & num%stepIncreaseCryst * crystallite_subStep(c,i,e)) - crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? - if (crystallite_todo(c,i,e)) then + todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? + if (todo(c,i,e)) then crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) @@ -415,12 +414,12 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) enddo ! cant restore dotState here, since not yet calculated in first cutback after initialization - crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) + todo(c,i,e) = crystallite_subStep(c,i,e) > num%subStepMinCryst ! still on track or already done (beyond repair) endif !-------------------------------------------------------------------------------------------------- ! prepare for integration - if (crystallite_todo(c,i,e)) then + if (todo(c,i,e)) then crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & + crystallite_subStep(c,i,e) *( crystallite_partionedF (1:3,1:3,c,i,e) & -crystallite_partionedF0(1:3,1:3,c,i,e)) @@ -438,9 +437,9 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) !-------------------------------------------------------------------------------------------------- ! integrate --- requires fully defined state array (basic + dependent state) - if (any(crystallite_todo)) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation + if (any(todo)) call integrateState(todo) ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further - crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation + todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation enddo cutbackLooping @@ -998,8 +997,9 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateStateFPI +subroutine integrateStateFPI(todo) + logical, dimension(:,:,:), intent(in) :: todo integer :: & NiterationState, & !< number of iterations in state loop e, & !< element index in element loop @@ -1023,7 +1023,7 @@ subroutine integrateStateFPI do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) @@ -1144,7 +1144,9 @@ end subroutine integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateEuler +subroutine integrateStateEuler(todo) + + logical, dimension(:,:,:), intent(in) :: todo integer :: & e, & !< element index in element loop @@ -1162,7 +1164,7 @@ subroutine integrateStateEuler do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) @@ -1210,7 +1212,9 @@ end subroutine integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -subroutine integrateStateAdaptiveEuler +subroutine integrateStateAdaptiveEuler(todo) + + logical, dimension(:,:,:), intent(in) :: todo integer :: & e, & ! element index in element loop @@ -1231,7 +1235,7 @@ subroutine integrateStateAdaptiveEuler do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) @@ -1309,7 +1313,9 @@ end subroutine integrateStateAdaptiveEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 4th order explicit Runge Kutta method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4 +subroutine integrateStateRK4(todo) + + logical, dimension(:,:,:), intent(in) :: todo real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1342,7 +1348,7 @@ subroutine integrateStateRK4 do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) @@ -1450,7 +1456,9 @@ end subroutine integrateStateRK4 !> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with !> adaptive step size (use 5th order solution to advance = "local extrapolation") !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRKCK45 +subroutine integrateStateRKCK45(todo) + + logical, dimension(:,:,:), intent(in) :: todo real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1492,7 +1500,7 @@ subroutine integrateStateRKCK45 do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(crystallite_todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) From 6eee8f34ac830f4d2db934eb92db66211a4eb886 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 09:01:03 +0200 Subject: [PATCH 06/47] homogeneous mesh --- src/crystallite.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 908d57ef8..4e475343b 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -269,7 +269,7 @@ subroutine crystallite_init #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then write(6,'(a42,1x,i10)') ' # of elements: ', eMax - write(6,'(a42,1x,i10)') 'max # of integration points/element: ', iMax + write(6,'(a42,1x,i10)') ' # of integration points/element: ', iMax write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) flush(6) From 079596f7cdcaa9ac69c7ccaa5260c51bb1725c30 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 09:23:43 +0200 Subject: [PATCH 07/47] unix standard 0/.false. = OK, 1/.true. = not OK --- src/crystallite.f90 | 46 ++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 4e475343b..4fe8c75dd 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -776,7 +776,7 @@ end subroutine crystallite_results !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -logical function integrateStress(ipc,ip,el,timeFraction) +function integrateStress(ipc,ip,el,timeFraction) result(broken) integer, intent(in):: el, & ! element index ip, & ! integration point index @@ -833,11 +833,11 @@ logical function integrateStress(ipc,ip,el,timeFraction) p, & jacoCounterLp, & jacoCounterLi ! counters to check for Jacobian update - logical :: error + logical :: error,broken external :: & dgesv - integrateStress = .false. + broken = .true. if (present(timeFraction)) then dt = crystallite_subdt(ipc,ip,el) * timeFraction @@ -981,7 +981,6 @@ logical function integrateStress(ipc,ip,el,timeFraction) call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error - integrateStress = .true. crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) crystallite_S (1:3,1:3,ipc,ip,el) = S crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess @@ -989,6 +988,7 @@ logical function integrateStress(ipc,ip,el,timeFraction) crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new crystallite_Fe (1:3,1:3,ipc,ip,el) = matmul(matmul(F,invFp_new),invFi_new) + broken = .false. end function integrateStress @@ -1060,7 +1060,7 @@ subroutine integrateStateFPI(todo) source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) enddo - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken) exit iteration broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1068,7 +1068,6 @@ subroutine integrateStateFPI(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken) exit iteration sizeDotState = plasticState(p)%sizeDotState @@ -1102,7 +1101,7 @@ subroutine integrateStateFPI(todo) enddo if(crystallite_converged(g,i,e)) then - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) exit iteration endif @@ -1189,12 +1188,12 @@ subroutine integrateStateEuler(todo) * crystallite_subdt(g,i,e) enddo - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. @@ -1263,12 +1262,12 @@ subroutine integrateStateAdaptiveEuler(todo) + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle @@ -1394,7 +1393,7 @@ subroutine integrateStateRK4(todo) * crystallite_subdt(g,i,e) enddo - broken = .not. integrateStress(g,i,e,CC(stage)) + broken = integrateStress(g,i,e,CC(stage)) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) exit @@ -1433,12 +1432,12 @@ subroutine integrateStateRK4(todo) * crystallite_subdt(g,i,e) enddo - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken @@ -1546,7 +1545,7 @@ subroutine integrateStateRKCK45(todo) * crystallite_subdt(g,i,e) enddo - broken = .not. integrateStress(g,i,e,CC(stage)) + broken = integrateStress(g,i,e,CC(stage)) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) exit @@ -1594,12 +1593,12 @@ subroutine integrateStateRKCK45(todo) nonlocalBroken = .true. if(broken) cycle - broken = .not. stateJump(g,i,e) + broken = stateJump(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. if(broken) cycle - broken = .not. integrateStress(g,i,e) + broken = integrateStress(g,i,e) if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken @@ -1645,7 +1644,7 @@ end function converged !> @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 !-------------------------------------------------------------------------------------------------- -logical function stateJump(ipc,ip,el) +function stateJump(ipc,ip,el) result(broken) integer, intent(in):: & el, & ! element index @@ -1658,15 +1657,16 @@ logical function stateJump(ipc,ip,el) mySource, & myOffset, & mySize + logical :: broken c = material_phaseMemberAt(ipc,ip,el) p = material_phaseAt(ipc,el) - stateJump = .not. constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & - crystallite_Fe(1:3,1:3,ipc,ip,el), & - crystallite_Fi(1:3,1:3,ipc,ip,el), & - ipc,ip,el) - if(.not. stateJump) return + broken = constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & + crystallite_Fe(1:3,1:3,ipc,ip,el), & + crystallite_Fi(1:3,1:3,ipc,ip,el), & + ipc,ip,el) + if(broken) return myOffset = plasticState(p)%offsetDeltaState mySize = plasticState(p)%sizeDeltaState From ce61606c0bd2563b0adeaa10d6d03017b290a5f7 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 09:34:49 +0200 Subject: [PATCH 08/47] not needed --- src/material.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/material.f90 b/src/material.f90 index aefe44878..2a447c036 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -11,7 +11,6 @@ module material use results use IO use debug - use numerics use rotations use discretization From b996b6c42ec1e0954b4e4896cd7950b99f0890e5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 09:56:59 +0200 Subject: [PATCH 09/47] merge functionality --- src/constitutive_plastic_disloUCLA.f90 | 2 +- src/constitutive_plastic_dislotwin.f90 | 2 +- src/constitutive_plastic_isotropic.f90 | 2 +- src/constitutive_plastic_kinehardening.f90 | 2 +- src/constitutive_plastic_none.f90 | 2 +- src/constitutive_plastic_nonlocal.f90 | 2 +- src/constitutive_plastic_phenopowerlaw.f90 | 2 +- src/material.f90 | 65 ++++++---------------- src/source_damage_anisoBrittle.f90 | 2 +- src/source_damage_anisoDuctile.f90 | 2 +- src/source_damage_isoBrittle.f90 | 2 +- src/source_damage_isoDuctile.f90 | 2 +- src/source_thermal_dissipation.f90 | 2 +- src/source_thermal_externalheat.f90 | 2 +- 14 files changed, 31 insertions(+), 60 deletions(-) diff --git a/src/constitutive_plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 index 6be86f266..90a933910 100644 --- a/src/constitutive_plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -209,7 +209,7 @@ module subroutine plastic_disloUCLA_init sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index d36e08846..abc46a45e 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -399,7 +399,7 @@ module subroutine plastic_dislotwin_init + size(['f_tr']) * prm%sum_N_tr sizeState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and atol diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index 94fc9817d..ecf029124 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -117,7 +117,7 @@ module subroutine plastic_isotropic_init sizeDotState = size(['xi ','accumulated_shear']) sizeState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index 5843f5b5e..36b1eedf9 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -164,7 +164,7 @@ module subroutine plastic_kinehardening_init sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl sizeState = sizeDotState + sizeDeltaState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization diff --git a/src/constitutive_plastic_none.f90 b/src/constitutive_plastic_none.f90 index 7ff1c76f7..667fe5638 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -29,7 +29,7 @@ module subroutine plastic_none_init if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - call material_allocatePlasticState(p,NipcMyPhase,0,0,0) + call material_allocateState(plasticState(p),NipcMyPhase,0,0,0) enddo diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 58218ac00..35cc7bdd1 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -384,7 +384,7 @@ module subroutine plastic_nonlocal_init 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure sizeDeltaState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) plasticState(p)%nonlocal = .true. plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index a980d6106..12a30478a 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -213,7 +213,7 @@ module subroutine plastic_phenopowerlaw_init + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw sizeState = sizeDotState - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) !-------------------------------------------------------------------------------------------------- ! state aliases and initialization diff --git a/src/material.f90 b/src/material.f90 index 2a447c036..749c9a3d8 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -173,8 +173,7 @@ module material public :: & material_init, & - material_allocatePlasticState, & - material_allocateSourceState, & + material_allocateState, & ELASTICITY_HOOKE_ID ,& PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & @@ -699,63 +698,35 @@ end subroutine material_parseTexture !-------------------------------------------------------------------------------------------------- -!> @brief allocates the plastic state of a phase +!> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- -subroutine material_allocatePlasticState(phase,NipcMyPhase,& - sizeState,sizeDotState,sizeDeltaState) +subroutine material_allocateState(state, & + NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + class(tState), intent(out) :: & + state integer, intent(in) :: & - phase, & NipcMyPhase, & sizeState, & sizeDotState, & sizeDeltaState - plasticState(phase)%sizeState = sizeState - plasticState(phase)%sizeDotState = sizeDotState - plasticState(phase)%sizeDeltaState = sizeDeltaState - plasticState(phase)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition + state%sizeState = sizeState + state%sizeDotState = sizeDotState + state%sizeDeltaState = sizeDeltaState + state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition - allocate(plasticState(phase)%atol (sizeState), source=0.0_pReal) - allocate(plasticState(phase)%state0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%partionedState0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%state (sizeState,NipcMyPhase), source=0.0_pReal) + allocate(state%atol (sizeState), source=0.0_pReal) + allocate(state%state0 (sizeState,NipcMyPhase), source=0.0_pReal) + allocate(state%partionedState0(sizeState,NipcMyPhase), source=0.0_pReal) + allocate(state%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) + allocate(state%state (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + allocate(state%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) - allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) + allocate(state%deltaState(sizeDeltaState,NipcMyPhase), source=0.0_pReal) -end subroutine material_allocatePlasticState +end subroutine material_allocateState -!-------------------------------------------------------------------------------------------------- -!> @brief allocates the source state of a phase -!-------------------------------------------------------------------------------------------------- -subroutine material_allocateSourceState(phase,of,NipcMyPhase,& - sizeState,sizeDotState,sizeDeltaState) - - integer, intent(in) :: & - phase, & - of, & - NipcMyPhase, & - sizeState, sizeDotState,sizeDeltaState - - sourceState(phase)%p(of)%sizeState = sizeState - sourceState(phase)%p(of)%sizeDotState = sizeDotState - sourceState(phase)%p(of)%sizeDeltaState = sizeDeltaState - sourceState(phase)%p(of)%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition - - allocate(sourceState(phase)%p(of)%atol (sizeState), source=0.0_pReal) - allocate(sourceState(phase)%p(of)%state0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(of)%partionedState0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(of)%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) - allocate(sourceState(phase)%p(of)%state (sizeState,NipcMyPhase), source=0.0_pReal) - - allocate(sourceState(phase)%p(of)%dotState (sizeDotState,NipcMyPhase),source=0.0_pReal) - - allocate(sourceState(phase)%p(of)%deltaState (sizeDeltaState,NipcMyPhase),source=0.0_pReal) - -end subroutine material_allocateSourceState - end module material diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 3978be959..b3af24f38 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -107,7 +107,7 @@ subroutine source_damage_anisoBrittle_init if (any(prm%critDisp < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_critDisp' NipcMyPhase = count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 2b818e2cf..79cc0c2f7 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -89,7 +89,7 @@ subroutine source_damage_anisoDuctile_init if (any(prm%critPlasticStrain < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_criticalplasticstrain' NipcMyPhase=count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) sourceState(p)%p(sourceOffset)%atol = config%getFloat('anisoductile_atol',defaultVal=1.0e-3_pReal) if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index ed6d89a89..9eacb4516 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -83,7 +83,7 @@ subroutine source_damage_isoBrittle_init if (prm%critStrainEnergy <= 0.0_pReal) extmsg = trim(extmsg)//' isobrittle_criticalstrainenergy' NipcMyPhase = count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,1) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,1) sourceState(p)%p(sourceOffset)%atol = config%getFloat('isobrittle_atol',defaultVal=1.0e-3_pReal) if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 7024e595a..96754725d 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -82,7 +82,7 @@ subroutine source_damage_isoDuctile_init if (prm%critPlasticStrain <= 0.0_pReal) extmsg = trim(extmsg)//' isoductile_criticalplasticstrain' NipcMyPhase=count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) sourceState(p)%p(sourceOffset)%atol = config%getFloat('isoductile_atol',defaultVal=1.0e-3_pReal) if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 521c79077..c323e68b5 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -67,7 +67,7 @@ subroutine source_thermal_dissipation_init prm%kappa = config%getFloat('dissipation_coldworkcoeff') NipcMyPhase = count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,0,0,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,0,0,0) end associate enddo diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index ade13bef2..06b8a5197 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -74,7 +74,7 @@ subroutine source_thermal_externalheat_init prm%heat_rate = config%getFloats('externalheat_rate',requiredSize = size(prm%time)) NipcMyPhase = count(material_phaseAt==p) * discretization_nIP - call material_allocateSourceState(p,sourceOffset,NipcMyPhase,1,1,0) + call material_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) end associate enddo From 79012c9ffb0c586ac916b42b73452acdd1d41849 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 10:00:12 +0200 Subject: [PATCH 10/47] not needed --- src/prec.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/prec.f90 b/src/prec.f90 index dd39ddc05..646f7dd69 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -53,8 +53,7 @@ module prec logical :: & nonlocal = .false. real(pReal), pointer, dimension(:,:) :: & - slipRate, & !< slip rate - accumulatedSlip !< accumulated plastic slip + slipRate !< slip rate end type type :: tSourceState From cde8c65bd18b9f41bd669f8db673f86066bcdc74 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 10:57:53 +0200 Subject: [PATCH 11/47] better store data locally --- src/crystallite.f90 | 12 +++++++++--- src/numerics.f90 | 23 ++++++++++------------- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 4fe8c75dd..7b3ee4115 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -15,7 +15,6 @@ module crystallite use DAMASK_interface use config use debug - use numerics use rotations use math use FEsolving @@ -83,7 +82,8 @@ module crystallite integer :: & iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp nState, & !< state loop limit - nStress !< stress loop limit + nStress, & !< stress loop limit + integrator !< integration scheme (ToDo: better use a string) real(pReal) :: & subStepMinCryst, & !< minimum (relative) size of sub-step allowed during cutback subStepSizeCryst, & !< size of first substep when cutback @@ -175,6 +175,8 @@ subroutine crystallite_init num%iJacoLpresiduum = config_numerics%getInt ('ijacolpresiduum', defaultVal=1) + num%integrator = config_numerics%getInt ('integrator', defaultVal=1) + num%nState = config_numerics%getInt ('nstate', defaultVal=20) num%nStress = config_numerics%getInt ('nstress', defaultVal=40) @@ -191,10 +193,14 @@ subroutine crystallite_init if(num%iJacoLpresiduum < 1) call IO_error(301,ext_msg='iJacoLpresiduum') + if(num%integrator < 1 .or. num%integrator > 5) & + call IO_error(301,ext_msg='integrator') + 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) + + select case(num%integrator) case(1) integrateState => integrateStateFPI case(2) diff --git a/src/numerics.f90 b/src/numerics.f90 index a29601322..8d242c71d 100644 --- a/src/numerics.f90 +++ b/src/numerics.f90 @@ -20,8 +20,7 @@ module numerics iJacoStiffness = 1, & !< frequency of stiffness update randomSeed = 0, & !< fixed seeding for pseudo-random number generator, Default 0: use random seed worldrank = 0, & !< MPI worldrank (/=0 for MPI simulations only) - worldsize = 1, & !< MPI worldsize (/=1 for MPI simulations only) - numerics_integrator = 1 !< method used for state integration Default 1: fix-point iteration + worldsize = 1 !< MPI worldsize (/=1 for MPI simulations only) integer(4), protected, public :: & DAMASK_NumThreadsInt = 0 !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive real(pReal), protected, public :: & @@ -134,8 +133,6 @@ subroutine numerics_init defgradTolerance = IO_floatValue(line,chunkPos,2) case ('ijacostiffness') iJacoStiffness = IO_intValue(line,chunkPos,2) - case ('integrator') - numerics_integrator = IO_intValue(line,chunkPos,2) case ('usepingpong') usepingpong = IO_intValue(line,chunkPos,2) > 0 case ('unitlength') @@ -176,6 +173,11 @@ subroutine numerics_init case ('maxstaggerediter') stagItMax = IO_intValue(line,chunkPos,2) +#ifdef PETSC + case ('petsc_options') + petsc_options = trim(line(chunkPos(4):)) +#endif + !-------------------------------------------------------------------------------------------------- ! spectral parameters #ifdef Grid @@ -187,8 +189,6 @@ subroutine numerics_init err_stress_tolrel = IO_floatValue(line,chunkPos,2) case ('err_stress_tolabs') err_stress_tolabs = IO_floatValue(line,chunkPos,2) - case ('petsc_options') - petsc_options = trim(line(chunkPos(4):)) case ('err_curl_tolabs') err_curl_tolAbs = IO_floatValue(line,chunkPos,2) case ('err_curl_tolrel') @@ -206,8 +206,6 @@ subroutine numerics_init integrationorder = IO_intValue(line,chunkPos,2) case ('structorder') structorder = IO_intValue(line,chunkPos,2) - case ('petsc_options') - petsc_options = trim(line(chunkPos(4):)) case ('bbarstabilisation') BBarStabilisation = IO_intValue(line,chunkPos,2) > 0 #endif @@ -223,7 +221,6 @@ subroutine numerics_init ! writing parameters to output write(6,'(a24,1x,es8.1)') ' defgradTolerance: ',defgradTolerance write(6,'(a24,1x,i8)') ' iJacoStiffness: ',iJacoStiffness - write(6,'(a24,1x,i8)') ' integrator: ',numerics_integrator write(6,'(a24,1x,L8)') ' use ping pong scheme: ',usepingpong write(6,'(a24,1x,es8.1,/)')' unitlength: ',numerics_unitlength @@ -266,7 +263,6 @@ subroutine numerics_init write(6,'(a24,1x,es8.1)') ' err_curl_tolRel: ',err_curl_tolRel write(6,'(a24,1x,es8.1)') ' polarAlpha: ',polarAlpha write(6,'(a24,1x,es8.1)') ' polarBeta: ',polarBeta - write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options) #endif !-------------------------------------------------------------------------------------------------- @@ -274,16 +270,17 @@ subroutine numerics_init #ifdef FEM write(6,'(a24,1x,i8)') ' integrationOrder: ',integrationOrder write(6,'(a24,1x,i8)') ' structOrder: ',structOrder - write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options) write(6,'(a24,1x,L8)') ' B-Bar stabilisation: ',BBarStabilisation #endif +#ifdef PETSC + write(6,'(a24,1x,a)') ' PETSc_options: ',trim(petsc_options) +#endif + !-------------------------------------------------------------------------------------------------- ! sanity checks if (defgradTolerance <= 0.0_pReal) call IO_error(301,ext_msg='defgradTolerance') if (iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness') - if (numerics_integrator <= 0 .or. numerics_integrator >= 6) & - call IO_error(301,ext_msg='integrator') if (numerics_unitlength <= 0.0_pReal) call IO_error(301,ext_msg='unitlength') if (residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (itmax <= 1) call IO_error(301,ext_msg='itmax') From 97e89f3f88ec4cb67e76160c5c51d1ac27e544fa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 11:09:30 +0200 Subject: [PATCH 12/47] nonlocal can run in local mode --- src/constitutive_plastic_nonlocal.f90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 35cc7bdd1..99ccf9ec4 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -320,6 +320,8 @@ module subroutine plastic_nonlocal_init prm%fEdgeMultiplication = config%getFloat('edgemultiplication') prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/') + plasticState(p)%nonlocal = config%KeyExists('/nonlocal/') + !-------------------------------------------------------------------------------------------------- ! sanity checks if (any(prm%burgers < 0.0_pReal)) extmsg = trim(extmsg)//' burgers' @@ -386,7 +388,6 @@ module subroutine plastic_nonlocal_init call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) - plasticState(p)%nonlocal = .true. plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) From 5af53f0be75617d1848f3135ddaf9ab534e922fa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 11:30:01 +0200 Subject: [PATCH 13/47] nonlocal is a property of the phase --- src/crystallite.f90 | 99 ++++++++++++++++++--------------------------- 1 file changed, 40 insertions(+), 59 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 7b3ee4115..611a7551d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1029,9 +1029,10 @@ subroutine integrateStateFPI(todo) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) + c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & @@ -1039,8 +1040,7 @@ subroutine integrateStateFPI(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1112,8 +1112,7 @@ subroutine integrateStateFPI(todo) endif enddo iteration - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. endif enddo; enddo; enddo @@ -1169,9 +1168,10 @@ subroutine integrateStateEuler(todo) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) + c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & @@ -1179,8 +1179,7 @@ subroutine integrateStateEuler(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1195,13 +1194,11 @@ subroutine integrateStateEuler(todo) enddo broken = stateJump(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken @@ -1240,9 +1237,10 @@ subroutine integrateStateAdaptiveEuler(todo) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) + c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & @@ -1250,8 +1248,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1269,13 +1266,11 @@ subroutine integrateStateAdaptiveEuler(todo) enddo broken = stateJump(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1284,8 +1279,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1322,16 +1316,16 @@ subroutine integrateStateRK4(todo) logical, dimension(:,:,:), intent(in) :: todo - real(pReal), dimension(3,3), parameter :: & - A = reshape([& + real(pReal), dimension(3,3), parameter :: & + A = reshape([& 0.5_pReal, 0.0_pReal, 0.0_pReal, & 0.0_pReal, 0.5_pReal, 0.0_pReal, & 0.0_pReal, 0.0_pReal, 1.0_pReal], & [3,3]) - real(pReal), dimension(3), parameter :: & - CC = [0.5_pReal, 0.5_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration - real(pReal), dimension(4), parameter :: & - B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] ! weight of slope used for Runge Kutta integration (final weight divided by 6) + real(pReal), dimension(3), parameter :: & + CC = [0.5_pReal, 0.5_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration + real(pReal), dimension(4), parameter :: & + B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] ! weight of slope used for Runge Kutta integration (final weight divided by 6) integer :: & e, & ! element index in element loop @@ -1353,9 +1347,11 @@ subroutine integrateStateRK4(todo) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + + c = material_phaseMemberAt(g,i,e) - p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & @@ -1363,8 +1359,7 @@ subroutine integrateStateRK4(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle do stage = 1,3 @@ -1400,8 +1395,6 @@ subroutine integrateStateRK4(todo) enddo broken = integrateStress(g,i,e,CC(stage)) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. if(broken) exit broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1409,13 +1402,10 @@ subroutine integrateStateRK4(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e)*CC(stage), g,i,e) - - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. if(broken) exit enddo - + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1439,13 +1429,11 @@ subroutine integrateStateRK4(todo) enddo broken = stateJump(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken endif @@ -1505,9 +1493,10 @@ subroutine integrateStateRKCK45(todo) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(todo(g,i,e) .and. (.not. nonlocalBroken .or. crystallite_localPlasticity(g,i,e)) ) then + p = material_phaseAt(g,e) + if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - p = material_phaseAt(g,e); c = material_phaseMemberAt(g,i,e) + c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & @@ -1515,8 +1504,7 @@ subroutine integrateStateRKCK45(todo) crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle do stage = 1,5 @@ -1552,8 +1540,6 @@ subroutine integrateStateRKCK45(todo) enddo broken = integrateStress(g,i,e,CC(stage)) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. if(broken) exit broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1561,12 +1547,10 @@ subroutine integrateStateRKCK45(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e)*CC(stage), g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. if(broken) exit enddo - + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1595,18 +1579,15 @@ subroutine integrateStateRKCK45(todo) sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = stateJump(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. .not. crystallite_localPlasticity(g,i,e)) & - nonlocalBroken = .true. + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken endif From b5efaa08a47ea7e0a662b7b7cbf3ca8e4762aa69 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 12:55:49 +0200 Subject: [PATCH 14/47] use already known mappings --- src/constitutive.f90 | 35 ++++++++++++++++++----------------- src/crystallite.f90 | 35 ++++++++++++++++------------------- 2 files changed, 34 insertions(+), 36 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 582222dae..e7f5628ed 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -709,12 +709,14 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) result(broken) +function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point - el !< element + el, & !< element + phase, & + of real(pReal), intent(in) :: & subdt !< timestep real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & @@ -727,17 +729,14 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el real(pReal), dimension(3,3) :: & Mp integer :: & - phase, & ho, & !< homogenization tme, & !< thermal member position i, & !< counter in source loop - instance, of + instance logical :: broken ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - of = material_phasememberAt(ipc,ip,el) - phase = material_phaseAt(ipc,el) instance = phase_plasticityInstance(phase) Mp = matmul(matmul(transpose(Fi),Fi),S) @@ -794,12 +793,14 @@ end function 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 !-------------------------------------------------------------------------------------------------- -function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) +function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point - el !< element + el, & !< element + phase, & + of real(pReal), intent(in), dimension(3,3) :: & S, & !< 2nd Piola Kirchhoff stress Fe, & !< elastic deformation gradient @@ -808,27 +809,28 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) Mp integer :: & i, & - instance, of, & - phase + instance logical :: & broken - Mp = matmul(matmul(transpose(Fi),Fi),S) - of = material_phasememberAt(ipc,ip,el) - phase = material_phaseAt(ipc,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + instance = phase_plasticityInstance(phase) plasticityType: select case (phase_plasticity(phase)) case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_deltaState(Mp,instance,of) + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + + case default + broken = .false. end select plasticityType - broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + sourceLoop: do i = 1, phase_Nsources(phase) @@ -837,11 +839,10 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) case (SOURCE_damage_isoBrittle_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & ipc, ip, el) + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) end select sourceType - broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) - enddo SourceLoop end function constitutive_deltaState diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 611a7551d..054bc2105 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1038,7 +1038,7 @@ subroutine integrateStateFPI(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1073,7 +1073,7 @@ subroutine integrateStateFPI(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken) exit iteration sizeDotState = plasticState(p)%sizeDotState @@ -1107,7 +1107,7 @@ subroutine integrateStateFPI(todo) enddo if(crystallite_converged(g,i,e)) then - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) exit iteration endif @@ -1177,7 +1177,7 @@ subroutine integrateStateEuler(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1193,7 +1193,7 @@ subroutine integrateStateEuler(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1246,7 +1246,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1265,7 +1265,7 @@ subroutine integrateStateAdaptiveEuler(todo) + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1277,7 +1277,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1357,7 +1357,7 @@ subroutine integrateStateRK4(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1401,7 +1401,7 @@ subroutine integrateStateRK4(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e) + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) if(broken) exit enddo @@ -1428,7 +1428,7 @@ subroutine integrateStateRK4(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1502,7 +1502,7 @@ subroutine integrateStateRKCK45(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1546,7 +1546,7 @@ subroutine integrateStateRKCK45(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e) + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) if(broken) exit enddo @@ -1582,7 +1582,7 @@ subroutine integrateStateRKCK45(todo) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1631,7 +1631,7 @@ end function converged !> @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 !-------------------------------------------------------------------------------------------------- -function stateJump(ipc,ip,el) result(broken) +function stateJump(ipc,ip,el,p,c) result(broken) integer, intent(in):: & el, & ! element index @@ -1646,13 +1646,10 @@ function stateJump(ipc,ip,el) result(broken) mySize logical :: broken - c = material_phaseMemberAt(ipc,ip,el) - p = material_phaseAt(ipc,el) - broken = constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & crystallite_Fe(1:3,1:3,ipc,ip,el), & crystallite_Fi(1:3,1:3,ipc,ip,el), & - ipc,ip,el) + ipc,ip,el,p,c) if(broken) return myOffset = plasticState(p)%offsetDeltaState From 70dd06c4ecca2e9df9fe3ce19c7618dd0cc1aacc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 14:42:38 +0200 Subject: [PATCH 15/47] constitutive should handle state jump alone --- src/constitutive.f90 | 21 ++++++++++++++- src/crystallite.f90 | 62 +++++++++++--------------------------------- 2 files changed, 35 insertions(+), 48 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index e7f5628ed..e2c9dbc05 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -809,7 +809,9 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke Mp integer :: & i, & - instance + instance, & + myOffset, & + mySize logical :: & broken @@ -831,6 +833,17 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke end select plasticityType + if(.not. broken) then + select case(phase_plasticity(phase)) + case (PLASTICITY_NONLOCAL_ID,PLASTICITY_KINEHARDENING_ID) + + myOffset = plasticState(phase)%offsetDeltaState + mySize = plasticState(phase)%sizeDeltaState + plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) = & + plasticState(phase)%state(myOffset + 1:myOffset + mySize,of) + plasticState(phase)%deltaState(1:mySize,of) + end select + endif + sourceLoop: do i = 1, phase_Nsources(phase) @@ -840,6 +853,12 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & ipc, ip, el) broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) + if(.not. broken) then + myOffset = sourceState(phase)%p(i)%offsetDeltaState + mySize = sourceState(phase)%p(i)%sizeDeltaState + sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) = & + sourceState(phase)%p(i)%state(myOffset + 1: myOffset + mySize,of) + sourceState(phase)%p(i)%deltaState(1:mySize,of) + endif end select sourceType diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 054bc2105..8bceff9a4 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1107,7 +1107,9 @@ subroutine integrateStateFPI(todo) enddo if(crystallite_converged(g,i,e)) then - broken = stateJump(g,i,e,p,c) + broken = constitutive_deltaState(crystallite_S(1:3,1:3,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,c) exit iteration endif @@ -1193,7 +1195,9 @@ subroutine integrateStateEuler(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e,p,c) + broken = constitutive_deltaState(crystallite_S(1:3,1:3,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,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1265,7 +1269,9 @@ subroutine integrateStateAdaptiveEuler(todo) + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e,p,c) + broken = constitutive_deltaState(crystallite_S(1:3,1:3,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,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1428,7 +1434,9 @@ subroutine integrateStateRK4(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e,p,c) + broken = constitutive_deltaState(crystallite_S(1:3,1:3,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,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1582,7 +1590,9 @@ subroutine integrateStateRKCK45(todo) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle - broken = stateJump(g,i,e,p,c) + broken = constitutive_deltaState(crystallite_S(1:3,1:3,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,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1627,48 +1637,6 @@ logical pure function converged(residuum,state,atol) end function converged -!-------------------------------------------------------------------------------------------------- -!> @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 -!-------------------------------------------------------------------------------------------------- -function stateJump(ipc,ip,el,p,c) result(broken) - - integer, intent(in):: & - el, & ! element index - ip, & ! integration point index - ipc ! grain index - - integer :: & - c, & - p, & - mySource, & - myOffset, & - mySize - logical :: broken - - broken = constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & - crystallite_Fe(1:3,1:3,ipc,ip,el), & - crystallite_Fi(1:3,1:3,ipc,ip,el), & - ipc,ip,el,p,c) - if(broken) return - - myOffset = plasticState(p)%offsetDeltaState - mySize = plasticState(p)%sizeDeltaState - - - plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = & - plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c) - - do mySource = 1, phase_Nsources(p) - myOffset = sourceState(p)%p(mySource)%offsetDeltaState - mySize = sourceState(p)%p(mySource)%sizeDeltaState - sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = & - sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c) - enddo - -end function stateJump - - !-------------------------------------------------------------------------------------------------- !> @brief Write current restart information (Field and constitutive data) to file. ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively From cf5fcf389b1502da09e86826684563c5391e4d40 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 17:53:07 +0200 Subject: [PATCH 16/47] phase is a property of the element and we have no homogenization for nonlocal --- src/crystallite.f90 | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 8bceff9a4..967a34dbb 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -615,14 +615,16 @@ subroutine crystallite_orientations enddo; enddo; enddo !$OMP END PARALLEL DO - nonlocalPresent: if (any(plasticState%nonLocal)) then + nonlocalPresent: if (any(plasticState%nonlocal)) then !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - if (plasticState(material_phaseAt(1,e))%nonLocal) & + if (plasticState(material_phaseAt(1,e))%nonlocal) then + do i = FEsolving_execIP(1),FEsolving_execIP(2) call plastic_nonlocal_updateCompatibility(crystallite_orientation, & phase_plasticityInstance(material_phaseAt(i,e)),i,e) - enddo; enddo + enddo + endif + enddo !$OMP END PARALLEL DO endif nonlocalPresent From d0d963a2cc8f368aca29ee2dcfd610600331051d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 17:58:50 +0200 Subject: [PATCH 17/47] set independent of number of slip systems --- src/constitutive_plastic_nonlocal.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 99ccf9ec4..d530f6651 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -320,7 +320,6 @@ module subroutine plastic_nonlocal_init prm%fEdgeMultiplication = config%getFloat('edgemultiplication') prm%shortRangeStressCorrection = config%keyExists('/shortrangestresscorrection/') - plasticState(p)%nonlocal = config%KeyExists('/nonlocal/') !-------------------------------------------------------------------------------------------------- ! sanity checks @@ -388,6 +387,7 @@ module subroutine plastic_nonlocal_init call material_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + plasticState(p)%nonlocal = config%KeyExists('/nonlocal/') plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:) From 3a4bb59057bcaeeac4f476a7ef992f64f0825bea Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 1 Apr 2020 18:58:48 +0200 Subject: [PATCH 18/47] no need to store the same information multiple times --- src/crystallite.f90 | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 967a34dbb..0d905aaa2 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -69,8 +69,7 @@ module crystallite logical, dimension(:,:,:), allocatable, public :: & crystallite_requested !< used by upper level (homogenization) to request crystallite calculation logical, dimension(:,:,:), allocatable :: & - crystallite_converged, & !< convergence flag - crystallite_localPlasticity !< indicates this grain to have purely local constitutive law + crystallite_converged !< convergence flag type :: tOutput !< new requested output (per phase) character(len=pStringLen), allocatable, dimension(:) :: & @@ -158,7 +157,6 @@ subroutine crystallite_init allocate(crystallite_orientation(cMax,iMax,eMax)) - allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) @@ -238,7 +236,6 @@ subroutine crystallite_init / 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)) crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), & crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) @@ -248,7 +245,7 @@ subroutine crystallite_init 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? + if(any(plasticState%nonlocal) .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 @@ -277,7 +274,6 @@ subroutine crystallite_init write(6,'(a42,1x,i10)') ' # of elements: ', eMax write(6,'(a42,1x,i10)') ' # of integration points/element: ', iMax write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax - write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) flush(6) endif @@ -1606,7 +1602,7 @@ subroutine integrateStateRKCK45(todo) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalBroken) call nonlocalConvergenceCheck + if(nonlocalBroken) call nonlocalConvergenceCheck end subroutine integrateStateRKCK45 @@ -1617,7 +1613,14 @@ end subroutine integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- subroutine nonlocalConvergenceCheck - where( .not. crystallite_localPlasticity) crystallite_converged = .false. + integer :: e,i,p + + do e = FEsolving_execElem(1),FEsolving_execElem(2) + p = material_phaseAt(1,e) + do i = FEsolving_execIP(1),FEsolving_execIP(2) + if(plasticState(p)%nonlocal) crystallite_converged(1,i,e) = .false. + enddo + enddo end subroutine nonlocalConvergenceCheck From d9806cb7f3ffe9a031b7946f3fbb33c193d4db26 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 2 Apr 2020 00:23:59 +0200 Subject: [PATCH 19/47] do not clutter with nonlocal checks --- src/crystallite.f90 | 30 ++++++------------------------ 1 file changed, 6 insertions(+), 24 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0d905aaa2..7fc51b1c9 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1037,7 +1037,6 @@ subroutine integrateStateFPI(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1113,7 +1112,6 @@ subroutine integrateStateFPI(todo) enddo iteration if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - endif enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1178,7 +1176,6 @@ subroutine integrateStateEuler(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1201,9 +1198,7 @@ subroutine integrateStateEuler(todo) broken = integrateStress(g,i,e) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. - crystallite_converged(g,i,e) = .not. broken - endif enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1239,6 +1234,7 @@ subroutine integrateStateAdaptiveEuler(todo) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) + broken = .false. p = material_phaseAt(g,e) if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then @@ -1249,8 +1245,6 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1270,11 +1264,9 @@ subroutine integrateStateAdaptiveEuler(todo) broken = constitutive_deltaState(crystallite_S(1:3,1:3,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,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1282,13 +1274,10 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & plasticState(p)%state(1:sizeDotState,c), & @@ -1296,7 +1285,6 @@ subroutine integrateStateAdaptiveEuler(todo) do s = 1, phase_Nsources(p) sizeDotState = sourceState(p)%p(s)%sizeDotState - crystallite_converged(g,i,e) = & crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & @@ -1305,6 +1293,7 @@ subroutine integrateStateAdaptiveEuler(todo) enddo endif + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1351,6 +1340,7 @@ subroutine integrateStateRK4(todo) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) + broken = .false. p = material_phaseAt(g,e) if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then @@ -1362,8 +1352,6 @@ subroutine integrateStateRK4(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle do stage = 1,3 @@ -1409,7 +1397,6 @@ subroutine integrateStateRK4(todo) if(broken) exit enddo - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1435,14 +1422,13 @@ subroutine integrateStateRK4(todo) broken = constitutive_deltaState(crystallite_S(1:3,1:3,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,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken endif + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1499,6 +1485,7 @@ subroutine integrateStateRKCK45(todo) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) + broken = .false. p = material_phaseAt(g,e) if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then @@ -1509,8 +1496,6 @@ subroutine integrateStateRKCK45(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle do stage = 1,5 @@ -1556,7 +1541,6 @@ subroutine integrateStateRKCK45(todo) if(broken) exit enddo - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1585,20 +1569,18 @@ subroutine integrateStateRKCK45(todo) sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = constitutive_deltaState(crystallite_S(1:3,1:3,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,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken endif + if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO From 565cf8239f077c8830dd195cf41bef3ae99fcd51 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 2 Apr 2020 09:29:58 +0200 Subject: [PATCH 20/47] can be done in parallel --- src/crystallite.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 7fc51b1c9..595947803 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1597,12 +1597,14 @@ subroutine nonlocalConvergenceCheck integer :: e,i,p + !$OMP PARALLEL DO PRIVATE(p) do e = FEsolving_execElem(1),FEsolving_execElem(2) p = material_phaseAt(1,e) do i = FEsolving_execIP(1),FEsolving_execIP(2) if(plasticState(p)%nonlocal) crystallite_converged(1,i,e) = .false. enddo enddo + !$OMP END PARALLEL DO end subroutine nonlocalConvergenceCheck From b375af83a3dcedc8dae6ef564ab0c1800d396a7a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Apr 2020 15:52:54 +0200 Subject: [PATCH 21/47] bugfix for issue introduced in last merge --- src/crystallite.f90 | 55 +-------------------------------------------- 1 file changed, 1 insertion(+), 54 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 743bcf73b..f0b5ecf08 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -301,7 +301,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) e, & !< counter in element loop startIP, endIP, & s - logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false in has of different Ngrains + logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & @@ -1626,59 +1626,6 @@ logical pure function converged(residuum,state,atol) end function converged -!-------------------------------------------------------------------------------------------------- -!> @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 -!-------------------------------------------------------------------------------------------------- -logical function stateJump(ipc,ip,el) - - integer, intent(in):: & - el, & ! element index - ip, & ! integration point index - ipc ! grain index - - integer :: & - c, & - p, & - mySource, & - myOffset, & - mySize - - c = material_phaseMemberAt(ipc,ip,el) - p = material_phaseAt(ipc,el) - - call constitutive_collectDeltaState(crystallite_S(1:3,1:3,ipc,ip,el), & - crystallite_Fe(1:3,1:3,ipc,ip,el), & - crystallite_Fi(1:3,1:3,ipc,ip,el), & - ipc,ip,el) - - myOffset = plasticState(p)%offsetDeltaState - mySize = plasticState(p)%sizeDeltaState - - if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySize,c)))) then - stateJump = .false. - return - endif - - plasticState(p)%state(myOffset + 1:myOffset + mySize,c) = & - plasticState(p)%state(myOffset + 1:myOffset + mySize,c) + plasticState(p)%deltaState(1:mySize,c) - - do mySource = 1, phase_Nsources(p) - 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 - stateJump = .false. - return - endif - sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) = & - sourceState(p)%p(mySource)%state(myOffset + 1: myOffset + mySize,c) + sourceState(p)%p(mySource)%deltaState(1:mySize,c) - enddo - - stateJump = .true. - -end function stateJump - - !-------------------------------------------------------------------------------------------------- !> @brief Write current restart information (Field and constitutive data) to file. ! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively From 399a0f1b66cd771b60f7f61c113227ea9018e255 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Apr 2020 18:53:59 +0200 Subject: [PATCH 22/47] store information about state size only once --- src/crystallite.f90 | 79 ++++++++++++++++++++++----------------------- 1 file changed, 39 insertions(+), 40 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index f0b5ecf08..f405fc930 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1012,7 +1012,9 @@ subroutine integrateStateFPI(todo) p, & c, & s, & - sizeDotState + size_pl + integer, dimension(maxval(phase_Nsources)) :: & + size_so real(pReal) :: & zeta real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & @@ -1023,7 +1025,7 @@ subroutine integrateStateFPI(todo) nonlocalBroken, broken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(sizeDotState,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState,broken) + !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1040,17 +1042,17 @@ subroutine integrateStateFPI(todo) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - plastic_dotState_p2 = 0.0_pReal * plasticState(p)%dotState (1:sizeDotState,c) ! ToDo can be done smarter/clearer + size_pl = plasticState(p)%sizeDotState + plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & + + plasticState(p)%dotState (1:size_pl,c) & + * crystallite_subdt(g,i,e) + plastic_dotState_p2 = 0.0_pReal * plasticState(p)%dotState (1:size_pl,c) do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - source_dotState(1:sizeDotState,2,s) = 0.0_pReal + size_so(s) = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & + + sourceState(p)%p(s)%dotState (1:size_so(s),c) & + * crystallite_subdt(g,i,e) + source_dotState(1:size_so(s),2,s) = 0.0_pReal enddo iteration: do NiterationState = 1, num%nState @@ -1058,49 +1060,46 @@ subroutine integrateStateFPI(todo) if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 plastic_dotState_p1 = plasticState(p)%dotState(:,c) do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - if(nIterationState > 1) source_dotState(1:sizeDotState,2,s) = source_dotState(1:sizeDotState,1,s) - source_dotState(1:sizeDotState,1,s) = sourceState(p)%p(s)%dotState(:,c) + if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) + source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) enddo broken = integrateStress(g,i,e) if(broken) exit iteration broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken) exit iteration - sizeDotState = plasticState(p)%sizeDotState zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + plastic_dotState_p1 * (1.0_pReal - zeta) - r(1:SizeDotState) = plasticState(p)%state (1:sizeDotState,c) & - - plasticState(p)%subState0(1:sizeDotState,c) & - - plasticState(p)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%state(1:sizeDotState,c) & - - r(1:sizeDotState) - crystallite_converged(g,i,e) = converged(r(1:sizeDotState), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & + - plasticState(p)%subState0(1:size_pl,c) & + - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) + plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & + - r(1:size_pl) + crystallite_converged(g,i,e) = converged(r(1:size_pl), & + plasticState(p)%state(1:size_pl,c), & + plasticState(p)%atol(1:size_pl)) do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState zeta = damper(sourceState(p)%p(s)%dotState(:,c), & - source_dotState(1:sizeDotState,1,s),& - source_dotState(1:sizeDotState,2,s)) + source_dotState(1:size_so(s),1,s),& + source_dotState(1:size_so(s),2,s)) sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState(1:sizeDotState,1,s)* (1.0_pReal - zeta) - r(1:sizeDotState) = sourceState(p)%p(s)%state (1:sizeDotState,c) & - - sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - - sourceState(p)%p(s)%dotState (1:sizeDotState,c) * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%state(1:sizeDotState,c) & - - r(1:sizeDotState) + + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) + r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & + - sourceState(p)%p(s)%subState0(1:size_so(s),c) & + - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & + - r(1:size_so(s)) crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(r(1:sizeDotState), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) + crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & + sourceState(p)%p(s)%state(1:size_so(s),c), & + sourceState(p)%p(s)%atol(1:size_so(s))) enddo if(crystallite_converged(g,i,e)) then From 1d5b1a17cd2795de3cb80eaafbbc6bafaebd0020 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 10 Apr 2020 19:36:29 +0200 Subject: [PATCH 23/47] same style for plastic and source --- src/crystallite.f90 | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index f405fc930..2d4b9785d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1019,13 +1019,14 @@ subroutine integrateStateFPI(todo) zeta real(pReal), dimension(max(constitutive_plasticity_maxSizeDotState,constitutive_source_maxSizeDotState)) :: & r ! state residuum - real(pReal), dimension(:), allocatable :: plastic_dotState_p1, plastic_dotState_p2 + real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: & + plastic_dotState real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & nonlocalBroken, broken nonlocalBroken = .false. - !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState_p1, plastic_dotState_p2,source_dotState,broken) + !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState,source_dotState,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) @@ -1046,7 +1047,7 @@ subroutine integrateStateFPI(todo) plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & + plasticState(p)%dotState (1:size_pl,c) & * crystallite_subdt(g,i,e) - plastic_dotState_p2 = 0.0_pReal * plasticState(p)%dotState (1:size_pl,c) + plastic_dotState(1:size_pl,2) = 0.0_pReal do s = 1, phase_Nsources(p) size_so(s) = sourceState(p)%p(s)%sizeDotState sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & @@ -1057,8 +1058,8 @@ subroutine integrateStateFPI(todo) iteration: do NiterationState = 1, num%nState - if(nIterationState > 1) plastic_dotState_p2 = plastic_dotState_p1 - plastic_dotState_p1 = plasticState(p)%dotState(:,c) + if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) + plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) do s = 1, phase_Nsources(p) if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) @@ -1074,9 +1075,10 @@ subroutine integrateStateFPI(todo) crystallite_subdt(g,i,e), g,i,e,p,c) if(broken) exit iteration - zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState_p1,plastic_dotState_p2) + zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& + plastic_dotState(1:size_pl,2)) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & - + plastic_dotState_p1 * (1.0_pReal - zeta) + + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & - plasticState(p)%subState0(1:size_pl,c) & - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) From 9979eb58f44174e1cf3c1333f5654dd31a896115 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 12:41:03 +0200 Subject: [PATCH 24/47] consistent layout for grid data has now always the shape ([x,y,z,...]) with x fastest. For conversion from or to linear layout ([x*y*z,...]), e.g. storage in ASCII table, reshape needs to have the 'F' option. Credits to Vitesh and Fran for pointing this out. --- processing/post/addCurl.py | 7 ++- processing/post/addDisplacement.py | 14 ++--- processing/post/addDivergence.py | 7 ++- processing/post/addGradient.py | 7 ++- python/damask/_geom.py | 7 +-- python/damask/grid_filters.py | 88 ++++++++++++++++-------------- python/tests/test_grid_filters.py | 14 ++--- 7 files changed, 75 insertions(+), 69 deletions(-) diff --git a/processing/post/addCurl.py b/processing/post/addCurl.py index 87d1ab2f6..17459a2df 100755 --- a/processing/post/addCurl.py +++ b/processing/post/addCurl.py @@ -49,9 +49,10 @@ for name in filenames: for label in options.labels: field = table.get(label) shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor - field = field.reshape(np.append(grid[::-1],shape)) + field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape) + curl = damask.grid_filters.curl(size,field) table.add('curlFFT({})'.format(label), - damask.grid_filters.curl(size[::-1],field).reshape(-1,np.prod(shape)), + curl.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape),order='F'), scriptID+' '+' '.join(sys.argv[1:])) - + table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py index f74d876bc..a9424f5e1 100755 --- a/processing/post/addDisplacement.py +++ b/processing/post/addDisplacement.py @@ -51,23 +51,23 @@ for name in filenames: table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos)) - - F = table.get(options.f).reshape(np.append(grid[::-1],(3,3))) + + F = table.get(options.f).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3)) if options.nodal: - table = damask.Table(damask.grid_filters.node_coord0(grid[::-1],size[::-1]).reshape(-1,3), + table = damask.Table(damask.grid_filters.node_coord0(grid,size).reshape(-1,3,order='F'), {'pos':(3,)}) table.add('avg({}).{}'.format(options.f,options.pos), - damask.grid_filters.node_displacement_avg(size[::-1],F).reshape(-1,3), + damask.grid_filters.node_displacement_avg(size,F).reshape(-1,3,order='F'), scriptID+' '+' '.join(sys.argv[1:])) table.add('fluct({}).{}'.format(options.f,options.pos), - damask.grid_filters.node_displacement_fluct(size[::-1],F).reshape(-1,3), + damask.grid_filters.node_displacement_fluct(size,F).reshape(-1,3,order='F'), scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'_nodal.txt') else: table.add('avg({}).{}'.format(options.f,options.pos), - damask.grid_filters.cell_displacement_avg(size[::-1],F).reshape(-1,3), + damask.grid_filters.cell_displacement_avg(size,F).reshape(-1,3,order='F'), scriptID+' '+' '.join(sys.argv[1:])) table.add('fluct({}).{}'.format(options.f,options.pos), - damask.grid_filters.cell_displacement_fluct(size[::-1],F).reshape(-1,3), + damask.grid_filters.cell_displacement_fluct(size,F).reshape(-1,3,order='F'), scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addDivergence.py b/processing/post/addDivergence.py index 2619bc499..50048b44e 100755 --- a/processing/post/addDivergence.py +++ b/processing/post/addDivergence.py @@ -49,9 +49,10 @@ for name in filenames: for label in options.labels: field = table.get(label) shape = (3,) if np.prod(field.shape)//np.prod(grid) == 3 else (3,3) # vector or tensor - field = field.reshape(np.append(grid[::-1],shape)) + field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape) + div = damask.grid_filters.divergence(size,field) table.add('divFFT({})'.format(label), - damask.grid_filters.divergence(size[::-1],field).reshape(-1,np.prod(shape)//3), + div.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)//3,order='F'), scriptID+' '+' '.join(sys.argv[1:])) - + table.to_ASCII(sys.stdout if name is None else name) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index 409b2ce6d..c6df0eacf 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -49,9 +49,10 @@ for name in filenames: for label in options.labels: field = table.get(label) shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,) # scalar or vector - field = field.reshape(np.append(grid[::-1],shape)) + field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape) + grad = damask.grid_filters.gradient(size,field) table.add('gradFFT({})'.format(label), - damask.grid_filters.gradient(size[::-1],field).reshape(-1,np.prod(shape)*3), + grad.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)*3,order='F'), scriptID+' '+' '.join(sys.argv[1:])) - + table.to_ASCII(sys.stdout if name is None else name) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index e2c2428fe..a33e855d8 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -357,7 +357,6 @@ class Geom: seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]]))) coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F') - else: weights_p = weights.flatten() seeds_p = seeds @@ -370,10 +369,10 @@ class Geom: microstructure = np.array(result.get()) if periodic: - microstructure = microstructure.reshape(grid*3) + microstructure = microstructure.reshape(grid*3,order='F') microstructure = microstructure[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] else: - microstructure = microstructure.reshape(grid) + microstructure = microstructure.reshape(grid,order='F') #comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version) return Geom(microstructure+1,size,homogenization=1) @@ -401,7 +400,7 @@ class Geom: devNull,microstructure = KDTree.query(coords) #comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version) - return Geom(microstructure.reshape(grid)+1,size,homogenization=1) + return Geom(microstructure.reshape(grid,order='F')+1,size,homogenization=1) def to_file(self,fname,pack=None): diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index d8b136a6b..0c1c3867a 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -7,7 +7,7 @@ def _ks(size,grid,first_order=False): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. """ @@ -19,8 +19,7 @@ def _ks(size,grid,first_order=False): k_si = _np.arange(grid[2]//2+1)/size[2] - kk, kj, ki = _np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij') - return _np.concatenate((ki[:,:,:,None],kj[:,:,:,None],kk[:,:,:,None]),axis = 3) + return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1) def curl(size,field): @@ -29,7 +28,7 @@ def curl(size,field): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. """ @@ -53,7 +52,7 @@ def divergence(size,field): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. """ @@ -73,7 +72,7 @@ def gradient(size,field): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. """ @@ -93,9 +92,9 @@ def cell_coord0(grid,size,origin=_np.zeros(3)): Parameters ---------- - grid : numpy.ndarray + grid : numpy.ndarray of shape (3) number of grid points. - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. origin : numpy.ndarray, optional physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. @@ -103,7 +102,11 @@ def cell_coord0(grid,size,origin=_np.zeros(3)): """ start = origin + size/grid*.5 end = origin + size - size/grid*.5 - return _np.mgrid[start[0]:end[0]:grid[0]*1j,start[1]:end[1]:grid[1]*1j,start[2]:end[2]:grid[2]*1j].T + + return _np.stack(_np.meshgrid(_np.linspace(start[0],end[0],grid[0]), + _np.linspace(start[1],end[1],grid[1]), + _np.linspace(start[2],end[2],grid[2]),indexing = 'ij'), + axis = -1) def cell_displacement_fluct(size,F): @@ -112,7 +115,7 @@ def cell_displacement_fluct(size,F): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. F : numpy.ndarray deformation gradient field. @@ -139,14 +142,14 @@ def cell_displacement_avg(size,F): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. F : numpy.ndarray deformation gradient field. """ F_avg = _np.average(F,axis=(0,1,2)) - return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),cell_coord0(F.shape[:3][::-1],size)) + return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),cell_coord0(F.shape[:3],size)) def cell_displacement(size,F): @@ -155,7 +158,7 @@ def cell_displacement(size,F): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. F : numpy.ndarray deformation gradient field. @@ -170,30 +173,30 @@ def cell_coord(size,F,origin=_np.zeros(3)): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. F : numpy.ndarray deformation gradient field. - origin : numpy.ndarray, optional + origin : numpy.ndarray of shape (3), optional physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. """ - return cell_coord0(F.shape[:3][::-1],size,origin) + cell_displacement(size,F) + return cell_coord0(F.shape[:3],size,origin) + cell_displacement(size,F) def cell_coord0_gridSizeOrigin(coord0,ordered=True): """ - Return grid 'DNA', i.e. grid, size, and origin from array of cell positions. + Return grid 'DNA', i.e. grid, size, and origin from 1D array of cell positions. Parameters ---------- - coord0 : numpy.ndarray - array of undeformed cell coordinates. + coord0 : numpy.ndarray of shape (:,3) + undeformed cell coordinates. ordered : bool, optional expect coord0 data to be ordered (x fast, z slow). """ - coords = [_np.unique(coord0[:,i]) for i in range(3)] + coords = [_np.unique(coord0[:,i]) for i in range(3)] # _np.unique(coord0, axis=1) mincorner = _np.array(list(map(min,coords))) maxcorner = _np.array(list(map(max,coords))) grid = _np.array(list(map(len,coords)),'i') @@ -216,7 +219,7 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True): _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2])): raise ValueError('Regular grid spacing violated.') - if ordered and not _np.allclose(coord0.reshape(tuple(grid[::-1])+(3,)),cell_coord0(grid,size,origin)): + if ordered and not _np.allclose(coord0.reshape(tuple(grid)+(3,),order='F'),cell_coord0(grid,size,origin)): raise ValueError('Input data is not a regular grid.') return (grid,size,origin) @@ -241,17 +244,18 @@ def node_coord0(grid,size,origin=_np.zeros(3)): Parameters ---------- - grid : numpy.ndarray + grid : numpy.ndarray of shape (3) number of grid points. - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. - origin : numpy.ndarray, optional + origin : numpy.ndarray of shape (3), optional physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. """ - return _np.mgrid[origin[0]:size[0]+origin[0]:(grid[0]+1)*1j, - origin[1]:size[1]+origin[1]:(grid[1]+1)*1j, - origin[2]:size[2]+origin[2]:(grid[2]+1)*1j].T + return _np.stack(_np.meshgrid(_np.linspace(origin[0],size[0]+origin[0],grid[0]+1), + _np.linspace(origin[1],size[1]+origin[1],grid[1]+1), + _np.linspace(origin[2],size[2]+origin[2],grid[2]+1),indexing = 'ij'), + axis = -1) def node_displacement_fluct(size,F): @@ -260,7 +264,7 @@ def node_displacement_fluct(size,F): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. F : numpy.ndarray deformation gradient field. @@ -275,14 +279,14 @@ def node_displacement_avg(size,F): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. F : numpy.ndarray deformation gradient field. """ F_avg = _np.average(F,axis=(0,1,2)) - return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),node_coord0(F.shape[:3][::-1],size)) + return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),node_coord0(F.shape[:3],size)) def node_displacement(size,F): @@ -291,7 +295,7 @@ def node_displacement(size,F): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. F : numpy.ndarray deformation gradient field. @@ -306,15 +310,15 @@ def node_coord(size,F,origin=_np.zeros(3)): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size of the periodic field. F : numpy.ndarray deformation gradient field. - origin : numpy.ndarray, optional + origin : numpy.ndarray of shape (3), optional physical origin of the periodic field. Defaults to [0.0,0.0,0.0]. """ - return node_coord0(F.shape[:3][::-1],size,origin) + node_displacement(size,F) + return node_coord0(F.shape[:3],size,origin) + node_displacement(size,F) def cell_2_node(cell_data): @@ -335,19 +339,19 @@ def node_2_cell(node_data): return c[:-1,:-1,:-1] -def node_coord0_gridSizeOrigin(coord0,ordered=False): +def node_coord0_gridSizeOrigin(coord0,ordered=True): """ - Return grid 'DNA', i.e. grid, size, and origin from array of nodal positions. + Return grid 'DNA', i.e. grid, size, and origin from 1D array of nodal positions. Parameters ---------- - coord0 : numpy.ndarray - array of undeformed nodal coordinates. + coord0 : numpy.ndarray of shape (:,3) + undeformed nodal coordinates. ordered : bool, optional expect coord0 data to be ordered (x fast, z slow). """ - coords = [_np.unique(coord0[:,i]) for i in range(3)] + coords = [_np.unique(coord0[:,i]) for i in range(3)] # _np.unique(coord0, axis=1) mincorner = _np.array(list(map(min,coords))) maxcorner = _np.array(list(map(max,coords))) grid = _np.array(list(map(len,coords)),'i') - 1 @@ -362,7 +366,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=False): _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1)): raise ValueError('Regular grid spacing violated.') - if ordered and not _np.allclose(coord0.reshape(tuple((grid+1)[::-1])+(3,)),node_coord0(grid,size,origin)): + if ordered and not _np.allclose(coord0.reshape(tuple(grid+1)+(3,),order='F'),node_coord0(grid,size,origin)): raise ValueError('Input data is not a regular grid.') return (grid,size,origin) @@ -374,7 +378,7 @@ def regrid(size,F,new_grid): Parameters ---------- - size : numpy.ndarray + size : numpy.ndarray of shape (3) physical size F : numpy.ndarray deformation gradient field @@ -382,7 +386,7 @@ def regrid(size,F,new_grid): new grid for undeformed coordinates """ - c = cell_coord0(F.shape[:3][::-1],size) \ + c = cell_coord0(F.shape[:3],size) \ + cell_displacement_avg(size,F) \ + cell_displacement_fluct(size,F) diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index acbdbf688..3152072b7 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -4,18 +4,18 @@ import numpy as np from damask import grid_filters class TestGridFilters: - + def test_cell_coord0(self): size = np.random.random(3) grid = np.random.randint(8,32,(3)) coord = grid_filters.cell_coord0(grid,size) - assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid[::-1]) + (3,) + assert np.allclose(coord[0,0,0],size/grid*.5) and coord.shape == tuple(grid) + (3,) def test_node_coord0(self): size = np.random.random(3) grid = np.random.randint(8,32,(3)) coord = grid_filters.node_coord0(grid,size) - assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid[::-1]+1) + (3,) + assert np.allclose(coord[-1,-1,-1],size) and coord.shape == tuple(grid+1) + (3,) def test_coord0(self): size = np.random.random(3) @@ -31,7 +31,7 @@ class TestGridFilters: size = np.random.random(3) origin = np.random.random(3) coord0 = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode)) # noqa - _grid,_size,_origin = eval('grid_filters.{}_coord0_gridSizeOrigin(coord0.reshape(-1,3))'.format(mode)) + _grid,_size,_origin = eval('grid_filters.{}_coord0_gridSizeOrigin(coord0.reshape(-1,3,order="F"))'.format(mode)) assert np.allclose(grid,_grid) and np.allclose(size,_size) and np.allclose(origin,_origin) def test_displacement_fluct_equivalence(self): @@ -57,9 +57,9 @@ class TestGridFilters: shifted = eval('grid_filters.{}_coord0(grid,size,origin)'.format(mode)) unshifted = eval('grid_filters.{}_coord0(grid,size)'.format(mode)) if mode == 'cell': - assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]) +(3,))) + assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid) +(3,))) elif mode == 'node': - assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid[::-1]+1)+(3,))) + assert np.allclose(shifted,unshifted+np.broadcast_to(origin,tuple(grid+1)+(3,))) @pytest.mark.parametrize('function',[grid_filters.cell_displacement_avg, grid_filters.node_displacement_avg]) @@ -83,5 +83,5 @@ class TestGridFilters: def test_regrid(self): size = np.random.random(3) grid = np.random.randint(8,32,(3)) - F = np.broadcast_to(np.eye(3), tuple(grid[::-1])+(3,3)) + F = np.broadcast_to(np.eye(3), tuple(grid)+(3,3)) assert all(grid_filters.regrid(size,F,grid) == np.arange(grid.prod())) From 02dde3c2556a4e8c7ccf83a4b21670e9cc713083 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 13:10:13 +0200 Subject: [PATCH 25/47] test also invalid operations --- python/tests/test_grid_filters.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index 3152072b7..9eba75298 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -80,6 +80,11 @@ class TestGridFilters: F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) assert np.allclose(function(size,F),0.0) + def test_invalid_coordinates(self): + invalid_coordinates = np.random.random((np.random.randint(12,52),3)) + with pytest.raises(ValueError): + grid_filters.coord0_check(invalid_coordinates) + def test_regrid(self): size = np.random.random(3) grid = np.random.randint(8,32,(3)) From f519e62cd5713d2ba5b80ec574e695ec358969cf Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 13:56:33 +0200 Subject: [PATCH 26/47] testing data layout checks --- python/damask/grid_filters.py | 20 +++++++++--------- python/tests/test_grid_filters.py | 34 +++++++++++++++++++++++++++++-- 2 files changed, 42 insertions(+), 12 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 0c1c3867a..2ee99d276 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -196,7 +196,7 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True): expect coord0 data to be ordered (x fast, z slow). """ - coords = [_np.unique(coord0[:,i]) for i in range(3)] # _np.unique(coord0, axis=1) + coords = [_np.unique(coord0[:,i]) for i in range(3)] mincorner = _np.array(list(map(min,coords))) maxcorner = _np.array(list(map(max,coords))) grid = _np.array(list(map(len,coords)),'i') @@ -214,13 +214,13 @@ def cell_coord0_gridSizeOrigin(coord0,ordered=True): start = origin + delta*.5 end = origin - delta*.5 + size - if not _np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0])) and \ - _np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1])) and \ - _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2])): + if not (_np.allclose(coords[0],_np.linspace(start[0],end[0],grid[0])) and \ + _np.allclose(coords[1],_np.linspace(start[1],end[1],grid[1])) and \ + _np.allclose(coords[2],_np.linspace(start[2],end[2],grid[2]))): raise ValueError('Regular grid spacing violated.') if ordered and not _np.allclose(coord0.reshape(tuple(grid)+(3,),order='F'),cell_coord0(grid,size,origin)): - raise ValueError('Input data is not a regular grid.') + raise ValueError('Input data is not ordered (x fast, z slow).') return (grid,size,origin) @@ -351,7 +351,7 @@ def node_coord0_gridSizeOrigin(coord0,ordered=True): expect coord0 data to be ordered (x fast, z slow). """ - coords = [_np.unique(coord0[:,i]) for i in range(3)] # _np.unique(coord0, axis=1) + coords = [_np.unique(coord0[:,i]) for i in range(3)] mincorner = _np.array(list(map(min,coords))) maxcorner = _np.array(list(map(max,coords))) grid = _np.array(list(map(len,coords)),'i') - 1 @@ -361,13 +361,13 @@ def node_coord0_gridSizeOrigin(coord0,ordered=True): if (grid+1).prod() != len(coord0): raise ValueError('Data count {} does not match grid {}.'.format(len(coord0),grid)) - if not _np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \ - _np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \ - _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1)): + if not (_np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],grid[0]+1)) and \ + _np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],grid[1]+1)) and \ + _np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],grid[2]+1))): raise ValueError('Regular grid spacing violated.') if ordered and not _np.allclose(coord0.reshape(tuple(grid+1)+(3,),order='F'),node_coord0(grid,size,origin)): - raise ValueError('Input data is not a regular grid.') + raise ValueError('Input data is not ordered (x fast, z slow).') return (grid,size,origin) diff --git a/python/tests/test_grid_filters.py b/python/tests/test_grid_filters.py index 9eba75298..eb359006a 100644 --- a/python/tests/test_grid_filters.py +++ b/python/tests/test_grid_filters.py @@ -80,10 +80,40 @@ class TestGridFilters: F = np.broadcast_to(np.random.random((3,3)), tuple(grid)+(3,3)) assert np.allclose(function(size,F),0.0) - def test_invalid_coordinates(self): + @pytest.mark.parametrize('function',[grid_filters.coord0_check, + grid_filters.node_coord0_gridSizeOrigin, + grid_filters.cell_coord0_gridSizeOrigin]) + def test_invalid_coordinates(self,function): invalid_coordinates = np.random.random((np.random.randint(12,52),3)) with pytest.raises(ValueError): - grid_filters.coord0_check(invalid_coordinates) + function(invalid_coordinates) + + @pytest.mark.parametrize('function',[grid_filters.node_coord0_gridSizeOrigin, + grid_filters.cell_coord0_gridSizeOrigin]) + def test_uneven_spaced_coordinates(self,function): + start = np.random.random(3) + end = np.random.random(3)*10. + start + grid = np.random.randint(8,32,(3)) + uneven = np.stack(np.meshgrid(np.logspace(start[0],end[0],grid[0]), + np.logspace(start[1],end[1],grid[1]), + np.logspace(start[2],end[2],grid[2]),indexing = 'ij'), + axis = -1).reshape((grid.prod(),3),order='F') + with pytest.raises(ValueError): + function(uneven) + + @pytest.mark.parametrize('mode',[True,False]) + @pytest.mark.parametrize('function',[grid_filters.node_coord0_gridSizeOrigin, + grid_filters.cell_coord0_gridSizeOrigin]) + def test_unordered_coordinates(self,function,mode): + origin = np.random.random(3) + size = np.random.random(3)*10.+origin + grid = np.random.randint(8,32,(3)) + unordered = grid_filters.node_coord0(grid,size,origin).reshape(-1,3) + if mode: + with pytest.raises(ValueError): + function(unordered,mode) + else: + function(unordered,mode) def test_regrid(self): size = np.random.random(3) From a367128436fb099c3fe4f8b920f61b31e169f1f4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 13:57:33 +0200 Subject: [PATCH 27/47] not used anymore --- processing/post/addDisplacement.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/processing/post/addDisplacement.py b/processing/post/addDisplacement.py index a9424f5e1..6f4a60192 100755 --- a/processing/post/addDisplacement.py +++ b/processing/post/addDisplacement.py @@ -5,8 +5,6 @@ import sys from io import StringIO from optparse import OptionParser -import numpy as np - import damask From 8465c3ab1b7b19a84b4aa00c56b7aada29131ea1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 19:13:45 +0200 Subject: [PATCH 28/47] adjusting to new data layout --- processing/pre/geom_fromVoronoiTessellation.py | 2 +- processing/pre/seeds_fromGeom.py | 8 ++++---- processing/pre/seeds_fromRandom.py | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py index df40176d8..2b4da133b 100755 --- a/processing/pre/geom_fromVoronoiTessellation.py +++ b/processing/pre/geom_fromVoronoiTessellation.py @@ -224,7 +224,7 @@ for name in filenames: header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\ + config_header - geom = damask.Geom(indices.reshape(grid),size,origin, + geom = damask.Geom(indices.reshape(grid,order='F'),size,origin, homogenization=options.homogenization,comments=header) damask.util.croak(geom) diff --git a/processing/pre/seeds_fromGeom.py b/processing/pre/seeds_fromGeom.py index 80dc0d6f5..0b741a077 100755 --- a/processing/pre/seeds_fromGeom.py +++ b/processing/pre/seeds_fromGeom.py @@ -45,7 +45,7 @@ options.blacklist = [int(i) for i in options.blacklist] for name in filenames: damask.util.report(scriptName,name) - + geom = damask.Geom.from_file(StringIO(''.join(sys.stdin.read())) if name is None else name) microstructure = geom.get_microstructure().reshape((-1,1),order='F') @@ -53,9 +53,9 @@ for name in filenames: np.full(geom.grid.prod(),True,dtype=bool), np.in1d(microstructure,options.blacklist,invert=True) if options.blacklist else \ np.full(geom.grid.prod(),True,dtype=bool)) - - seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3) - + + seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F') + comments = geom.comments \ + [scriptID + ' ' + ' '.join(sys.argv[1:]), 'grid\ta {}\tb {}\tc {}'.format(*geom.grid), diff --git a/processing/pre/seeds_fromRandom.py b/processing/pre/seeds_fromRandom.py index 2de513c2c..be690713c 100755 --- a/processing/pre/seeds_fromRandom.py +++ b/processing/pre/seeds_fromRandom.py @@ -128,7 +128,7 @@ for name in filenames: if not options.selective: - coords = damask.grid_filters.cell_coord0(grid,size).reshape(-1,3) + coords = damask.grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') seeds = coords[np.random.choice(np.prod(grid), options.N, replace=False)] \ + np.broadcast_to(size/grid,(options.N,3))*(np.random.rand(options.N,3)*.5-.25) # wobble without leaving grid else: From 4db91fff03ab68f68cf26ac2c5311fd10c82b8c8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 20:14:27 +0200 Subject: [PATCH 29/47] simplified --- python/damask/_geom.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index a33e855d8..f01277f95 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -322,11 +322,10 @@ class Geom: if i != grid.prod(): raise TypeError('Invalid file: expected {} entries, found {}'.format(grid.prod(),i)) - microstructure = microstructure.reshape(grid,order='F') - if not np.any(np.mod(microstructure.flatten(),1) != 0.0): # no float present + if not np.any(np.mod(microstructure,1) != 0.0): # no float present microstructure = microstructure.astype('int') - return Geom(microstructure.reshape(grid),size,origin,homogenization,comments) + return Geom(microstructure.reshape(grid,order='F'),size,origin,homogenization,comments) @staticmethod From 1ccda22ae354de44c78da730543ccd7b747cb906 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 20:16:25 +0200 Subject: [PATCH 30/47] nothing to flatten here --- python/damask/_geom.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index f01277f95..8010cd167 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -351,13 +351,13 @@ class Geom: """ if periodic: - weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3) + weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3) seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]]))) coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F') else: - weights_p = weights.flatten() + weights_p = weights seeds_p = seeds coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') From 25d29bf438d00390f105bbdff29c6c85ddc3472d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 20:24:55 +0200 Subject: [PATCH 31/47] no need for Fortran order forward/backward --- python/damask/_geom.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 8010cd167..892000b7c 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -355,11 +355,11 @@ class Geom: seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]]))) - coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3,order='F') + coords = grid_filters.cell_coord0(grid*3,size*3,-size).reshape(-1,3) else: weights_p = weights seeds_p = seeds - coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') + coords = grid_filters.cell_coord0(grid,size).reshape(-1,3) pool = multiprocessing.Pool(processes = int(Environment().options['DAMASK_NUM_THREADS'])) result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords]) @@ -368,10 +368,10 @@ class Geom: microstructure = np.array(result.get()) if periodic: - microstructure = microstructure.reshape(grid*3,order='F') + microstructure = microstructure.reshape(grid*3) microstructure = microstructure[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] else: - microstructure = microstructure.reshape(grid,order='F') + microstructure = microstructure.reshape(grid) #comments = 'geom.py:from_Laguerre_tessellation v{}'.format(version) return Geom(microstructure+1,size,homogenization=1) @@ -394,12 +394,12 @@ class Geom: perform a periodic tessellation. Defaults to True. """ - coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F') + coords = grid_filters.cell_coord0(grid,size).reshape(-1,3) KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) devNull,microstructure = KDTree.query(coords) #comments = 'geom.py:from_Voronoi_tessellation v{}'.format(version) - return Geom(microstructure.reshape(grid,order='F')+1,size,homogenization=1) + return Geom(microstructure.reshape(grid)+1,size,homogenization=1) def to_file(self,fname,pack=None): From a8b75a23c5f746fec18fd771db2c02dbe30fb320 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 20:40:47 +0200 Subject: [PATCH 32/47] simplified --- processing/pre/geom_fromVoronoiTessellation.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/processing/pre/geom_fromVoronoiTessellation.py b/processing/pre/geom_fromVoronoiTessellation.py index 2b4da133b..e6755c28e 100755 --- a/processing/pre/geom_fromVoronoiTessellation.py +++ b/processing/pre/geom_fromVoronoiTessellation.py @@ -24,22 +24,22 @@ def findClosestSeed(seeds, weights, point): def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), periodic = True, cpus = 2): if periodic: - weights_p = np.tile(weights,27).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3) + weights_p = np.tile(weights.squeeze(),27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3) seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.]))) seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]]))) - coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3,order='F') + coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3) else: - weights_p = weights.flatten() + weights_p = weights.squeeze() seeds_p = seeds - coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F') + coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3) if cpus > 1: pool = multiprocessing.Pool(processes = cpus) result = pool.map_async(partial(findClosestSeed,seeds_p,weights_p), [coord for coord in coords]) pool.close() pool.join() - closest_seed = np.array(result.get()) + closest_seed = np.array(result.get()).reshape(-1,3) else: closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords]) @@ -52,7 +52,7 @@ def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), peri def Voronoi_tessellation(grid, size, seeds, origin = np.zeros(3), periodic = True): - coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3,order='F') + coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3) KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds) devNull,closest_seed = KDTree.query(coords) @@ -224,7 +224,7 @@ for name in filenames: header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\ + config_header - geom = damask.Geom(indices.reshape(grid,order='F'),size,origin, + geom = damask.Geom(indices.reshape(grid),size,origin, homogenization=options.homogenization,comments=header) damask.util.croak(geom) From 9a54c326e23a5e876cd18febb154452591f7904a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 21:28:23 +0200 Subject: [PATCH 33/47] adopted to new data layout --- processing/post/averageDown.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/processing/post/averageDown.py b/processing/post/averageDown.py index 39d925d19..341cc748d 100755 --- a/processing/post/averageDown.py +++ b/processing/post/averageDown.py @@ -61,7 +61,7 @@ if any(shift != 0): prefix += 'shift{:+}{:+}{:+}_'.format(*shift) for name in filenames: damask.util.report(scriptName,name) - + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) if (options.grid is None or options.size is None): @@ -87,11 +87,11 @@ for name in filenames: origin = list(-(packing//2)) + [0])\ [::packing[0],::packing[1],::packing[2],:].reshape((packedGrid.prod(),-1),order = 'F') - + table = damask.Table(averagedDown,table.shapes,table.comments) coords = damask.grid_filters.cell_coord0(packedGrid,size,shift/packedGrid*size+origin) - table.set(options.pos, coords.reshape(-1,3)) + table.set(options.pos, coords.reshape(-1,3,order='F')) outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name)) From 4278ba32ae2ddd71a4522bfe390e0f7efec518a8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 21:44:38 +0200 Subject: [PATCH 34/47] adopting to standard data layout --- processing/post/addCompatibilityMismatch.py | 55 ++++++++++----------- 1 file changed, 27 insertions(+), 28 deletions(-) diff --git a/processing/post/addCompatibilityMismatch.py b/processing/post/addCompatibilityMismatch.py index 0678e78ba..7cab9c2ca 100755 --- a/processing/post/addCompatibilityMismatch.py +++ b/processing/post/addCompatibilityMismatch.py @@ -16,8 +16,8 @@ scriptID = ' '.join([scriptName,damask.version]) def volTetrahedron(coords): """ Return the volume of the tetrahedron with given vertices or sides. - - Ifvertices are given they must be in a NumPy array with shape (4,3): the + + If vertices are given they must be in a NumPy array with shape (4,3): the position vectors of the 4 vertices in 3 dimensions; if the six sides are given, they must be an array of length 6. If both are given, the sides will be used in the calculation. @@ -62,19 +62,18 @@ def volTetrahedron(coords): def volumeMismatch(size,F,nodes): """ Calculates the volume mismatch. - - volume mismatch is defined as the difference between volume of reconstructed + + volume mismatch is defined as the difference between volume of reconstructed (compatible) cube and determinant of deformation gradient at Fourier point. """ coords = np.empty([8,3]) - vMismatch = np.empty(grid[::-1]) - volInitial = size.prod()/grid.prod() - + vMismatch = np.empty(F.shape[:3]) + #-------------------------------------------------------------------------------------------------- # calculate actual volume and volume resulting from deformation gradient - for k in range(grid[2]): + for k in range(grid[0]): for j in range(grid[1]): - for i in range(grid[0]): + for i in range(grid[2]): coords[0,0:3] = nodes[k, j, i ,0:3] coords[1,0:3] = nodes[k ,j, i+1,0:3] coords[2,0:3] = nodes[k ,j+1,i+1,0:3] @@ -91,21 +90,21 @@ def volumeMismatch(size,F,nodes): + abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[5,0:3]])) \ + abs(volTetrahedron([coords[6,0:3],coords[4,0:3],coords[1,0:3],coords[0,0:3]]))) \ /np.linalg.det(F[k,j,i,0:3,0:3]) - return vMismatch/volInitial + return vMismatch/(size.prod()/grid.prod()) def shapeMismatch(size,F,nodes,centres): """ Routine to calculate the shape mismatch. - + shape mismatch is defined as difference between the vectors from the central point to the corners of reconstructed (combatible) volume element and the vectors calculated by deforming the initial volume element with the current deformation gradient. """ coordsInitial = np.empty([8,3]) - sMismatch = np.empty(grid[::-1]) - + sMismatch = np.empty(F.shape[:3]) + #-------------------------------------------------------------------------------------------------- # initial positions coordsInitial[0,0:3] = [-size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]] @@ -117,21 +116,21 @@ def shapeMismatch(size,F,nodes,centres): coordsInitial[6,0:3] = [+size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]] coordsInitial[7,0:3] = [-size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]] coordsInitial = coordsInitial/2.0 - + #-------------------------------------------------------------------------------------------------- # compare deformed original and deformed positions to actual positions - for k in range(grid[2]): + for k in range(grid[0]): for j in range(grid[1]): - for i in range(grid[0]): + for i in range(grid[2]): sMismatch[k,j,i] = \ + np.linalg.norm(nodes[k, j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[0,0:3]))\ - + np.linalg.norm(nodes[k, j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\ - + np.linalg.norm(nodes[k, j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\ + + np.linalg.norm(nodes[k+1,j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[1,0:3]))\ + + np.linalg.norm(nodes[k+1,j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[2,0:3]))\ + np.linalg.norm(nodes[k, j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[3,0:3]))\ - + np.linalg.norm(nodes[k+1,j, i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\ + + np.linalg.norm(nodes[k, j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[4,0:3]))\ + np.linalg.norm(nodes[k+1,j, i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[5,0:3]))\ + np.linalg.norm(nodes[k+1,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[6,0:3]))\ - + np.linalg.norm(nodes[k+1,j+1,i ,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3])) + + np.linalg.norm(nodes[k ,j+1,i+1,0:3] - centres[k,j,i,0:3] - np.dot(F[k,j,i,:,:], coordsInitial[7,0:3])) return sMismatch @@ -174,24 +173,24 @@ if filenames == []: filenames = [None] for name in filenames: damask.util.report(scriptName,name) - + table = damask.Table.from_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name) grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos)) - F = table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3) + F = table.get(options.defgrad).reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+(3,3)) nodes = damask.grid_filters.node_coord(size,F) - + if options.shape: centers = damask.grid_filters.cell_coord(size,F) - shapeMismatch = shapeMismatch( size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes,centers) + shapeMismatch = shapeMismatch(size,F,nodes,centers) table.add('shapeMismatch(({}))'.format(options.defgrad), - shapeMismatch.reshape(-1,1), + shapeMismatch.reshape(-1,1,order='F'), scriptID+' '+' '.join(sys.argv[1:])) - + if options.volume: - volumeMismatch = volumeMismatch(size,table.get(options.defgrad).reshape(grid[2],grid[1],grid[0],3,3),nodes) + volumeMismatch = volumeMismatch(size,F,nodes) table.add('volMismatch(({}))'.format(options.defgrad), - volumeMismatch.reshape(-1,1), + volumeMismatch.reshape(-1,1,order='F'), scriptID+' '+' '.join(sys.argv[1:])) table.to_ASCII(sys.stdout if name is None else name) From 85f831d5b44594db2c6460e693715a0d92d77164 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 22:14:57 +0200 Subject: [PATCH 35/47] simplified --- processing/post/addCompatibilityMismatch.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/processing/post/addCompatibilityMismatch.py b/processing/post/addCompatibilityMismatch.py index 7cab9c2ca..59eafc077 100755 --- a/processing/post/addCompatibilityMismatch.py +++ b/processing/post/addCompatibilityMismatch.py @@ -93,7 +93,6 @@ def volumeMismatch(size,F,nodes): return vMismatch/(size.prod()/grid.prod()) - def shapeMismatch(size,F,nodes,centres): """ Routine to calculate the shape mismatch. @@ -102,20 +101,19 @@ def shapeMismatch(size,F,nodes,centres): the corners of reconstructed (combatible) volume element and the vectors calculated by deforming the initial volume element with the current deformation gradient. """ - coordsInitial = np.empty([8,3]) sMismatch = np.empty(F.shape[:3]) #-------------------------------------------------------------------------------------------------- # initial positions - coordsInitial[0,0:3] = [-size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]] - coordsInitial[1,0:3] = [+size[0]/grid[0],-size[1]/grid[1],-size[2]/grid[2]] - coordsInitial[2,0:3] = [+size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]] - coordsInitial[3,0:3] = [-size[0]/grid[0],+size[1]/grid[1],-size[2]/grid[2]] - coordsInitial[4,0:3] = [-size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]] - coordsInitial[5,0:3] = [+size[0]/grid[0],-size[1]/grid[1],+size[2]/grid[2]] - coordsInitial[6,0:3] = [+size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]] - coordsInitial[7,0:3] = [-size[0]/grid[0],+size[1]/grid[1],+size[2]/grid[2]] - coordsInitial = coordsInitial/2.0 + delta = size/grid*.5 + coordsInitial = np.vstack((delta * np.array((-1,-1,-1)), + delta * np.array((+1,-1,-1)), + delta * np.array((+1,+1,-1)), + delta * np.array((-1,+1,-1)), + delta * np.array((-1,-1,+1)), + delta * np.array((+1,-1,+1)), + delta * np.array((+1,+1,+1)), + delta * np.array((-1,+1,+1)))) #-------------------------------------------------------------------------------------------------- # compare deformed original and deformed positions to actual positions From 93db1b37cacc4ccf6b0ac3d3fbaa7660f19d76c5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 22:26:29 +0200 Subject: [PATCH 36/47] some useful information --- python/damask/grid_filters.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 2ee99d276..3925c8db3 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -1,3 +1,17 @@ +""" +Filters for operations on regular grids. + +Notes +----- +The grids are defined as (x,y,z,...) where x is fastest and z is slowest. +This convention is consistent with the geom file format. +When converting to/from a plain list (e.g. storage in ASCII table), +the following operations are required for tensorial: + +D3 = D1.reshape(grid+(-1,),order='F').reshape(grid+(3,3)) +D1 = D3.reshape(grid+(-1,)).reshape(-1,9,order='F') + +""" from scipy import spatial as _spatial import numpy as _np From 00cb17ccfaf867a6f7e4f321a7153f89eec34098 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 22:47:19 +0200 Subject: [PATCH 37/47] commenting --- python/damask/grid_filters.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 3925c8db3..f74bc92ac 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -6,7 +6,7 @@ Notes The grids are defined as (x,y,z,...) where x is fastest and z is slowest. This convention is consistent with the geom file format. When converting to/from a plain list (e.g. storage in ASCII table), -the following operations are required for tensorial: +the following operations are required for tensorial data: D3 = D1.reshape(grid+(-1,),order='F').reshape(grid+(3,3)) D1 = D3.reshape(grid+(-1,)).reshape(-1,9,order='F') @@ -23,6 +23,10 @@ def _ks(size,grid,first_order=False): ---------- size : numpy.ndarray of shape (3) physical size of the periodic field. + grid : numpy.ndarray of shape (3) + number of grid points. + first_order : bool, optional + correction for first order derivatives, defaults to False. """ k_sk = _np.where(_np.arange(grid[0])>grid[0]//2,_np.arange(grid[0])-grid[0],_np.arange(grid[0]))/size[0] @@ -44,6 +48,8 @@ def curl(size,field): ---------- size : numpy.ndarray of shape (3) physical size of the periodic field. + field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3) + periodic field of which the curl is calculated. """ n = _np.prod(field.shape[3:]) @@ -68,6 +74,8 @@ def divergence(size,field): ---------- size : numpy.ndarray of shape (3) physical size of the periodic field. + field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3) + periodic field of which the divergence is calculated. """ n = _np.prod(field.shape[3:]) @@ -82,12 +90,14 @@ def divergence(size,field): def gradient(size,field): """ - Calculate gradient of a vector or scalar field in Fourier space. + Calculate gradient of a scalar or tensor field in Fourier space. Parameters ---------- size : numpy.ndarray of shape (3) physical size of the periodic field. + field : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3) + periodic field of which the gradient is calculated. """ n = _np.prod(field.shape[3:]) @@ -394,9 +404,9 @@ def regrid(size,F,new_grid): ---------- size : numpy.ndarray of shape (3) physical size - F : numpy.ndarray + F : numpy.ndarray of shape (:,:,:,3,3) deformation gradient field - new_grid : numpy.ndarray + new_grid : numpy.ndarray of shape (3) new grid for undeformed coordinates """ From 528378eba22536d7b377db1de9f47d68923f51ce Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 22:51:51 +0200 Subject: [PATCH 38/47] polishing/bugfix * HDF5 data has x fast, therefore reshape needs to be aware of 'F' order * Datatype for True/False is 'bool' in python, not 'Boolean' --- python/damask/_result.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index ee7d696a9..65ac7a01b 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -105,7 +105,7 @@ class Result: select from 'set', 'add', and 'del' what : str attribute to change (must be from self.selection) - datasets : list of str or Boolean + datasets : list of str or bool name of datasets as list, supports ? and * wildcards. True is equivalent to [*], False is equivalent to [] @@ -197,7 +197,7 @@ class Result: ---------- what : str attribute to change (must be from self.selection) - datasets : list of str or Boolean + datasets : list of str or bool name of datasets as list, supports ? and * wildcards. True is equivalent to [*], False is equivalent to [] @@ -213,7 +213,7 @@ class Result: ---------- what : str attribute to change (must be from self.selection) - datasets : list of str or Boolean + datasets : list of str or bool name of datasets as list, supports ? and * wildcards. True is equivalent to [*], False is equivalent to [] @@ -229,7 +229,7 @@ class Result: ---------- what : str attribute to change (must be from self.selection) - datasets : list of str or Boolean + datasets : list of str or bool name of datasets as list, supports ? and * wildcards. True is equivalent to [*], False is equivalent to [] @@ -256,10 +256,10 @@ class Result: datasets : iterable or str component : int homogenization component to consider for constituent data - tagged : Boolean + tagged : bool tag Table.column name with '#component' defaults to False - split : Boolean + split : bool split Table by increment and return dictionary of Tables defaults to True @@ -320,7 +320,7 @@ class Result: Parameters ---------- - datasets : iterable or str or Boolean + datasets : iterable or str or bool Examples -------- @@ -454,7 +454,7 @@ class Result: def cell_coordinates(self): """Return initial coordinates of the cell centers.""" if self.structured: - return grid_filters.cell_coord0(self.grid,self.size,self.origin).reshape(-1,3) + return grid_filters.cell_coord0(self.grid,self.size,self.origin).reshape(-1,3,order='F') else: with h5py.File(self.fname,'r') as f: return f['geometry/x_c'][()] From 7084fc9616d5833fc554faf471c2d2d3bb3754d4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 20 Apr 2020 23:19:38 +0200 Subject: [PATCH 39/47] adjusting to consistent data layout --- processing/post/blowUp.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/processing/post/blowUp.py b/processing/post/blowUp.py index 316b74753..23c6b2ef2 100755 --- a/processing/post/blowUp.py +++ b/processing/post/blowUp.py @@ -59,13 +59,13 @@ for name in filenames: packing = np.array(options.packing,'i') outSize = grid*packing - data = table.data.values.reshape(tuple(grid)+(-1,)) - blownUp = ndimage.interpolation.zoom(data,tuple(packing)+(1,),order=0,mode='nearest').reshape(outSize.prod(),-1) + data = table.data.values.reshape(tuple(grid)+(-1,),order='F') + blownUp = ndimage.interpolation.zoom(data,tuple(packing)+(1,),order=0,mode='nearest').reshape(outSize.prod(),-1,order='F') table = damask.Table(blownUp,table.shapes,table.comments) coords = damask.grid_filters.cell_coord0(outSize,size,origin) - table.set(options.pos,coords.reshape(-1,3)) + table.set(options.pos,coords.reshape(-1,3,order='F')) table.set('elem',np.arange(1,outSize.prod()+1)) outname = os.path.join(os.path.dirname(name),prefix+os.path.basename(name)) From 4a05cdbbbb0b6fa977ecf083c6ca8c0bddff1e93 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Apr 2020 03:26:26 +0200 Subject: [PATCH 40/47] copy and paste error --- python/damask/grid_filters.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index f74bc92ac..eaee42924 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -90,7 +90,7 @@ def divergence(size,field): def gradient(size,field): """ - Calculate gradient of a scalar or tensor field in Fourier space. + Calculate gradient of a scalar or vector field in Fourier space. Parameters ---------- From b26bc78a44f42e9507e017c9f1d8aa7cb5a61c74 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Apr 2020 18:49:50 +0200 Subject: [PATCH 41/47] adopting to new grid layout --- processing/post/DADF5_postResults.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/DADF5_postResults.py b/processing/post/DADF5_postResults.py index 70c883aeb..2230ba563 100755 --- a/processing/post/DADF5_postResults.py +++ b/processing/post/DADF5_postResults.py @@ -33,7 +33,7 @@ for filename in options.filenames: results = damask.Result(filename) if not results.structured: continue - coords = damask.grid_filters.cell_coord0(results.grid,results.size,results.origin) + coords = damask.grid_filters.cell_coord0(results.grid,results.size,results.origin).reshape(-1,3,order='F') N_digits = int(np.floor(np.log10(int(results.increments[-1][3:]))))+1 N_digits = 5 # hack to keep test intact From 7452ebc727af4d5e71f46036b4d6166aec034fb2 Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 21 Apr 2020 19:37:07 +0200 Subject: [PATCH 42/47] [skip ci] updated version information after successful test of v2.0.3-2311-gf1afc159 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 7be7a7f6c..7cb1691ed 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2303-g2a6132b7 +v2.0.3-2311-gf1afc159 From 936ce6a160020e2b038805e8091757665362ce8e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 22 Apr 2020 18:26:10 +0200 Subject: [PATCH 43/47] reshape for scalar/vector not needed --- processing/post/addGradient.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/addGradient.py b/processing/post/addGradient.py index c6df0eacf..f63f24789 100755 --- a/processing/post/addGradient.py +++ b/processing/post/addGradient.py @@ -49,7 +49,7 @@ for name in filenames: for label in options.labels: field = table.get(label) shape = (1,) if np.prod(field.shape)//np.prod(grid) == 1 else (3,) # scalar or vector - field = field.reshape(tuple(grid)+(-1,),order='F').reshape(tuple(grid)+shape) + field = field.reshape(tuple(grid)+(-1,),order='F') grad = damask.grid_filters.gradient(size,field) table.add('gradFFT({})'.format(label), grad.reshape(tuple(grid)+(-1,)).reshape(-1,np.prod(shape)*3,order='F'), From f0e2437d751b23f98f7cfa8cceb8805c7cf5b2e0 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 22 Apr 2020 21:36:05 -0400 Subject: [PATCH 44/47] [skip ci] added "list(map..." required since Python3 --- processing/post/permuteData.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/processing/post/permuteData.py b/processing/post/permuteData.py index 81af71adb..72759f2ce 100755 --- a/processing/post/permuteData.py +++ b/processing/post/permuteData.py @@ -60,7 +60,7 @@ for name in filenames: table.head_read() -# ------------------------------------------ process labels --------------------------------------- +# ------------------------------------------ process labels --------------------------------------- errors = [] remarks = [] @@ -80,7 +80,7 @@ for name in filenames: damask.util.croak(errors) table.close(dismiss = True) continue - + # ------------------------------------------ assemble header --------------------------------------- randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed # random seed per file @@ -97,17 +97,17 @@ for name in filenames: for col,dim in zip(columns,dims): if options.unique: s = set(map(tuple,table.data[:,col:col+dim])) # generate set of (unique) values - uniques = np.array(map(np.array,s)) # translate set to np.array + uniques = np.array(list(map(np.array,s))) # translate set to np.array shuffler = dict(zip(s,np.random.permutation(len(s)))) # random permutation - table.data[:,col:col+dim] = uniques[np.array(map(lambda x: shuffler[tuple(x)], - table.data[:,col:col+dim]))] # fill table with mapped uniques + table.data[:,col:col+dim] = uniques[np.array(list(map(lambda x: shuffler[tuple(x)], + table.data[:,col:col+dim])))] # fill table with mapped uniques else: np.random.shuffle(table.data[:,col:col+dim]) # independently shuffle every row -# ------------------------------------------ output result ----------------------------------------- +# ------------------------------------------ output result ----------------------------------------- table.data_writeArray() -# ------------------------------------------ output finalization ----------------------------------- +# ------------------------------------------ output finalization ----------------------------------- table.close() # close ASCII tables From b2f8e8847fc6395e03fab058a6be6ad7f530b44e Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 23 Apr 2020 21:34:11 +0200 Subject: [PATCH 45/47] [skip ci] updated version information after successful test of v2.0.3-2339-g28980887 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 7cb1691ed..10a6e2bdb 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2311-gf1afc159 +v2.0.3-2339-g28980887 From 62f7363a29cac549a5853185ede1dfd78c27cf45 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 23 Apr 2020 22:11:57 +0200 Subject: [PATCH 46/47] orientation conversion test in PRIVATE has the tendency to fail new tolerances lead more often to a transfer from phi_2 to phi_1 for PHI=180 or PHI=0. The pytest takes care of that, the deprecated test in PRIVATE not --- .gitlab-ci.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index dd20c8baa..501b9e193 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -203,7 +203,6 @@ Post_OrientationConversion: stage: postprocessing script: - OrientationConversion/test.py - - OrientationConversion/test2.py except: - master - release From 05fd943cf480868d6a1d486281808185ed012bba Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 24 Apr 2020 00:45:37 +0200 Subject: [PATCH 47/47] [skip ci] updated version information after successful test of v2.0.3-2364-g62f7363a --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 10a6e2bdb..6e24519bd 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-2339-g28980887 +v2.0.3-2364-g62f7363a