Merge commit 'v3.0.0-alpha2-158-gf4e3c872a'

src/constitutive.f90

@@ -64,10 +64,6 @@ module constitutive
   real(pReal),               dimension(:,:,:,:,:),    allocatable, public :: &
     crystallite_partitionedF                                                                        !< def grad to be reached at end of homog inc
 
-  logical,                    dimension(:,:,:),       allocatable :: &
-    crystallite_converged                                                                           !< convergence flag
-
-
   type :: tTensorContainer
     real(pReal), dimension(:,:,:), allocatable :: data
   end type
@@ -185,10 +181,10 @@ module constitutive
 
 ! == cleaned:end ===================================================================================
 
-    module function crystallite_stress(dt,co,ip,el)
+    module function crystallite_stress(dt,co,ip,el) result(converged_)
       real(pReal), intent(in) :: dt
       integer, intent(in) :: co, ip, el
-      logical :: crystallite_stress
+      logical :: converged_
     end function crystallite_stress
 
     module function constitutive_homogenizedC(co,ip,el) result(C)
@@ -872,10 +868,8 @@ subroutine crystallite_init
            source = crystallite_partitionedF)
 
   allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal)
-
   allocate(crystallite_orientation(cMax,iMax,eMax))
 
-  allocate(crystallite_converged(cMax,iMax,eMax),             source=.true.)
 
   num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
 
@@ -1253,7 +1247,7 @@ end function crystallite_push33ToRef
 !> @brief integrate stress, state with adaptive 1st order explicit Euler method
 !> using Fixed Point Iteration to adapt the stepsize
 !--------------------------------------------------------------------------------------------------
-subroutine integrateSourceState(co,ip,el)
+function integrateSourceState(co,ip,el) result(broken)
 
   integer, intent(in) :: &
     el, &                                                                                            !< element index in element loop
@@ -1273,12 +1267,13 @@ subroutine integrateSourceState(co,ip,el)
     r                                                                                               ! state residuum
   real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
   logical :: &
-    broken
+    broken, converged_
 
 
   ph = material_phaseAt(co,el)
   me = material_phaseMemberAt(co,ip,el)
 
+  converged_ = .true.
   broken = constitutive_thermal_collectDotState(ph,me)
   broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me)
   if(broken) return
@@ -1313,19 +1308,20 @@ subroutine integrateSourceState(co,ip,el)
                       - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * crystallite_subdt(co,ip,el)
       sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) &
                                                 - r(1:size_so(so))
-      crystallite_converged(co,ip,el) = &
+      converged_ = converged_  .and. converged(r(1:size_so(so)), &
-      crystallite_converged(co,ip,el) .and. converged(r(1:size_so(so)), &
+                                               sourceState(ph)%p(so)%state(1:size_so(so),me), &
-                                                   sourceState(ph)%p(so)%state(1:size_so(so),me), &
+                                               sourceState(ph)%p(so)%atol(1:size_so(so)))
-                                                   sourceState(ph)%p(so)%atol(1:size_so(so)))
     enddo
 
-    if(crystallite_converged(co,ip,el)) then
+    if(converged_) then
       broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me)
       exit iteration
     endif
 
   enddo iteration
 
+  broken = broken .or. .not. converged_
+
 
   contains
 
@@ -1349,7 +1345,7 @@ subroutine integrateSourceState(co,ip,el)
 
   end function damper
 
-end subroutine integrateSourceState
+end function integrateSourceState
 
 
 !--------------------------------------------------------------------------------------------------

src/constitutive_mech.f90

@@ -951,7 +951,7 @@ end function integrateStress
 !> @brief integrate stress, state with adaptive 1st order explicit Euler method
 !> using Fixed Point Iteration to adapt the stepsize
 !--------------------------------------------------------------------------------------------------
-subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
+function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F
   real(pReal), intent(in) :: Delta_t
@@ -1004,11 +1004,7 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
                  - plasticState(ph)%dotState (1:size_pl,me) * Delta_t
     plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%state(1:size_pl,me) &
                                          - r(1:size_pl)
-    crystallite_converged(co,ip,el) = converged(r(1:size_pl), &
+    if (converged(r(1:size_pl),plasticState(ph)%state(1:size_pl,me),plasticState(ph)%atol(1:size_pl))) then
-                                             plasticState(ph)%state(1:size_pl,me), &
-                                             plasticState(ph)%atol(1:size_pl))
-
-    if(crystallite_converged(co,ip,el)) then
       broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
                                        constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
       exit iteration
@@ -1016,7 +1012,6 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
 
   enddo iteration
 
-
   contains
 
   !--------------------------------------------------------------------------------------------------
@@ -1039,13 +1034,13 @@ subroutine integrateStateFPI(F_0,F,Delta_t,co,ip,el)
 
   end function damper
 
-end subroutine integrateStateFPI
+end function integrateStateFPI
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief integrate state with 1st order explicit Euler method
 !--------------------------------------------------------------------------------------------------
-subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el)
+function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F
   real(pReal), intent(in) :: Delta_t
@@ -1075,15 +1070,14 @@ subroutine integrateStateEuler(F_0,F,Delta_t,co,ip,el)
   if(broken) return
 
   broken = integrateStress(F,Delta_t,co,ip,el)
-  crystallite_converged(co,ip,el) = .not. broken
 
-end subroutine integrateStateEuler
+end function integrateStateEuler
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief integrate stress, state with 1st order Euler method with adaptive step size
 !--------------------------------------------------------------------------------------------------
-subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el)
+function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F
   real(pReal), intent(in) :: Delta_t
@@ -1123,24 +1117,22 @@ subroutine integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el)
   broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
   if(broken) return
 
+  broken = .not. converged(residuum_plastic(1:sizeDotState) + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, &
+                           plasticState(ph)%state(1:sizeDotState,me), &
+                           plasticState(ph)%atol(1:sizeDotState))
 
-  sizeDotState = plasticState(ph)%sizeDotState
+end function integrateStateAdaptiveEuler
-  crystallite_converged(co,ip,el) = converged(residuum_plastic(1:sizeDotState) &
-                                           + 0.5_pReal * plasticState(ph)%dotState(:,me) * Delta_t, &
-                                           plasticState(ph)%state(1:sizeDotState,me), &
-                                           plasticState(ph)%atol(1:sizeDotState))
-
-end subroutine integrateStateAdaptiveEuler
 
 
 !---------------------------------------------------------------------------------------------------
 !> @brief Integrate state (including stress integration) with the classic Runge Kutta method
 !---------------------------------------------------------------------------------------------------
-subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el)
+function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: co,ip,el
+  logical :: broken
 
   real(pReal), dimension(3,3), parameter :: &
     A = reshape([&
@@ -1153,19 +1145,20 @@ subroutine integrateStateRK4(F_0,F,Delta_t,co,ip,el)
   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]
 
-  call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C)
+  broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C)
 
-end subroutine integrateStateRK4
+end function integrateStateRK4
 
 
 !---------------------------------------------------------------------------------------------------
 !> @brief Integrate state (including stress integration) with the Cash-Carp method
 !---------------------------------------------------------------------------------------------------
-subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el)
+function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: co,ip,el
+  logical :: broken
 
   real(pReal), dimension(5,5), parameter :: &
     A = reshape([&
@@ -1185,16 +1178,16 @@ subroutine integrateStateRKCK45(F_0,F,Delta_t,co,ip,el)
       [2825.0_pReal/27648.0_pReal,    .0_pReal,                18575.0_pReal/48384.0_pReal,&
       13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal,  1._pReal/4._pReal]
 
-  call integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
+  broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
 
-end subroutine integrateStateRKCK45
+end function integrateStateRKCK45
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
 !! embedded explicit Runge-Kutta method
 !--------------------------------------------------------------------------------------------------
-subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
+function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken)
 
   real(pReal), intent(in),dimension(3,3) :: F_0,F
   real(pReal), intent(in) :: Delta_t
@@ -1205,15 +1198,14 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
     el, &                                                                                            !< element index in element loop
     ip, &                                                                                            !< integration point index in ip loop
     co                                                                                               !< grain index in grain loop
+  logical :: broken
 
- integer :: &
+  integer :: &
     stage, &                                                                                        ! stage index in integration stage loop
     n, &
     ph, &
     me, &
     sizeDotState
-  logical :: &
-    broken
   real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
 
 
@@ -1266,10 +1258,8 @@ subroutine integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
   if(broken) return
 
   broken = integrateStress(F,Delta_t,co,ip,el)
-  crystallite_converged(co,ip,el) = .not. broken
 
-
+end function integrateStateRK
-end subroutine integrateStateRK
 
 
 !--------------------------------------------------------------------------------------------------
@@ -1479,15 +1469,14 @@ end function constitutive_homogenizedC
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress (P)
 !--------------------------------------------------------------------------------------------------
-module function crystallite_stress(dt,co,ip,el)
+module function crystallite_stress(dt,co,ip,el) result(converged_)
 
   real(pReal), intent(in) :: dt
   integer, intent(in) :: &
     co, &
     ip, &
     el
-
+  logical :: converged_
-  logical :: crystallite_stress
 
   real(pReal) :: &
     formerSubStep
@@ -1519,7 +1508,7 @@ module function crystallite_stress(dt,co,ip,el)
   subFrac = 0.0_pReal
   subStep = 1.0_pReal/num%subStepSizeCryst
   todo = .true.
-  crystallite_converged(co,ip,el) = .false.                                                      ! pretend failed step of 1/subStepSizeCryst
+  converged_ = .false.                                                      ! pretend failed step of 1/subStepSizeCryst
 
   todo = .true.
   NiterationCrystallite = 0
@@ -1528,7 +1517,7 @@ module function crystallite_stress(dt,co,ip,el)
 
 !--------------------------------------------------------------------------------------------------
 !  wind forward
-    if (crystallite_converged(co,ip,el)) then
+    if (converged_) then
       formerSubStep = subStep
       subFrac = subFrac + subStep
       subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep)
@@ -1579,17 +1568,13 @@ module function crystallite_stress(dt,co,ip,el)
                                              math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
                                                                constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
       crystallite_subdt(co,ip,el) = subStep * dt
-      crystallite_converged(co,ip,el) = .false.
+      converged_ = .not. integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),&
-      call integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),&
+                                        crystallite_subdt(co,ip,el),co,ip,el)
-                          crystallite_subdt(co,ip,el),co,ip,el)
+      converged_ = converged_ .and. .not. integrateSourceState(co,ip,el)
-      call integrateSourceState(co,ip,el)
     endif
 
   enddo cutbackLooping
 
-! return whether converged or not
-  crystallite_stress = crystallite_converged(co,ip,el)
-
 end function crystallite_stress
 
 end submodule constitutive_mech