From b58465415b01ba44379bad7e7d2ee15bc1fc580e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 18:20:47 +0100 Subject: [PATCH 01/12] store damage parameter like temperature --- src/grid/grid_damage_spectral.f90 | 2 -- src/homogenization.f90 | 26 ++++++++++++---- src/homogenization_damage.f90 | 49 ++++++++++++++++++++++++++++--- 3 files changed, 66 insertions(+), 11 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index a891f070d..dc78d9e44 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -198,7 +198,6 @@ function grid_damage_spectral_solution(timeinc) result(solution) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) - homogenization_phi(cell) = phi_current(i,j,k) enddo; enddo; enddo call VecMin(solution_vec,devNull,phi_min,ierr); CHKERRQ(ierr) @@ -236,7 +235,6 @@ subroutine grid_damage_spectral_forward(cutBack) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 call damage_nonlocal_putNonLocalDamage(phi_current(i,j,k),1,cell) - homogenization_phi(cell) = phi_current(i,j,k) enddo; enddo; enddo else phi_lastInc = phi_current diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 64c2aade7..048bf46df 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -25,9 +25,6 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point - real(pReal), dimension(:), allocatable, public :: & - homogenization_phi, & - homogenization_dot_phi real(pReal), dimension(:,:,:), allocatable, public :: & homogenization_F0, & !< def grad of IP at start of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment @@ -75,8 +72,7 @@ module homogenization integer, intent(in) :: ce end subroutine thermal_partition - module subroutine damage_partition(phi,ce) - real(pReal), intent(in) :: phi + module subroutine damage_partition(ce) integer, intent(in) :: ce end subroutine damage_partition @@ -330,6 +326,26 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE enddo !$OMP END DO + !$OMP DO PRIVATE(ho,ph,ce) + do el = FEsolving_execElem(1),FEsolving_execElem(2) + if (terminallyIll) continue + ho = material_homogenizationAt(el) + do ip = FEsolving_execIP(1),FEsolving_execIP(2) + ce = (el-1)*discretization_nIPs + ip + call damage_partition(ce) + ! do co = 1, homogenization_Nconstituents(ho) + ! ph = material_phaseAt(co,el) + ! if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then + ! if (.not. terminallyIll) & ! so first signals terminally ill... + ! print*, ' Integration point ', ip,' at element ', el, ' terminally ill' + ! terminallyIll = .true. ! ...and kills all others + ! endif + ! call thermal_homogenize(ip,el) + ! enddo + enddo + enddo + !$OMP END DO + !$OMP DO PRIVATE(ho) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 479318340..ad029cce8 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -3,19 +3,58 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_damage + type :: tDataContainer + real(pReal), dimension(:), allocatable :: phi + end type tDataContainer + + type(tDataContainer), dimension(:), allocatable :: current + + type :: tParameters + character(len=pStringLen), allocatable, dimension(:) :: & + output + end type tParameters + + type(tparameters), dimension(:), allocatable :: & + param contains + !-------------------------------------------------------------------------------------------------- !> @brief Allocate variables and set parameters. !-------------------------------------------------------------------------------------------------- module subroutine damage_init() + class(tNode), pointer :: & + configHomogenizations, & + configHomogenization, & + configHomogenizationDamage + integer :: ho + print'(/,a)', ' <<<+- homogenization_damage init -+>>>' - allocate(homogenization_phi(discretization_nIPs*discretization_Nelems)) - allocate(homogenization_dot_phi(discretization_nIPs*discretization_Nelems)) + + configHomogenizations => config_material%get('homogenization') + allocate(param(configHomogenizations%length)) + allocate(current(configHomogenizations%length)) + + do ho = 1, configHomogenizations%length + allocate(current(ho)%phi(count(material_homogenizationAt2==ho)), source=1.0_pReal) + configHomogenization => configHomogenizations%get(ho) + associate(prm => param(ho)) + if (configHomogenization%contains('damage')) then + configHomogenizationDamage => configHomogenization%get('damage') +#if defined (__GFORTRAN__) + prm%output = output_asStrings(configHomogenizationDamage) +#else + prm%output = configHomogenizationDamage%get_asStrings('output',defaultVal=emptyStringArray) +#endif + else + prm%output = emptyStringArray + endif + end associate + enddo end subroutine damage_init @@ -23,13 +62,15 @@ end subroutine damage_init !-------------------------------------------------------------------------------------------------- !> @brief Partition temperature onto the individual constituents. !-------------------------------------------------------------------------------------------------- -module subroutine damage_partition(phi,ce) +module subroutine damage_partition(ce) - real(pReal), intent(in) :: phi + real(pReal) :: phi integer, intent(in) :: ce integer :: co + + phi = current(material_homogenizationAt2(ce))%phi(material_homogenizationMemberAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) call constitutive_damage_set_phi(phi,co,ce) enddo From 1f94a64ca6b1bee8396e7d347de27bc90cc16f88 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 18:47:19 +0100 Subject: [PATCH 02/12] part of homogenization --- src/damage_nonlocal.f90 | 28 ++-------------------------- src/homogenization.f90 | 10 ++++++++++ src/homogenization_damage.f90 | 26 ++++++++++++++++++++++++++ 3 files changed, 38 insertions(+), 26 deletions(-) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 4423c1e3a..1ae327833 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -23,7 +23,7 @@ module damage_nonlocal real(pReal) :: & charLength !< characteristic length scale for gradient problems end type tNumerics - + type(tparameters), dimension(:), allocatable :: & param type(tNumerics), private :: & @@ -33,7 +33,6 @@ module damage_nonlocal damage_nonlocal_init, & damage_nonlocal_getSourceAndItsTangent, & damage_nonlocal_getDiffusion, & - damage_nonlocal_getMobility, & damage_nonlocal_putNonLocalDamage, & damage_nonlocal_results @@ -104,7 +103,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - + call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) @@ -139,29 +138,6 @@ function damage_nonlocal_getDiffusion(ip,el) end function damage_nonlocal_getDiffusion -!-------------------------------------------------------------------------------------------------- -!> @brief Returns homogenized nonlocal damage mobility -!-------------------------------------------------------------------------------------------------- -real(pReal) function damage_nonlocal_getMobility(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - integer :: & - co - - damage_nonlocal_getMobility = 0.0_pReal - - do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) - damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(co,el)) - enddo - - damage_nonlocal_getMobility = damage_nonlocal_getMobility/& - real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end function damage_nonlocal_getMobility - - !-------------------------------------------------------------------------------------------------- !> @brief updated nonlocal damage field with solution from damage phase field PDE !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 048bf46df..816da4fa6 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -145,6 +145,15 @@ module homogenization real(pReal), intent(out) :: Tdot end subroutine thermal_conduction_getSource + module function damage_nonlocal_getMobility(ip,el) result(M) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + integer :: & + co + real(pReal) :: M + end function damage_nonlocal_getMobility + end interface public :: & @@ -154,6 +163,7 @@ module homogenization thermal_conduction_getConductivity, & thermal_conduction_getMassDensity, & thermal_conduction_getSource, & + damage_nonlocal_getMobility, & homogenization_thermal_setfield, & homogenization_thermal_T, & homogenization_forward, & diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index ad029cce8..3a6758914 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -3,6 +3,8 @@ !-------------------------------------------------------------------------------------------------- submodule(homogenization) homogenization_damage + use lattice + type :: tDataContainer real(pReal), dimension(:), allocatable :: phi end type tDataContainer @@ -78,4 +80,28 @@ module subroutine damage_partition(ce) end subroutine damage_partition + +!-------------------------------------------------------------------------------------------------- +!> @brief Returns homogenized nonlocal damage mobility +!-------------------------------------------------------------------------------------------------- +module function damage_nonlocal_getMobility(ip,el) result(M) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + integer :: & + co + real(pReal) :: M + + M = 0.0_pReal + + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + M = M + lattice_M(material_phaseAt(co,el)) + enddo + + M = M/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + +end function damage_nonlocal_getMobility + + end submodule homogenization_damage From 5a35c5ebc38ac6f2bf341ad12ffb3379d324b0a9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 24 Jan 2021 22:44:47 +0100 Subject: [PATCH 03/12] homogenization functionality --- src/damage_nonlocal.f90 | 24 ------------------------ src/homogenization.f90 | 11 +++++++++++ src/homogenization_damage.f90 | 22 ++++++++++++++++++++++ 3 files changed, 33 insertions(+), 24 deletions(-) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 1ae327833..fe11b0db7 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -31,7 +31,6 @@ module damage_nonlocal public :: & damage_nonlocal_init, & - damage_nonlocal_getSourceAndItsTangent, & damage_nonlocal_getDiffusion, & damage_nonlocal_putNonLocalDamage, & damage_nonlocal_results @@ -88,29 +87,6 @@ subroutine damage_nonlocal_init end subroutine damage_nonlocal_init -!-------------------------------------------------------------------------------------------------- -!> @brief calculates homogenized damage driving forces -!-------------------------------------------------------------------------------------------------- -subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - real(pReal) :: & - phiDot, dPhiDot_dPhi - - phiDot = 0.0_pReal - dPhiDot_dPhi = 0.0_pReal - - call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) - phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) - -end subroutine damage_nonlocal_getSourceAndItsTangent - - !-------------------------------------------------------------------------------------------------- !> @brief returns homogenized non local damage diffusion tensor in reference configuration !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 816da4fa6..64169382a 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -154,6 +154,16 @@ module homogenization real(pReal) :: M end function damage_nonlocal_getMobility +module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal) :: & + phiDot, dPhiDot_dPhi +end subroutine damage_nonlocal_getSourceAndItsTangent end interface public :: & @@ -164,6 +174,7 @@ module homogenization thermal_conduction_getMassDensity, & thermal_conduction_getSource, & damage_nonlocal_getMobility, & + damage_nonlocal_getSourceAndItsTangent, & homogenization_thermal_setfield, & homogenization_thermal_T, & homogenization_forward, & diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 3a6758914..630327e8d 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -104,4 +104,26 @@ module function damage_nonlocal_getMobility(ip,el) result(M) end function damage_nonlocal_getMobility +!-------------------------------------------------------------------------------------------------- +!> @brief calculates homogenized damage driving forces +!-------------------------------------------------------------------------------------------------- +module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal) :: & + phiDot, dPhiDot_dPhi + + phiDot = 0.0_pReal + dPhiDot_dPhi = 0.0_pReal + + call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal) + +end subroutine damage_nonlocal_getSourceAndItsTangent + end submodule homogenization_damage From 2533e0616d6fd7dc1777869b61cbe6018f23d073 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jan 2021 07:28:47 +0100 Subject: [PATCH 04/12] don't discriminate by number of constituents/grains --- src/homogenization.f90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 64169382a..17ead5d68 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -274,8 +274,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State(:,me) endif steppingNeeded - elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite - num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep + elseif (subStep <= 1.0 ) then ! cutback makes no sense if (.not. terminallyIll) & ! so first signals terminally ill... print*, ' Integration point ', ip,' at element ', el, ' terminally ill' From 6cca9202ad309c63e9667f6c22d475f57263dd63 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jan 2021 07:42:12 +0100 Subject: [PATCH 05/12] substepping seems to be inactive even in RGC test --- src/homogenization.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 17ead5d68..abc5aa219 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -266,7 +266,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration steppingNeeded: if (subStep > num%subStepMinHomog) then - + error stop ! wind forward grain starting point call constitutive_windForward(ip,el) @@ -294,6 +294,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & .and. NiterationMPstate < num%nMPstate) NiterationMPstate = NiterationMPstate + 1 + if (subStep /= 1.0 .or. subFrac /= 0.0) error stop + !-------------------------------------------------------------------------------------------------- ! deformation partitioning From 2f5c988a8979eeb1aada043112e7cc59b72d9f78 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jan 2021 13:34:06 +0100 Subject: [PATCH 06/12] removed dead code at least in the tests it was not used ... --- src/homogenization.f90 | 23 +++-------------------- 1 file changed, 3 insertions(+), 20 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index abc5aa219..1886e87b6 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -265,15 +265,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE subFrac = subFrac + subStep subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration - steppingNeeded: if (subStep > num%subStepMinHomog) then - error stop - ! wind forward grain starting point - call constitutive_windForward(ip,el) - - if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me) - if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State(:,me) - - endif steppingNeeded elseif (subStep <= 1.0 ) then ! cutback makes no sense if (.not. terminallyIll) & ! so first signals terminally ill... @@ -294,29 +285,21 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & .and. NiterationMPstate < num%nMPstate) NiterationMPstate = NiterationMPstate + 1 - if (subStep /= 1.0 .or. subFrac /= 0.0) error stop - !-------------------------------------------------------------------------------------------------- ! deformation partitioning if (.not. doneAndHappy(1)) then - call mech_partition( homogenization_F0(1:3,1:3,ce) & - + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce))*(subStep+subFrac), & - ip,el) + call mech_partition(homogenization_F(1:3,1:3,ce),ip,el) converged = .true. do co = 1, myNgrains - converged = converged .and. crystallite_stress(dt*subStep,co,ip,el) + converged = converged .and. crystallite_stress(dt,co,ip,el) enddo if (.not. converged) then doneAndHappy = [.true.,.false.] else - doneAndHappy = mech_updateState(dt*subStep, & - homogenization_F0(1:3,1:3,ce) & - + (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) & - *(subStep+subFrac), & - ip,el) + doneAndHappy = mech_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el) converged = all(doneAndHappy) endif endif From 5592f5aa34f96860c27e0fc26132854b8eeae028 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jan 2021 13:49:29 +0100 Subject: [PATCH 07/12] simplified --- src/homogenization.f90 | 75 +++++++++++++++--------------------------- 1 file changed, 27 insertions(+), 48 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 1886e87b6..977bd4078 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -242,7 +242,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE doneAndHappy !$OMP PARALLEL - !$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy) + !$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) myNgrains = homogenization_Nconstituents(ho) @@ -252,60 +252,39 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call constitutive_initializeRestorationPoints(ip,el) - subFrac = 0.0_pReal - converged = .false. ! pretend failed step ... - subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation + call constitutive_restore(ce,.false.) ! wrong name (is more a forward function) - if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me) - if (damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State0(:,me) + if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%State0(:,me) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%State0(:,me) - cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) + doneAndHappy = [.false.,.true.] - if (converged) then - subFrac = subFrac + subStep - subStep = min(1.0_pReal-subFrac,num%stepIncreaseHomog*subStep) ! introduce flexibility for step increase/acceleration + NiterationMPstate = 0 + convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & + .and. NiterationMPstate < num%nMPstate) + NiterationMPstate = NiterationMPstate + 1 - elseif (subStep <= 1.0 ) then - ! cutback makes no sense - if (.not. terminallyIll) & ! so first signals terminally ill... - print*, ' Integration point ', ip,' at element ', el, ' terminally ill' - terminallyIll = .true. ! ...and kills all others - else ! cutback makes sense - subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback - call constitutive_restore(ce,subStep < 1.0_pReal) + if (.not. doneAndHappy(1)) then + call mech_partition(homogenization_F(1:3,1:3,ce),ip,el) + converged = .true. + do co = 1, myNgrains + converged = converged .and. crystallite_stress(dt,co,ip,el) + enddo - if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me) - if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%subState0(:,me) + if (.not. converged) then + doneAndHappy = [.true.,.false.] + else + doneAndHappy = mech_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el) + converged = all(doneAndHappy) + endif endif - if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] - - NiterationMPstate = 0 - convergenceLooping: do while (.not. (terminallyIll .or. doneAndHappy(1)) & - .and. NiterationMPstate < num%nMPstate) - NiterationMPstate = NiterationMPstate + 1 - -!-------------------------------------------------------------------------------------------------- -! deformation partitioning - - if (.not. doneAndHappy(1)) then - call mech_partition(homogenization_F(1:3,1:3,ce),ip,el) - converged = .true. - do co = 1, myNgrains - converged = converged .and. crystallite_stress(dt,co,ip,el) - enddo - - if (.not. converged) then - doneAndHappy = [.true.,.false.] - else - doneAndHappy = mech_updateState(dt,homogenization_F(1:3,1:3,ce),ip,el) - converged = all(doneAndHappy) - endif - endif - - enddo convergenceLooping - enddo cutBackLooping + enddo convergenceLooping + if (.not. converged) then + if (.not. terminallyIll) print*, ' Integration point ', ip,' at element ', el, ' terminally ill' + terminallyIll = .true. + endif enddo enddo !$OMP END DO @@ -337,7 +316,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ho = material_homogenizationAt(el) do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip - call damage_partition(ce) + ! call damage_partition(ce) ! do co = 1, homogenization_Nconstituents(ho) ! ph = material_phaseAt(co,el) ! if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then From 6292094a6591e103b9c69fdee67d8fc738aa6fed Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jan 2021 14:03:49 +0100 Subject: [PATCH 08/12] not needed partitionedState was only for RGC-triggered cutback. subState is for internal iteration, no need to store for all points --- src/constitutive.f90 | 45 ++----------- src/constitutive_mech.f90 | 93 +++++++-------------------- src/constitutive_thermal.f90 | 7 +- src/damage_none.f90 | 1 - src/damage_nonlocal.f90 | 1 - src/homogenization.f90 | 5 -- src/homogenization_mech_RGC.f90 | 1 - src/homogenization_mech_isostrain.f90 | 1 - src/homogenization_mech_none.f90 | 1 - src/prec.f90 | 1 - 10 files changed, 31 insertions(+), 125 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 149d09a49..630756ec4 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -115,12 +115,6 @@ module constitutive integer, intent(in) :: ph end subroutine damage_results - - module subroutine mech_initializeRestorationPoints(ph,me) - integer, intent(in) :: ph, me - end subroutine mech_initializeRestorationPoints - - module subroutine mech_windForward(ph,me) integer, intent(in) :: ph, me end subroutine mech_windForward @@ -359,7 +353,6 @@ module constitutive constitutive_mech_getP, & constitutive_mech_setF, & constitutive_mech_getF, & - constitutive_initializeRestorationPoints, & constitutive_windForward, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & @@ -404,11 +397,9 @@ subroutine constitutive_init PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- ! partition and initialize state - plasticState(ph)%partitionedState0 = plasticState(ph)%state0 - plasticState(ph)%state = plasticState(ph)%partitionedState0 + plasticState(ph)%state = plasticState(ph)%state0 forall(so = 1:phase_Nsources(ph)) - damageState(ph)%p(so)%partitionedState0 = damageState(ph)%p(so)%state0 - damageState(ph)%p(so)%state = damageState(ph)%p(so)%partitionedState0 + damageState(ph)%p(so)%state = damageState(ph)%p(so)%state0 end forall constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & @@ -473,7 +464,6 @@ subroutine constitutive_allocateState(state, & allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal) - allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal) allocate(state%state (sizeState,Nconstituents), source=0.0_pReal) allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal) @@ -500,7 +490,7 @@ subroutine constitutive_restore(ce,includeL) do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do so = 1, phase_Nsources(material_phaseAt2(co,ce)) damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & - damageState(material_phaseAt2(co,ce))%p(so)%partitionedState0(:,material_phasememberAt2(co,ce)) + damageState(material_phaseAt2(co,ce))%p(so)%state0(:,material_phasememberAt2(co,ce)) enddo enddo @@ -651,33 +641,6 @@ subroutine crystallite_init() end subroutine crystallite_init -!-------------------------------------------------------------------------------------------------- -!> @brief Backup data for homog cutback. -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_initializeRestorationPoints(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - integer :: & - co, & !< constituent number - so,ph, me - - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - call mech_initializeRestorationPoints(ph,me) - - do so = 1, size(damageState(ph)%p) - damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me) - enddo - - enddo - -end subroutine constitutive_initializeRestorationPoints - - !-------------------------------------------------------------------------------------------------- !> @brief Wind homog inc forward. !-------------------------------------------------------------------------------------------------- @@ -699,7 +662,7 @@ subroutine constitutive_windForward(ip,el) call mech_windForward(ph,me) do so = 1, phase_Nsources(material_phaseAt(co,el)) - damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state(:,me) + damageState(ph)%p(so)%state0(:,me) = damageState(ph)%p(so)%state(:,me) enddo enddo diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 56e03ac6f..be933ef4f 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -39,16 +39,7 @@ submodule(constitutive) constitutive_mech constitutive_mech_F0, & constitutive_mech_Li0, & constitutive_mech_Lp0, & - constitutive_mech_S0, & - ! converged value at end of last homogenization increment (RGC only) - constitutive_mech_partitionedFi0, & - constitutive_mech_partitionedFp0, & - constitutive_mech_partitionedF0, & - constitutive_mech_partitionedLi0, & - constitutive_mech_partitionedLp0, & - constitutive_mech_partitionedS0 - - + constitutive_mech_S0 integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & @@ -361,23 +352,17 @@ module subroutine mech_init(phases) allocate(constitutive_mech_Fe(phases%length)) allocate(constitutive_mech_Fi(phases%length)) allocate(constitutive_mech_Fi0(phases%length)) - allocate(constitutive_mech_partitionedFi0(phases%length)) allocate(constitutive_mech_Fp(phases%length)) allocate(constitutive_mech_Fp0(phases%length)) - allocate(constitutive_mech_partitionedFp0(phases%length)) allocate(constitutive_mech_F(phases%length)) allocate(constitutive_mech_F0(phases%length)) - allocate(constitutive_mech_partitionedF0(phases%length)) allocate(constitutive_mech_Li(phases%length)) allocate(constitutive_mech_Li0(phases%length)) - allocate(constitutive_mech_partitionedLi0(phases%length)) - allocate(constitutive_mech_partitionedLp0(phases%length)) allocate(constitutive_mech_Lp0(phases%length)) allocate(constitutive_mech_Lp(phases%length)) allocate(constitutive_mech_S(phases%length)) allocate(constitutive_mech_P(phases%length)) allocate(constitutive_mech_S0(phases%length)) - allocate(constitutive_mech_partitionedS0(phases%length)) do ph = 1, phases%length Nconstituents = count(material_phaseAt == ph) * discretization_nIPs @@ -385,23 +370,17 @@ module subroutine mech_init(phases) allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Fe(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedFi0(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedFp0(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedLi0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedLp0(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Lp0(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Lp(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(constitutive_mech_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(constitutive_mech_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) - allocate(constitutive_mech_partitionedS0(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(constitutive_mech_F(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_F0(ph)%data(3,3,Nconstituents)) - allocate(constitutive_mech_partitionedF0(ph)%data(3,3,Nconstituents)) phase => phases%get(ph) mech => phase%get('mechanics') @@ -454,10 +433,6 @@ module subroutine mech_init(phases) constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) constitutive_mech_F(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) - enddo enddo; enddo !$OMP END PARALLEL DO @@ -1464,26 +1439,6 @@ subroutine crystallite_results(group,ph) end subroutine crystallite_results -!-------------------------------------------------------------------------------------------------- -!> @brief Backup data for homog cutback. -!-------------------------------------------------------------------------------------------------- -module subroutine mech_initializeRestorationPoints(ph,me) - - integer, intent(in) :: ph, me - - - constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) - - plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state0(:,me) - -end subroutine mech_initializeRestorationPoints - - !-------------------------------------------------------------------------------------------------- !> @brief Wind homog inc forward. !-------------------------------------------------------------------------------------------------- @@ -1492,14 +1447,14 @@ module subroutine mech_windForward(ph,me) integer, intent(in) :: ph, me - constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) = constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) = constitutive_mech_S(ph)%data(1:3,1:3,me) + constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp(ph)%data(1:3,1:3,me) + constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi(ph)%data(1:3,1:3,me) + constitutive_mech_F0(ph)%data(1:3,1:3,me) = constitutive_mech_F(ph)%data(1:3,1:3,me) + constitutive_mech_Li0(ph)%data(1:3,1:3,me) = constitutive_mech_Li(ph)%data(1:3,1:3,me) + constitutive_mech_Lp0(ph)%data(1:3,1:3,me) = constitutive_mech_Lp(ph)%data(1:3,1:3,me) + constitutive_mech_S0(ph)%data(1:3,1:3,me) = constitutive_mech_S(ph)%data(1:3,1:3,me) - plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state(:,me) + plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me) end subroutine mech_windForward @@ -1578,17 +1533,17 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) me = material_phaseMemberAt(co,ip,el) sizeDotState = plasticState(ph)%sizeDotState - subLi0 = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) - subLp0 = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) - subState0 = plasticState(ph)%partitionedState0(:,me) + subLi0 = constitutive_mech_Li0(ph)%data(1:3,1:3,me) + subLp0 = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) + subState0 = plasticState(ph)%State0(:,me) do so = 1, phase_Nsources(ph) - damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%partitionedState0(:,me) + damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state0(:,me) enddo - subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) - subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) - subF0 = constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me) + subFp0 = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) + subFi0 = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) + subF0 = constitutive_mech_F0(ph)%data(1:3,1:3,me) subFrac = 0.0_pReal subStep = 1.0_pReal/num%subStepSizeCryst todo = .true. @@ -1638,7 +1593,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) ! prepare for integration if (todo) then subF = subF0 & - + subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_partitionedF0(ph)%data(1:3,1:3,me)) + + subStep * (constitutive_mech_F(ph)%data(1:3,1:3,me) - constitutive_mech_F0(ph)%data(1:3,1:3,me)) constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) @@ -1667,15 +1622,15 @@ module subroutine mech_restore(ce,includeL) ph = material_phaseAt2(co,ce) me = material_phaseMemberAt2(co,ce) if (includeL) then - constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLp0(ph)%data(1:3,1:3,me) - constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedLi0(ph)%data(1:3,1:3,me) + constitutive_mech_Lp(ph)%data(1:3,1:3,me) = constitutive_mech_Lp0(ph)%data(1:3,1:3,me) + constitutive_mech_Li(ph)%data(1:3,1:3,me) = constitutive_mech_Li0(ph)%data(1:3,1:3,me) endif ! maybe protecting everything from overwriting makes more sense - constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) - constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) - constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) + constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) + constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me) + constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_S0(ph)%data(1:3,1:3,me) - plasticState(ph)%state(:,me) = plasticState(ph)%partitionedState0(:,me) + plasticState(ph)%state(:,me) = plasticState(ph)%State0(:,me) enddo end subroutine mech_restore @@ -1727,8 +1682,8 @@ module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me)) invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me)) - invSubFp0 = math_inv33(constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)) - invSubFi0 = math_inv33(constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)) + invSubFp0 = math_inv33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me)) + invSubFi0 = math_inv33(constitutive_mech_Fi0(ph)%data(1:3,1:3,me)) if (sum(abs(dLidS)) < tol_math_check) then dFidS = 0.0_pReal diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 57b8a3117..831fd236d 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -116,12 +116,11 @@ module subroutine thermal_init(phases) PhaseLoop2:do ph = 1,phases%length do so = 1,thermal_Nsources(ph) - deallocate(thermalState(ph)%p(so)%partitionedState0) - thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0 + thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%state0 enddo - thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & - maxval(thermalState(ph)%p%sizeDotState)) + thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & + maxval(thermalState(ph)%p%sizeDotState)) enddo PhaseLoop2 !-------------------------------------------------------------------------------------------------- diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 078d42af7..680110d55 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -27,7 +27,6 @@ subroutine damage_none_init Nmaterialpoints = count(material_homogenizationAt == h) damageState_h(h)%sizeState = 0 allocate(damageState_h(h)%state0 (0,Nmaterialpoints)) - allocate(damageState_h(h)%subState0(0,Nmaterialpoints)) allocate(damageState_h(h)%state (0,Nmaterialpoints)) allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index fe11b0db7..6e46d8799 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -76,7 +76,6 @@ subroutine damage_nonlocal_init Nmaterialpoints = count(material_homogenizationAt == h) damageState_h(h)%sizeState = 1 allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState_h(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal) damage(h)%p => damageState_h(h)%state(1,:) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 977bd4078..5523f38ac 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -233,9 +233,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ip, & !< integration point number el, & !< element number myNgrains, co, ce, ho, me, ph - real(pReal) :: & - subFrac, & - subStep logical :: & converged logical, dimension(2) :: & @@ -250,8 +247,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE ce = (el-1)*discretization_nIPs + ip me = material_homogenizationMemberAt2(ce) - call constitutive_initializeRestorationPoints(ip,el) - call constitutive_restore(ce,.false.) ! wrong name (is more a forward function) if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%State0(:,me) diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 24d14e059..2c7a5a8cb 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -171,7 +171,6 @@ module subroutine mech_RGC_init(num_homogMech) homogState(h)%sizeState = sizeState allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal) - allocate(homogState(h)%subState0(sizeState,Nmaterialpoints), source=0.0_pReal) allocate(homogState(h)%state (sizeState,Nmaterialpoints), source=0.0_pReal) stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:) diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index a56104647..df8f5fc9d 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -64,7 +64,6 @@ module subroutine mech_isostrain_init Nmaterialpoints = count(material_homogenizationAt == h) homogState(h)%sizeState = 0 allocate(homogState(h)%state0 (0,Nmaterialpoints)) - allocate(homogState(h)%subState0(0,Nmaterialpoints)) allocate(homogState(h)%state (0,Nmaterialpoints)) end associate diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90 index d434d1ca0..767dbf557 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mech_none.f90 @@ -32,7 +32,6 @@ module subroutine mech_none_init Nmaterialpoints = count(material_homogenizationAt == h) homogState(h)%sizeState = 0 allocate(homogState(h)%state0 (0,Nmaterialpoints)) - allocate(homogState(h)%subState0(0,Nmaterialpoints)) allocate(homogState(h)%state (0,Nmaterialpoints)) enddo diff --git a/src/prec.f90 b/src/prec.f90 index 4d73462c4..eabeecaf4 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -47,7 +47,6 @@ module prec dotState, & !< rate of state change deltaState !< increment of state change real(pReal), allocatable, dimension(:,:) :: & - partitionedState0, & subState0 end type From fdc48a7987757c2b5d9bde33827c8694949eae41 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jan 2021 14:43:09 +0100 Subject: [PATCH 09/12] not used --- src/constitutive.f90 | 3 +-- src/homogenization.f90 | 11 ----------- 2 files changed, 1 insertion(+), 13 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 630756ec4..1412f5259 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -489,7 +489,7 @@ subroutine constitutive_restore(ce,includeL) do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do so = 1, phase_Nsources(material_phaseAt2(co,ce)) - damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & + damageState(material_phaseAt2(co,ce))%p(so)%state( :,material_phasememberAt2(co,ce)) = & damageState(material_phaseAt2(co,ce))%p(so)%state0(:,material_phasememberAt2(co,ce)) enddo enddo @@ -559,7 +559,6 @@ subroutine crystallite_init() iMax, & !< maximum number of integration points eMax !< maximum number of elements - class(tNode), pointer :: & num_crystallite, & debug_crystallite, & ! pointer to debug options for crystallite diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 5523f38ac..9d5ea5d4f 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -38,10 +38,6 @@ module homogenization type :: tNumerics integer :: & nMPstate !< materialpoint state loop limit - real(pReal) :: & - subStepMinHomog, & !< minimum (relative) size of sub-step allowed during cutback in homogenization - subStepSizeHomog, & !< size of first substep when cutback in homogenization - stepIncreaseHomog !< increase of next substep size when previous substep converged in homogenization end type tNumerics type(tNumerics) :: num @@ -200,14 +196,7 @@ subroutine homogenization_init() num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict) num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10) - num%subStepMinHomog = num_homogGeneric%get_asFloat('subStepMin', defaultVal=1.0e-3_pReal) - num%subStepSizeHomog = num_homogGeneric%get_asFloat('subStepSize', defaultVal=0.25_pReal) - num%stepIncreaseHomog = num_homogGeneric%get_asFloat('stepIncrease', defaultVal=1.5_pReal) - if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate') - if (num%subStepMinHomog <= 0.0_pReal) call IO_error(301,ext_msg='subStepMinHomog') - if (num%subStepSizeHomog <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeHomog') - if (num%stepIncreaseHomog <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseHomog') call mech_init(num_homog) From 8b11af0d8477bf64e209dbcfbc5c286f872b175d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jan 2021 15:13:17 +0100 Subject: [PATCH 10/12] more sensible location --- src/damage_nonlocal.f90 | 69 +---------------------------------- src/homogenization.f90 | 37 ++++++++++++++----- src/homogenization_damage.f90 | 44 ++++++++++++++++++++++ 3 files changed, 74 insertions(+), 76 deletions(-) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index fe11b0db7..c4abaccea 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -14,26 +14,17 @@ module damage_nonlocal implicit none private - type :: tParameters - character(len=pStringLen), allocatable, dimension(:) :: & - output - end type tParameters - type, private :: tNumerics real(pReal) :: & charLength !< characteristic length scale for gradient problems end type tNumerics - type(tparameters), dimension(:), allocatable :: & - param type(tNumerics), private :: & num public :: & damage_nonlocal_init, & - damage_nonlocal_getDiffusion, & - damage_nonlocal_putNonLocalDamage, & - damage_nonlocal_results + damage_nonlocal_getDiffusion contains @@ -46,9 +37,7 @@ subroutine damage_nonlocal_init integer :: Ninstances,Nmaterialpoints,h class(tNode), pointer :: & num_generic, & - material_homogenization, & - homog, & - homogDamage + material_homogenization print'(/,a)', ' <<<+- damage_nonlocal init -+>>>'; flush(6) @@ -58,20 +47,10 @@ subroutine damage_nonlocal_init num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) Ninstances = count(damage_type == DAMAGE_nonlocal_ID) - allocate(param(Ninstances)) material_homogenization => config_material%get('homogenization') do h = 1, material_homogenization%length if (damage_type(h) /= DAMAGE_NONLOCAL_ID) cycle - homog => material_homogenization%get(h) - homogDamage => homog%get('damage') - associate(prm => param(damage_typeInstance(h))) - -#if defined (__GFORTRAN__) - prm%output = output_asStrings(homogDamage) -#else - prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray) -#endif Nmaterialpoints = count(material_homogenizationAt == h) damageState_h(h)%sizeState = 1 @@ -81,7 +60,6 @@ subroutine damage_nonlocal_init damage(h)%p => damageState_h(h)%state(1,:) - end associate enddo end subroutine damage_nonlocal_init @@ -114,47 +92,4 @@ function damage_nonlocal_getDiffusion(ip,el) end function damage_nonlocal_getDiffusion -!-------------------------------------------------------------------------------------------------- -!> @brief updated nonlocal damage field with solution from damage phase field PDE -!-------------------------------------------------------------------------------------------------- -subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - integer :: & - homog, & - offset - - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - damage(homog)%p(offset) = phi - -end subroutine damage_nonlocal_putNonLocalDamage - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine damage_nonlocal_results(homog,group) - - integer, intent(in) :: homog - character(len=*), intent(in) :: group - - integer :: o - - associate(prm => param(damage_typeInstance(homog))) - outputsLoop: do o = 1,size(prm%output) - select case(prm%output(o)) - case ('phi') - call results_writeDataset(group,damage(homog)%p,prm%output(o),& - 'damage indicator','-') - end select - enddo outputsLoop - end associate - -end subroutine damage_nonlocal_results - end module damage_nonlocal diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 17ead5d68..9d517dd1f 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -154,16 +154,34 @@ module homogenization real(pReal) :: M end function damage_nonlocal_getMobility -module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) + module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el) - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - real(pReal), intent(in) :: & - phi - real(pReal) :: & - phiDot, dPhiDot_dPhi -end subroutine damage_nonlocal_getSourceAndItsTangent + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal) :: & + phiDot, dPhiDot_dPhi + end subroutine damage_nonlocal_getSourceAndItsTangent + + + module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + + end subroutine damage_nonlocal_putNonLocalDamage + + module subroutine damage_nonlocal_results(homog,group) + + integer, intent(in) :: homog + character(len=*), intent(in) :: group + + end subroutine damage_nonlocal_results end interface public :: & @@ -175,6 +193,7 @@ end subroutine damage_nonlocal_getSourceAndItsTangent thermal_conduction_getSource, & damage_nonlocal_getMobility, & damage_nonlocal_getSourceAndItsTangent, & + damage_nonlocal_putNonLocalDamage, & homogenization_thermal_setfield, & homogenization_thermal_T, & homogenization_forward, & diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 630327e8d..3115e6a6f 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -126,4 +126,48 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p end subroutine damage_nonlocal_getSourceAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief updated nonlocal damage field with solution from damage phase field PDE +!-------------------------------------------------------------------------------------------------- +module subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + integer :: & + homog, & + offset + + homog = material_homogenizationAt(el) + offset = material_homogenizationMemberAt(ip,el) + damage(homog)%p(offset) = phi + +end subroutine damage_nonlocal_putNonLocalDamage + + +!-------------------------------------------------------------------------------------------------- +!> @brief writes results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +module subroutine damage_nonlocal_results(homog,group) + + integer, intent(in) :: homog + character(len=*), intent(in) :: group + + integer :: o + + associate(prm => param(damage_typeInstance(homog))) + outputsLoop: do o = 1,size(prm%output) + select case(prm%output(o)) + case ('phi') + call results_writeDataset(group,damage(homog)%p,prm%output(o),& + 'damage indicator','-') + end select + enddo outputsLoop + end associate + +end subroutine damage_nonlocal_results + end submodule homogenization_damage From 546fcba93d0941751150eedf8b1f93a76c42acdc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 25 Jan 2021 19:53:05 +0100 Subject: [PATCH 11/12] polishing --- src/prec.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/prec.f90 b/src/prec.f90 index eabeecaf4..1a96c75a9 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -32,13 +32,13 @@ module prec real(pReal), dimension(:), pointer :: p end type group_float - ! http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array type :: tState integer :: & sizeState = 0, & !< size of state sizeDotState = 0, & !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates offsetDeltaState = 0, & !< index offset of delta state sizeDeltaState = 0 !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments + ! http://stackoverflow.com/questions/3948210 real(pReal), pointer, dimension(:), contiguous :: & atol real(pReal), pointer, dimension(:,:), contiguous :: & ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous From 7a4ad53b8c12ad9f33cd3843b9a6a88b19766646 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 26 Jan 2021 16:55:19 +0100 Subject: [PATCH 12/12] don't rely on compiler to remove dead loops --- src/homogenization.f90 | 38 +++++++++++++++++++------------------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 48919b46d..9b8c33c2e 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -313,25 +313,25 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE enddo !$OMP END DO - !$OMP DO PRIVATE(ho,ph,ce) - do el = FEsolving_execElem(1),FEsolving_execElem(2) - if (terminallyIll) continue - ho = material_homogenizationAt(el) - do ip = FEsolving_execIP(1),FEsolving_execIP(2) - ce = (el-1)*discretization_nIPs + ip - ! call damage_partition(ce) - ! do co = 1, homogenization_Nconstituents(ho) - ! ph = material_phaseAt(co,el) - ! if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then - ! if (.not. terminallyIll) & ! so first signals terminally ill... - ! print*, ' Integration point ', ip,' at element ', el, ' terminally ill' - ! terminallyIll = .true. ! ...and kills all others - ! endif - ! call thermal_homogenize(ip,el) - ! enddo - enddo - enddo - !$OMP END DO +! !$OMP DO PRIVATE(ho,ph,ce) +! do el = FEsolving_execElem(1),FEsolving_execElem(2) +! if (terminallyIll) continue +! ho = material_homogenizationAt(el) +! do ip = FEsolving_execIP(1),FEsolving_execIP(2) +! ce = (el-1)*discretization_nIPs + ip +! call damage_partition(ce) +! do co = 1, homogenization_Nconstituents(ho) +! ph = material_phaseAt(co,el) +! if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then +! if (.not. terminallyIll) & ! so first signals terminally ill... +! print*, ' Integration point ', ip,' at element ', el, ' terminally ill' +! terminallyIll = .true. ! ...and kills all others +! endif +! call thermal_homogenize(ip,el) +! enddo +! enddo +! enddo +! !$OMP END DO !$OMP DO PRIVATE(ho) elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)