diff --git a/code/crystallite.f90 b/code/crystallite.f90
index 703d297e3..bb491c7e1 100644
--- a/code/crystallite.f90
+++ b/code/crystallite.f90
@@ -16,107 +16,112 @@
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see .
!
-!##############################################################
-!* $Id$
-!***************************************
-!* Module: CRYSTALLITE *
-!***************************************
-!* contains: *
-!* - _init *
-!* - materialpoint_stressAndItsTangent *
-!* - _partitionDeformation *
-!* - _updateState *
-!* - _stressAndItsTangent *
-!* - _postResults *
-!***************************************
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief crystallite state integration functions and reporting of results
+!--------------------------------------------------------------------------------------------------
module crystallite
-
use prec, only: pReal, pInt
implicit none
- private :: crystallite_integrateStateFPI, &
- crystallite_integrateStateEuler, &
- crystallite_integrateStateAdaptiveEuler, &
- crystallite_integrateStateRK4, &
- crystallite_integrateStateRKCK45, &
- crystallite_integrateStress, &
- crystallite_stateJump
- external :: dgesv
-
-! ****************************************************************
-! *** General variables for the crystallite calculation ***
-! ****************************************************************
- integer(pInt) crystallite_maxSizePostResults
- integer(pInt), dimension(:), allocatable :: crystallite_sizePostResults
- integer(pInt), dimension(:,:), allocatable :: crystallite_sizePostResult
- character(len=64), dimension(:,:), allocatable :: crystallite_output !< name of each post result output
- integer(pInt), dimension (:,:,:), allocatable :: &
- crystallite_symmetryID !< crystallographic symmetry 1=cubic 2=hexagonal, needed in all orientation calcs
+ private
+ character(len=64), dimension(:,:), allocatable, private :: &
+ crystallite_output !< name of each post result output
+ integer(pInt), public, protected :: &
+ crystallite_maxSizePostResults
+ integer(pInt), dimension(:), allocatable, public, protected :: &
+ crystallite_sizePostResults
+ integer(pInt), dimension(:,:), allocatable, private :: &
+ crystallite_sizePostResult
+ integer(pInt), dimension(:,:,:), allocatable, private :: &
+ crystallite_symmetryID !< crystallographic symmetry 1=cubic 2=hexagonal, needed in all orientation calcs
- real(pReal), dimension (:,:,:), allocatable :: &
- crystallite_dt, & !< requested time increment of each grain
- crystallite_subdt, & !< substepped time increment of each grain
- crystallite_subFrac, & !< already calculated fraction of increment
- crystallite_subStep, & !< size of next integration step
- crystallite_Temperature, & !< Temp of each grain
- crystallite_partionedTemperature0, & !< Temp of each grain at start of homog inc
- crystallite_subTemperature0, & !< Temp of each grain at start of crystallite inc
- crystallite_dotTemperature !< evolution of Temperature of each grain
- real(pReal), dimension (:,:,:,:), allocatable :: &
- crystallite_Tstar_v, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
- crystallite_Tstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
- crystallite_partionedTstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of homog inc
- crystallite_subTstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of crystallite inc
- crystallite_orientation, & !< orientation as quaternion
- crystallite_orientation0, & !< initial orientation as quaternion
- crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame
- real(pReal), dimension (:,:,:,:,:), allocatable :: &
- crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
- crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc
- crystallite_Fp, & !< current plastic def grad (end of converged time step)
- crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step)
- crystallite_Fp0, & !< plastic def grad at start of FE inc
- crystallite_partionedFp0,& !< plastic def grad at start of homog inc
- crystallite_subFp0,& !< plastic def grad at start of crystallite inc
- crystallite_F0, & !< def grad at start of FE inc
- crystallite_partionedF, & !< def grad to be reached at end of homog inc
- crystallite_partionedF0, & !< def grad at start of homog inc
- crystallite_subF, & !< def grad to be reached at end of crystallite inc
- crystallite_subF0, & !< def grad at start of crystallite inc
- crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step)
- crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
- crystallite_partionedLp0,& !< plastic velocity grad at start of homog inc
- crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc
- crystallite_P, & !< 1st Piola-Kirchhoff stress per grain
- crystallite_disorientation !< disorientation between two neighboring ips (only calculated for single grain IPs)
- real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: &
- crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step)
- crystallite_dPdF0, & !< individual dPdF per grain at start of FE inc
- crystallite_partioneddPdF0, & !< individual dPdF per grain at start of homog inc
- crystallite_fallbackdPdF !< dPdF fallback for non-converged grains (elastic prediction)
- logical, dimension (:,:,:), allocatable :: &
- crystallite_localPlasticity, & !< indicates this grain to have purely local constitutive law
- crystallite_requested, & !< flag to request crystallite calculation
- crystallite_todo, & !< flag to indicate need for further computation
- crystallite_converged !< convergence flag
- logical, dimension (:,:), allocatable :: &
- crystallite_clearToWindForward, &
- crystallite_clearToCutback, &
- crystallite_syncSubFrac, &
- crystallite_syncSubFracCompleted, &
- crystallite_neighborEnforcedCutback
-
+ real(pReal), dimension (:,:,:), allocatable, public :: &
+ crystallite_dt, & !< requested time increment of each grain
+ crystallite_Temperature, & !< Temp of each grain
+ crystallite_partionedTemperature0 !< Temp of each grain at start of homog inc
+ real(pReal), dimension (:,:,:), allocatable, private :: &
+ crystallite_subdt, & !< substepped time increment of each grain
+ crystallite_subFrac, & !< already calculated fraction of increment
+ crystallite_subStep, & !< size of next integration step
+ crystallite_subTemperature0, & !< Temp of each grain at start of crystallite inc
+ crystallite_dotTemperature !< evolution of Temperature of each grain
+ real(pReal), dimension (:,:,:,:), allocatable, public :: &
+ crystallite_Tstar_v, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
+ crystallite_Tstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
+ crystallite_partionedTstar0_v !< 2nd Piola-Kirchhoff stress vector at start of homog inc
+ real(pReal), dimension (:,:,:,:), allocatable, private :: &
+ crystallite_subTstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of crystallite inc
+ crystallite_orientation, & !< orientation as quaternion
+ crystallite_orientation0, & !< initial orientation as quaternion
+ crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame
+ real(pReal), dimension (:,:,:,:,:), allocatable, public :: &
+ crystallite_Fp, & !< current plastic def grad (end of converged time step)
+ crystallite_Fp0, & !< plastic def grad at start of FE inc
+ crystallite_partionedFp0,& !< plastic def grad at start of homog inc
+ crystallite_F0, & !< def grad at start of FE inc
+ crystallite_partionedF, & !< def grad to be reached at end of homog inc
+ crystallite_partionedF0, & !< def grad at start of homog inc
+ crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step)
+ crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
+ crystallite_partionedLp0,& !< plastic velocity grad at start of homog inc
+ crystallite_P !< 1st Piola-Kirchhoff stress per grain
+ real(pReal), dimension (:,:,:,:,:), allocatable, private :: &
+ crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
+ crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc
+ crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step)
+ crystallite_subFp0,& !< plastic def grad at start of crystallite inc
+ crystallite_subF, & !< def grad to be reached at end of crystallite inc
+ crystallite_subF0, & !< def grad at start of crystallite inc
+ crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc
+ crystallite_disorientation !< disorientation between two neighboring ips (only calculated for single grain IPs)
+ real(pReal), dimension (:,:,:,:,:,:,:), allocatable, public :: &
+ crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step)
+ crystallite_dPdF0, & !< individual dPdF per grain at start of FE inc
+ crystallite_partioneddPdF0 !< individual dPdF per grain at start of homog inc
+ real(pReal), dimension (:,:,:,:,:,:,:), allocatable, private :: &
+ crystallite_fallbackdPdF !< dPdF fallback for non-converged grains (elastic prediction)
+ logical, dimension (:,:,:), allocatable, public :: &
+ crystallite_requested !< flag to request crystallite calculation
+ logical, dimension (:,:,:), allocatable, public, protected :: &
+ crystallite_converged !< convergence flag
+ logical, dimension (:,:,:), allocatable, private :: &
+ crystallite_localPlasticity, & !< indicates this grain to have purely local constitutive law
+ crystallite_todo !< flag to indicate need for further computation
+ logical, dimension (:,:), allocatable, private :: &
+ crystallite_clearToWindForward, &
+ crystallite_clearToCutback, &
+ crystallite_syncSubFrac, &
+ crystallite_syncSubFracCompleted, &
+ crystallite_neighborEnforcedCutback
+
+ public :: &
+ crystallite_init, &
+ crystallite_stressAndItsTangent, &
+ crystallite_orientations, &
+ crystallite_postResults
+ private :: &
+ crystallite_integrateStateFPI, &
+ crystallite_integrateStateEuler, &
+ crystallite_integrateStateAdaptiveEuler, &
+ crystallite_integrateStateRK4, &
+ crystallite_integrateStateRKCK45, &
+ crystallite_integrateStress, &
+ crystallite_stateJump
+
contains
-!********************************************************************
-! allocate and initialize per grain variables
-!********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief allocates and initialize per grain variables
+!--------------------------------------------------------------------------------------------------
subroutine crystallite_init(Temperature)
-
- !*** variables and functions from other modules ***!
- use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use debug, only: debug_info, &
debug_reset, &
debug_level, &
@@ -149,11 +154,9 @@ subroutine crystallite_init(Temperature)
constitutive_nonlocal_structureName
implicit none
+ real(pReal), intent(in) :: Temperature
integer(pInt), parameter :: myFile = 200_pInt, &
maxNchunks = 2_pInt
-
- !*** input variables ***!
- real(pReal) Temperature
!*** local variables ***!
integer(pInt), dimension(1+2*maxNchunks) :: positions
@@ -175,13 +178,9 @@ subroutine crystallite_init(Temperature)
character(len=64) tag
character(len=1024) line
-
-!$OMP CRITICAL (write2out)
- write(6,*)
- write(6,*) '<<<+- crystallite init -+>>>'
- write(6,*) '$Id$'
+ write(6,'(/,a)') ' <<<+- crystallite init -+>>>'
+ write(6,'(a)') ' $Id$'
#include "compilation_info.f90"
-!$OMP END CRITICAL (write2out)
gMax = homogenization_maxNgrains
@@ -243,27 +242,27 @@ subroutine crystallite_init(Temperature)
material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt
- if (.not. IO_open_jobFile_stat(myFile,material_localFileExt)) then ! no local material configuration present...
- call IO_open_file(myFile,material_configFile) ! ...open material.config file
+ if (.not. IO_open_jobFile_stat(myFile,material_localFileExt)) then ! no local material configuration present...
+ call IO_open_file(myFile,material_configFile) ! ...open material.config file
endif
line = ''
section = 0_pInt
- do while (IO_lc(IO_getTag(line,'<','>')) /= material_partCrystallite) ! wind forward to
+ do while (IO_lc(IO_getTag(line,'<','>')) /= material_partCrystallite) ! wind forward to
read(myFile,'(a1024)',END=100) line
enddo
- do ! read thru sections of phase part
+ do ! read thru sections of phase part
read(myFile,'(a1024)',END=100) line
- if (IO_isBlank(line)) cycle ! skip empty lines
- if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
- if (IO_getTag(line,'[',']') /= '') then ! next section
+ if (IO_isBlank(line)) cycle ! skip empty lines
+ if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
+ if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt
- output = 0_pInt ! reset output counter
+ output = 0_pInt ! reset output counter
endif
if (section > 0_pInt) then
positions = IO_stringPos(line,maxNchunks)
- tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
+ tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('(output)')
output = output + 1_pInt
@@ -274,7 +273,7 @@ subroutine crystallite_init(Temperature)
100 close(myFile)
- do i = 1_pInt,material_Ncrystallite ! sanity checks
+ do i = 1_pInt,material_Ncrystallite ! sanity checks
enddo
do i = 1_pInt,material_Ncrystallite
@@ -282,9 +281,9 @@ subroutine crystallite_init(Temperature)
select case(crystallite_output(j,i))
case('phase','texture','volume','grainrotationx','grainrotationy','grainrotationz')
mySize = 1_pInt
- case('orientation','grainrotation') ! orientation as quaternion, or deviation from initial grain orientation in axis-angle form (angle in degrees)
+ case('orientation','grainrotation') ! orientation as quaternion, or deviation from initial grain orientation in axis-angle form (angle in degrees)
mySize = 4_pInt
- case('eulerangles') ! Bunge (3-1-3) Euler angles
+ case('eulerangles') ! Bunge (3-1-3) Euler angles
mySize = 3_pInt
case('defgrad','f','fe','fp','lp','e','ee','p','firstpiola','1stpiola','s','tstar','secondpiola','2ndpiola')
mySize = 9_pInt
@@ -294,7 +293,7 @@ subroutine crystallite_init(Temperature)
mySize = 0_pInt
end select
- if (mySize > 0_pInt) then ! any meaningful output found
+ if (mySize > 0_pInt) then ! any meaningful output found
crystallite_sizePostResult(j,i) = mySize
crystallite_sizePostResults(i) = crystallite_sizePostResults(i) + mySize
endif
@@ -384,9 +383,8 @@ crystallite_orientation0 = crystallite_orientation ! Store initial o
call crystallite_stressAndItsTangent(.true.,.false.) ! request elastic answers
crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback
-! *** Output to MARC output file ***
+! *** Output ***
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- !$OMP CRITICAL (write2out)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Temperature: ', shape(crystallite_Temperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_dotTemperature: ', shape(crystallite_dotTemperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe)
@@ -433,7 +431,6 @@ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,*)
write(6,*) 'Number of nonlocal grains: ',count(.not. crystallite_localPlasticity)
flush(6)
- !$OMP END CRITICAL (write2out)
endif
call debug_info
@@ -442,13 +439,10 @@ endif
end subroutine crystallite_init
-
-!********************************************************************
-! calculate stress (P) and tangent (dPdF) for crystallites
-!********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief calculate stress (P) and tangent (dPdF) for crystallites
+!--------------------------------------------------------------------------------------------------
subroutine crystallite_stressAndItsTangent(updateJaco,rate_sensitivity)
-
-!*** variables and functions from other modules ***!
use numerics, only: subStepMinCryst, &
subStepSizeCryst, &
stepIncreaseCryst, &
@@ -502,9 +496,7 @@ use constitutive, only: constitutive_sizeState, &
implicit none
-!*** input variables ***!
logical, intent(in) :: updateJaco, rate_sensitivity ! flag indicating wehther we want to update the Jacobian (stiffness) or not
-!*** local variables ***!
real(pReal) myPert, & ! perturbation with correct sign
formerSubStep, &
subFracIntermediate
@@ -1026,12 +1018,12 @@ if(updateJaco) then
! --- CALCULATE STATE AND STRESS FOR PERTURBATION ---
- dPdF_perturbation1 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
- dPdF_perturbation2 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
- do perturbation = 1,2 ! forward and backward perturbation
- if (iand(pert_method,perturbation) > 0_pInt) then ! mask for desired direction
- myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step
- do k = 1,3; do l = 1,3 ! ...alter individual components
+ dPdF_perturbation1 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
+ dPdF_perturbation2 = crystallite_dPdF0 ! initialize stiffness with known good values from last increment
+ do perturbation = 1,2 ! forward and backward perturbation
+ if (iand(pert_method,perturbation) > 0_pInt) then ! mask for desired direction
+ myPert = -pert_Fg * (-1.0_pReal)**perturbation ! set perturbation step
+ do k = 1,3; do l = 1,3 ! ...alter individual components
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(a,2(1x,i1),1x,a)') '<< CRYST >> [[[[[[ Stiffness perturbation',k,l,']]]]]]'
@@ -1085,7 +1077,7 @@ if(updateJaco) then
crystallite_converged = convergenceFlag_backup
crystallite_todo = crystallite_requested .and. crystallite_converged
- where (crystallite_todo) crystallite_converged = .false. ! start out non-converged
+ where (crystallite_todo) crystallite_converged = .false. ! start out non-converged
select case(numerics_integrator(numerics_integrationMode))
case(1_pInt)
@@ -1105,16 +1097,16 @@ if(updateJaco) then
select case(perturbation)
case(1_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
- crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) & ! converged state warrants stiffness update
+ crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) & ! converged state warrants stiffness update
dPdF_perturbation1(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl
case(2_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
- crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) & ! converged state warrants stiffness update
+ crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) & ! converged state warrants stiffness update
dPdF_perturbation2(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl
end select
enddo
- enddo; enddo ! k,l component perturbation loop
+ enddo; enddo ! k,l component perturbation loop
endif
enddo ! perturbation direction
@@ -1127,21 +1119,21 @@ if(updateJaco) then
select case(pert_method)
case(1_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
- crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) & ! perturbation mode 1: central solution converged
+ crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) & ! perturbation mode 1: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e)
case(2_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
- crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) & ! perturbation mode 2: central solution converged
+ crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) & ! perturbation mode 2: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e)
case(3_pInt)
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
- crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) & ! perturbation mode 3: central solution converged
+ crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) & ! perturbation mode 3: central solution converged
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) &
+ dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e))
end select
forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains, &
- crystallite_requested(g,i,e) .and. .not. convergenceFlag_backup(g,i,e)) & ! for any pertubation mode: if central solution did not converge...
- crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! ...use (elastic) fallback
+ crystallite_requested(g,i,e) .and. .not. convergenceFlag_backup(g,i,e)) & ! for any pertubation mode: if central solution did not converge...
+ crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = crystallite_fallbackdPdF(1:3,1:3,1:3,1:3,g,i,e) ! ...use (elastic) fallback
enddo elementLooping
@@ -1166,25 +1158,25 @@ if(updateJaco) then
crystallite_P = P_backup
crystallite_converged = convergenceFlag_backup
- else ! Calculate Jacobian using analytical expression
+ else ! Calculate Jacobian using analytical expression
! --- CALCULATE ANALYTIC dPdF ---
!$OMP PARALLEL DO PRIVATE(dFedF,dSdF,dSdFe,myNgrains)
- do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
+ do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
- do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
+ do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1_pInt,myNgrains
dFedF = 0.0_pReal
do p=1_pInt,3_pInt; do o=1_pInt,3_pInt
- dFedF(p,o,o,1:3) = crystallite_invFp(1:3,p,g,i,e) ! dFe^T_ij/dF_kl = delta_jk * (Fp current^-1)_li
+ dFedF(p,o,o,1:3) = crystallite_invFp(1:3,p,g,i,e) ! dFe^T_ij/dF_kl = delta_jk * (Fp current^-1)_li
enddo; enddo
- call constitutive_TandItsTangent(junk,dSdFe,crystallite_subFe0(1:3,1:3,g,i,e),g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative
- dSdF = math_mul3333xx3333(dSdFe,dFedF) ! dS/dF = dS/dFe * dFe/dF
+ call constitutive_TandItsTangent(junk,dSdFe,crystallite_subFe0(1:3,1:3,g,i,e),g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative
+ dSdF = math_mul3333xx3333(dSdFe,dFedF) ! dS/dF = dS/dFe * dFe/dF
do p=1_pInt,3_pInt; do o=1_pInt,3_pInt
crystallite_dPdF(1:3,1:3,o,p,g,i,e) = math_mul33x33(math_mul33x33(dFedF(1:3,1:3,o,p),&
math_Mandel6to33(crystallite_Tstar_v)),math_transpose33(&
- crystallite_invFp(1:3,1:3,g,i,e))) & ! dP/dF = dFe/dF * S * Fp^-T...
+ crystallite_invFp(1:3,1:3,g,i,e))) & ! dP/dF = dFe/dF * S * Fp^-T...
+ math_mul33x33(crystallite_subFe0(1:3,1:3,g,i,e),&
math_mul33x33(dSdF(1:3,1:3,o,p),math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))) ! + Fe * dS/dF * Fp^-T
enddo; enddo
@@ -1237,2255 +1229,2206 @@ endif
end subroutine crystallite_stressAndItsTangent
-
-!********************************************************************
-! integrate stress, state and Temperature with
-! 4h order explicit Runge Kutta method
-!********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief integrate stress, state and Temperature with 4th order explicit Runge Kutta method
+!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRK4(gg,ii,ee)
-
-!*** variables and functions from other modules ***!
-use prec, only: pInt, &
- pReal
-use numerics, only: numerics_integrationMode
-use debug, only: debug_level, &
- debug_crystallite, &
- debug_levelBasic, &
- debug_levelExtensive, &
- debug_levelSelective, &
- debug_e, &
- debug_i, &
- debug_g, &
- debug_StateLoopDistribution
-use FEsolving, only: FEsolving_execElem, &
- FEsolving_execIP
-use mesh, only: mesh_element, &
- mesh_NcpElems, &
- mesh_maxNips
-use material, only: homogenization_Ngrains, &
- homogenization_maxNgrains
-use constitutive, only: constitutive_sizeDotState, &
- constitutive_state, &
- constitutive_subState0, &
- constitutive_dotState, &
- constitutive_RK4dotState, &
- constitutive_collectDotState, &
- constitutive_deltaState, &
- constitutive_collectDeltaState, &
- constitutive_dotTemperature, &
- constitutive_microstructure
-
-implicit none
-
-real(pReal), dimension(4), parameter :: timeStepFraction = [0.5_pReal, 0.5_pReal, 1.0_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
-real(pReal), dimension(4), parameter :: weight = [1.0_pReal, 2.0_pReal, 2.0_pReal, 1.0_pReal] ! weight of slope used for Runge Kutta integration
-
-!*** input variables ***!
-integer(pInt), optional, intent(in):: ee, & ! element index
- ii, & ! integration point index
- gg ! grain index
-
-!*** output variables ***!
-
-!*** local variables ***!
-integer(pInt) e, & ! element index in element loop
- i, & ! integration point index in ip loop
- g, & ! grain index in grain loop
- n, &
- mySizeDotState
-integer(pInt), dimension(2) :: eIter ! bounds for element iteration
-integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
- gIter ! bounds for grain iteration
-real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
- RK4dotTemperature ! evolution of Temperature of each grain for Runge Kutta integration
-logical singleRun ! flag indicating computation for single (g,i,e) triple
-
-
-if (present(ee) .and. present(ii) .and. present(gg)) then
- eIter = ee
- iIter(1:2,ee) = ii
- gIter(1:2,ee) = gg
- singleRun = .true.
-else
- eIter = FEsolving_execElem(1:2)
- do e = eIter(1),eIter(2)
- iIter(1:2,e) = FEsolving_execIP(1:2,e)
- gIter(1:2,e) = [ 1_pInt,homogenization_Ngrains(mesh_element(3,e))]
- enddo
- singleRun = .false.
-endif
-
-
-! --- FIRST RUNGE KUTTA STEP ---
-
-RK4dotTemperature = 0.0_pReal ! initialize Runge-Kutta dotTemperature
-!$OMP PARALLEL PRIVATE(mySizeDotState)
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- constitutive_RK4dotState(g,i,e)%p = 0.0_pReal ! initialize Runge-Kutta dotState
- if (crystallite_todo(g,i,e)) then
- call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
- crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
- crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
- crystallite_Temperature(g,i,e),g,i,e)
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
- .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
- if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- else ! if broken local...
- crystallite_todo(g,i,e) = .false. ! ... skip this one next time
- endif
- endif
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-!$OMP END PARALLEL
-
-
-! --- SECOND TO FOURTH RUNGE KUTTA STEP PLUS FINAL INTEGRATION ---
-
-do n = 1_pInt,4_pInt
-
- ! --- state update ---
-
- !$OMP PARALLEL PRIVATE(mySizeDotState)
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- if (n < 4) then
- constitutive_RK4dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p + weight(n)*constitutive_dotState(g,i,e)%p
- RK4dotTemperature(g,i,e) = RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e)
- elseif (n == 4) then
- constitutive_dotState(g,i,e)%p = (constitutive_RK4dotState(g,i,e)%p + &
- weight(n)*constitutive_dotState(g,i,e)%p) / 6.0_pReal ! use weighted RKdotState for final integration
- crystallite_dotTemperature(g,i,e) = (RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e)) / 6.0_pReal
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
- + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) * timeStepFraction(n)
- crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
- + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) * timeStepFraction(n)
- if (n == 4) then ! final integration step
+ use prec, only: pInt, &
+ pReal
+ use numerics, only: numerics_integrationMode
+ use debug, only: debug_level, &
+ debug_crystallite, &
+ debug_levelBasic, &
+ debug_levelExtensive, &
+ debug_levelSelective, &
+ debug_e, &
+ debug_i, &
+ debug_g, &
+ debug_StateLoopDistribution
+ use FEsolving, only: FEsolving_execElem, &
+ FEsolving_execIP
+ use mesh, only: mesh_element, &
+ mesh_NcpElems, &
+ mesh_maxNips
+ use material, only: homogenization_Ngrains, &
+ homogenization_maxNgrains
+ use constitutive, only: constitutive_sizeDotState, &
+ constitutive_state, &
+ constitutive_subState0, &
+ constitutive_dotState, &
+ constitutive_RK4dotState, &
+ constitutive_collectDotState, &
+ constitutive_deltaState, &
+ constitutive_collectDeltaState, &
+ constitutive_dotTemperature, &
+ constitutive_microstructure
+
+ implicit none
+
+ real(pReal), dimension(4), parameter :: timeStepFraction = [0.5_pReal, 0.5_pReal, 1.0_pReal, 1.0_pReal] ! factor giving the fraction of the original timestep used for Runge Kutta Integration
+ real(pReal), dimension(4), parameter :: weight = [1.0_pReal, 2.0_pReal, 2.0_pReal, 1.0_pReal] ! weight of slope used for Runge Kutta integration
+
+ !*** input variables ***!
+ integer(pInt), optional, intent(in):: ee, & ! element index
+ ii, & ! integration point index
+ gg ! grain index
+
+ !*** output variables ***!
+
+ !*** local variables ***!
+ integer(pInt) e, & ! element index in element loop
+ i, & ! integration point index in ip loop
+ g, & ! grain index in grain loop
+ n, &
+ mySizeDotState
+ integer(pInt), dimension(2) :: eIter ! bounds for element iteration
+ integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
+ gIter ! bounds for grain iteration
+ real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+ RK4dotTemperature ! evolution of Temperature of each grain for Runge Kutta integration
+ logical singleRun ! flag indicating computation for single (g,i,e) triple
+
+
+ if (present(ee) .and. present(ii) .and. present(gg)) then
+ eIter = ee
+ iIter(1:2,ee) = ii
+ gIter(1:2,ee) = gg
+ singleRun = .true.
+ else
+ eIter = FEsolving_execElem(1:2)
+ do e = eIter(1),eIter(2)
+ iIter(1:2,e) = FEsolving_execIP(1:2,e)
+ gIter(1:2,e) = [ 1_pInt,homogenization_Ngrains(mesh_element(3,e))]
+ enddo
+ singleRun = .false.
+ endif
+
+
+ ! --- FIRST RUNGE KUTTA STEP ---
+
+ RK4dotTemperature = 0.0_pReal ! initialize Runge-Kutta dotTemperature
+ !$OMP PARALLEL PRIVATE(mySizeDotState)
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ constitutive_RK4dotState(g,i,e)%p = 0.0_pReal ! initialize Runge-Kutta dotState
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
+ crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
+ crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
+ crystallite_Temperature(g,i,e),g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
+ .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
+ if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ else ! if broken local...
+ crystallite_todo(g,i,e) = .false. ! ... skip this one next time
+ endif
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP END PARALLEL
+
+
+ ! --- SECOND TO FOURTH RUNGE KUTTA STEP PLUS FINAL INTEGRATION ---
+
+ do n = 1_pInt,4_pInt
+
+ ! --- state update ---
+
+ !$OMP PARALLEL PRIVATE(mySizeDotState)
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ if (n < 4) then
+ constitutive_RK4dotState(g,i,e)%p = constitutive_RK4dotState(g,i,e)%p + weight(n)*constitutive_dotState(g,i,e)%p
+ RK4dotTemperature(g,i,e) = RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e)
+ elseif (n == 4) then
+ constitutive_dotState(g,i,e)%p = (constitutive_RK4dotState(g,i,e)%p + &
+ weight(n)*constitutive_dotState(g,i,e)%p) / 6.0_pReal ! use weighted RKdotState for final integration
+ crystallite_dotTemperature(g,i,e) = (RK4dotTemperature(g,i,e) + weight(n)*crystallite_dotTemperature(g,i,e)) / 6.0_pReal
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e) * timeStepFraction(n)
+ crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
+ + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e) * timeStepFraction(n)
+ if (n == 4) then ! final integration step
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
- write(6,*)
- write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ endif
#endif
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- state jump ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- update dependent states ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
- crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- stress integration ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e,timeStepFraction(n)) ! fraction of original times step
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- dot state and RK dot state---
-
- if (n < 4) then
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
- crystallite_Temperature(g,i,e), timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
- crystallite_subFrac, g,i,e)
- crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
- crystallite_Temperature(g,i,e),g,i,e)
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
- .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
- if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- else ! if broken local...
- crystallite_todo(g,i,e) = .false. ! ... skip this one next time
- endif
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- endif
- !$OMP END PARALLEL
-
-enddo
-
-
-! --- SET CONVERGENCE FLAG ---
-
-do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- crystallite_converged(g,i,e) = .true. ! if still "to do" then converged per definitionem
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- !$OMP CRITICAL (distributionState)
- debug_StateLoopDistribution(4,numerics_integrationMode) = &
- debug_StateLoopDistribution(4,numerics_integrationMode) + 1_pInt
- !$OMP END CRITICAL (distributionState)
- endif
- endif
-enddo; enddo; enddo
-
-
-! --- CHECK NONLOCAL CONVERGENCE ---
-
-if (.not. singleRun) then ! if not requesting Integration of just a single IP
- if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) then ! any non-local not yet converged (or broken)...
- crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
- endif
-endif
-
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- state jump ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- update dependent states ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
+ crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- stress integration ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e,timeStepFraction(n)) ! fraction of original times step
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- dot state and RK dot state---
+
+ if (n < 4) then
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
+ crystallite_Temperature(g,i,e), timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
+ crystallite_subFrac, g,i,e)
+ crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
+ crystallite_Temperature(g,i,e),g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
+ .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
+ if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ else ! if broken local...
+ crystallite_todo(g,i,e) = .false. ! ... skip this one next time
+ endif
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ endif
+ !$OMP END PARALLEL
+
+ enddo
+
+
+ ! --- SET CONVERGENCE FLAG ---
+
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ crystallite_converged(g,i,e) = .true. ! if still "to do" then converged per definitionem
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ !$OMP CRITICAL (distributionState)
+ debug_StateLoopDistribution(4,numerics_integrationMode) = &
+ debug_StateLoopDistribution(4,numerics_integrationMode) + 1_pInt
+ !$OMP END CRITICAL (distributionState)
+ endif
+ endif
+ enddo; enddo; enddo
+
+
+ ! --- CHECK NONLOCAL CONVERGENCE ---
+
+ if (.not. singleRun) then ! if not requesting Integration of just a single IP
+ if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) then ! any non-local not yet converged (or broken)...
+ crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
+ endif
+ endif
+
end subroutine crystallite_integrateStateRK4
-
-!********************************************************************
-! integrate stress, state and Temperature with
-! 5th order Runge-Kutta Cash-Karp method with adaptive step size
-! (use 5th order solution to advance = "local extrapolation")
-!********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief integrate stress, state and Temperature with 5th order Runge-Kutta Cash-Karp method with
+!> adaptive step size (use 5th order solution to advance = "local extrapolation")
+!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateRKCK45(gg,ii,ee)
-
-!*** variables and functions from other modules ***!
-use debug, only: debug_level, &
- debug_crystallite, &
- debug_levelBasic, &
- debug_levelExtensive, &
- debug_levelSelective, &
- debug_e, &
- debug_i, &
- debug_g, &
- debug_StateLoopDistribution
-use numerics, only: rTol_crystalliteState, &
- rTol_crystalliteTemperature, &
- numerics_integrationMode
-use FEsolving, only: FEsolving_execElem, &
- FEsolving_execIP
-use mesh, only: mesh_element, &
- mesh_NcpElems, &
- mesh_maxNips
-use material, only: homogenization_Ngrains, &
- homogenization_maxNgrains
-use constitutive, only: constitutive_sizeDotState, &
- constitutive_maxSizeDotState, &
- constitutive_state, &
- constitutive_aTolState, &
- constitutive_subState0, &
- constitutive_dotState, &
- constitutive_RKCK45dotState, &
- constitutive_collectDotState, &
- constitutive_deltaState, &
- constitutive_collectDeltaState, &
- constitutive_dotTemperature, &
- constitutive_microstructure
-
-implicit none
-!*** input variables ***!
-integer(pInt), optional, intent(in):: ee, & ! element index
- ii, & ! integration point index
- gg ! grain index
-!*** local variables ***!
-integer(pInt) e, & ! element index in element loop
- i, & ! integration point index in ip loop
- g, & ! grain index in grain loop
- n, & ! stage index in integration stage loop
- mySizeDotState, & ! size of dot State
- s ! state index
-integer(pInt), dimension(2) :: eIter ! bounds for element iteration
-integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
- gIter ! bounds for grain iteration
-real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
- RKCK45dotTemperature ! evolution of Temperature of each grain for Runge Kutta Cash Karp integration
-real(pReal), dimension(5,5) :: a ! coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
-real(pReal), dimension(6) :: b, db ! coefficients in Butcher tableau (used for final integration and error estimate)
-real(pReal), dimension(5) :: c ! coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
-real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
- stateResiduum, & ! residuum from evolution in micrstructure
- relStateResiduum ! relative residuum from evolution in microstructure
-real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
- temperatureResiduum, & ! residuum from evolution in temperature
- relTemperatureResiduum ! relative residuum from evolution in temperature
-logical singleRun ! flag indicating computation for single (g,i,e) triple
-
-
-! --- FILL BUTCHER TABLEAU ---
-
-a = 0.0_pReal
-b = 0.0_pReal
-db = 0.0_pReal
-c = 0.0_pReal
-
-a(1,1) = 0.2_pReal
-a(1,2) = 0.075_pReal
-a(2,2) = 0.225_pReal
-a(1,3) = 0.3_pReal
-a(2,3) = -0.9_pReal
-a(3,3) = 1.2_pReal
-a(1,4) = -11.0_pReal / 54.0_pReal
-a(2,4) = 2.5_pReal
-a(3,4) = -70.0_pReal / 27.0_pReal
-a(4,4) = 35.0_pReal / 27.0_pReal
-a(1,5) = 1631.0_pReal / 55296.0_pReal
-a(2,5) = 175.0_pReal / 512.0_pReal
-a(3,5) = 575.0_pReal / 13824.0_pReal
-a(4,5) = 44275.0_pReal / 110592.0_pReal
-a(5,5) = 253.0_pReal / 4096.0_pReal
-
-b(1) = 37.0_pReal / 378.0_pReal
-b(3) = 250.0_pReal / 621.0_pReal
-b(4) = 125.0_pReal / 594.0_pReal
-b(6) = 512.0_pReal / 1771.0_pReal
-
-db(1) = b(1) - 2825.0_pReal / 27648.0_pReal
-db(3) = b(3) - 18575.0_pReal / 48384.0_pReal
-db(4) = b(4) - 13525.0_pReal / 55296.0_pReal
-db(5) = - 277.0_pReal / 14336.0_pReal
-db(6) = b(6) - 0.25_pReal
-
-c(1) = 0.2_pReal
-c(2) = 0.3_pReal
-c(3) = 0.6_pReal
-c(4) = 1.0_pReal
-c(5) = 0.875_pReal
-
-
-! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
-
-if (present(ee) .and. present(ii) .and. present(gg)) then
- eIter = ee
- iIter(1:2,ee) = ii
- gIter(1:2,ee) = gg
- singleRun = .true.
-else
- eIter = FEsolving_execElem(1:2)
- do e = eIter(1),eIter(2)
- iIter(1:2,e) = FEsolving_execIP(1:2,e)
- gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
- enddo
- singleRun = .false.
-endif
-
-
-
-! --- FIRST RUNGE KUTTA STEP ---
-
-if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
- !$OMP CRITICAL (write2out)
- write(6,'(a,1x,i1)') '<< CRYST >> RUNGE KUTTA STEP',1
- !$OMP END CRITICAL (write2out)
-endif
-!$OMP PARALLEL PRIVATE(mySizeDotState)
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
- crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
- crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
- crystallite_Temperature(g,i,e),g,i,e)
- endif
+ use debug, only: debug_level, &
+ debug_crystallite, &
+ debug_levelBasic, &
+ debug_levelExtensive, &
+ debug_levelSelective, &
+ debug_e, &
+ debug_i, &
+ debug_g, &
+ debug_StateLoopDistribution
+ use numerics, only: rTol_crystalliteState, &
+ rTol_crystalliteTemperature, &
+ numerics_integrationMode
+ use FEsolving, only: FEsolving_execElem, &
+ FEsolving_execIP
+ use mesh, only: mesh_element, &
+ mesh_NcpElems, &
+ mesh_maxNips
+ use material, only: homogenization_Ngrains, &
+ homogenization_maxNgrains
+ use constitutive, only: constitutive_sizeDotState, &
+ constitutive_maxSizeDotState, &
+ constitutive_state, &
+ constitutive_aTolState, &
+ constitutive_subState0, &
+ constitutive_dotState, &
+ constitutive_RKCK45dotState, &
+ constitutive_collectDotState, &
+ constitutive_deltaState, &
+ constitutive_collectDeltaState, &
+ constitutive_dotTemperature, &
+ constitutive_microstructure
+
+ implicit none
+ !*** input variables ***!
+ integer(pInt), optional, intent(in):: ee, & ! element index
+ ii, & ! integration point index
+ gg ! grain index
+ !*** local variables ***!
+ integer(pInt) e, & ! element index in element loop
+ i, & ! integration point index in ip loop
+ g, & ! grain index in grain loop
+ n, & ! stage index in integration stage loop
+ mySizeDotState, & ! size of dot State
+ s ! state index
+ integer(pInt), dimension(2) :: eIter ! bounds for element iteration
+ integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
+ gIter ! bounds for grain iteration
+ real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+ RKCK45dotTemperature ! evolution of Temperature of each grain for Runge Kutta Cash Karp integration
+ real(pReal), dimension(5,5) :: a ! coefficients in Butcher tableau (used for preliminary integration in stages 2 to 6)
+ real(pReal), dimension(6) :: b, db ! coefficients in Butcher tableau (used for final integration and error estimate)
+ real(pReal), dimension(5) :: c ! coefficients in Butcher tableau (fractions of original time step in stages 2 to 6)
+ real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+ stateResiduum, & ! residuum from evolution in micrstructure
+ relStateResiduum ! relative residuum from evolution in microstructure
+ real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+ temperatureResiduum, & ! residuum from evolution in temperature
+ relTemperatureResiduum ! relative residuum from evolution in temperature
+ logical singleRun ! flag indicating computation for single (g,i,e) triple
+
+
+ ! --- FILL BUTCHER TABLEAU ---
+
+ a = 0.0_pReal
+ b = 0.0_pReal
+ db = 0.0_pReal
+ c = 0.0_pReal
+
+ a(1,1) = 0.2_pReal
+ a(1,2) = 0.075_pReal
+ a(2,2) = 0.225_pReal
+ a(1,3) = 0.3_pReal
+ a(2,3) = -0.9_pReal
+ a(3,3) = 1.2_pReal
+ a(1,4) = -11.0_pReal / 54.0_pReal
+ a(2,4) = 2.5_pReal
+ a(3,4) = -70.0_pReal / 27.0_pReal
+ a(4,4) = 35.0_pReal / 27.0_pReal
+ a(1,5) = 1631.0_pReal / 55296.0_pReal
+ a(2,5) = 175.0_pReal / 512.0_pReal
+ a(3,5) = 575.0_pReal / 13824.0_pReal
+ a(4,5) = 44275.0_pReal / 110592.0_pReal
+ a(5,5) = 253.0_pReal / 4096.0_pReal
+
+ b(1) = 37.0_pReal / 378.0_pReal
+ b(3) = 250.0_pReal / 621.0_pReal
+ b(4) = 125.0_pReal / 594.0_pReal
+ b(6) = 512.0_pReal / 1771.0_pReal
+
+ db(1) = b(1) - 2825.0_pReal / 27648.0_pReal
+ db(3) = b(3) - 18575.0_pReal / 48384.0_pReal
+ db(4) = b(4) - 13525.0_pReal / 55296.0_pReal
+ db(5) = - 277.0_pReal / 14336.0_pReal
+ db(6) = b(6) - 0.25_pReal
+
+ c(1) = 0.2_pReal
+ c(2) = 0.3_pReal
+ c(3) = 0.6_pReal
+ c(4) = 1.0_pReal
+ c(5) = 0.875_pReal
+
+
+ ! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
+
+ if (present(ee) .and. present(ii) .and. present(gg)) then
+ eIter = ee
+ iIter(1:2,ee) = ii
+ gIter(1:2,ee) = gg
+ singleRun = .true.
+ else
+ eIter = FEsolving_execElem(1:2)
+ do e = eIter(1),eIter(2)
+ iIter(1:2,e) = FEsolving_execIP(1:2,e)
+ gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
+ enddo
+ singleRun = .false.
+ endif
+
+
+
+ ! --- FIRST RUNGE KUTTA STEP ---
+
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
+ !$OMP CRITICAL (write2out)
+ write(6,'(a,1x,i1)') '<< CRYST >> RUNGE KUTTA STEP',1
+ !$OMP END CRITICAL (write2out)
+ endif
+ !$OMP PARALLEL PRIVATE(mySizeDotState)
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
+ crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
+ crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
+ crystallite_Temperature(g,i,e),g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
+ .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
+ if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ else ! if broken local...
+ crystallite_todo(g,i,e) = .false. ! ... skip this one next time
+ endif
+ endif
+ endif
enddo; enddo; enddo
-!$OMP ENDDO
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
- .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
- if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- else ! if broken local...
- crystallite_todo(g,i,e) = .false. ! ... skip this one next time
- endif
- endif
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-!$OMP END PARALLEL
-
-
-! --- SECOND TO SIXTH RUNGE KUTTA STEP ---
-
-do n = 1_pInt,5_pInt
-
- ! --- state update ---
-
- !$OMP PARALLEL PRIVATE(mySizeDotState)
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- constitutive_RKCK45dotState(n,g,i,e)%p = constitutive_dotState(g,i,e)%p ! store Runge-Kutta dotState
- RKCK45dotTemperature(n,g,i,e) = crystallite_dotTemperature(g,i,e) ! store Runge-Kutta dotTemperature
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- if (n == 1) then ! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION (CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE)
- constitutive_dotState(g,i,e)%p = a(1,1) * constitutive_RKCK45dotState(1,g,i,e)%p
- crystallite_dotTemperature(g,i,e) = a(1,1) * RKCK45dotTemperature(1,g,i,e)
- elseif (n == 2) then
- constitutive_dotState(g,i,e)%p = a(1,2) * constitutive_RKCK45dotState(1,g,i,e)%p &
- + a(2,2) * constitutive_RKCK45dotState(2,g,i,e)%p
- crystallite_dotTemperature(g,i,e) = a(1,2) * RKCK45dotTemperature(1,g,i,e) &
- + a(2,2) * RKCK45dotTemperature(2,g,i,e)
- elseif (n == 3) then
- constitutive_dotState(g,i,e)%p = a(1,3) * constitutive_RKCK45dotState(1,g,i,e)%p &
- + a(2,3) * constitutive_RKCK45dotState(2,g,i,e)%p &
- + a(3,3) * constitutive_RKCK45dotState(3,g,i,e)%p
- crystallite_dotTemperature(g,i,e) = a(1,3) * RKCK45dotTemperature(1,g,i,e) &
- + a(2,3) * RKCK45dotTemperature(2,g,i,e) &
- + a(3,3) * RKCK45dotTemperature(3,g,i,e)
- elseif (n == 4) then
- constitutive_dotState(g,i,e)%p = a(1,4) * constitutive_RKCK45dotState(1,g,i,e)%p &
- + a(2,4) * constitutive_RKCK45dotState(2,g,i,e)%p &
- + a(3,4) * constitutive_RKCK45dotState(3,g,i,e)%p &
- + a(4,4) * constitutive_RKCK45dotState(4,g,i,e)%p
- crystallite_dotTemperature(g,i,e) = a(1,4) * RKCK45dotTemperature(1,g,i,e) &
- + a(2,4) * RKCK45dotTemperature(2,g,i,e) &
- + a(3,4) * RKCK45dotTemperature(3,g,i,e) &
- + a(4,4) * RKCK45dotTemperature(4,g,i,e)
- elseif (n == 5) then
- constitutive_dotState(g,i,e)%p = a(1,5) * constitutive_RKCK45dotState(1,g,i,e)%p &
- + a(2,5) * constitutive_RKCK45dotState(2,g,i,e)%p &
- + a(3,5) * constitutive_RKCK45dotState(3,g,i,e)%p &
- + a(4,5) * constitutive_RKCK45dotState(4,g,i,e)%p &
- + a(5,5) * constitutive_RKCK45dotState(5,g,i,e)%p
- crystallite_dotTemperature(g,i,e) = a(1,5) * RKCK45dotTemperature(1,g,i,e) &
- + a(2,5) * RKCK45dotTemperature(2,g,i,e) &
- + a(3,5) * RKCK45dotTemperature(3,g,i,e) &
- + a(4,5) * RKCK45dotTemperature(4,g,i,e) &
- + a(5,5) * RKCK45dotTemperature(5,g,i,e)
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
- + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
- crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
- + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- state jump ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- update dependent states ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
- crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- stress integration ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e,c(n)) ! fraction of original time step
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- dot state and RK dot state---
+ !$OMP ENDDO
+ !$OMP END PARALLEL
+
+
+ ! --- SECOND TO SIXTH RUNGE KUTTA STEP ---
+
+ do n = 1_pInt,5_pInt
+
+ ! --- state update ---
+
+ !$OMP PARALLEL PRIVATE(mySizeDotState)
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ constitutive_RKCK45dotState(n,g,i,e)%p = constitutive_dotState(g,i,e)%p ! store Runge-Kutta dotState
+ RKCK45dotTemperature(n,g,i,e) = crystallite_dotTemperature(g,i,e) ! store Runge-Kutta dotTemperature
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ if (n == 1) then ! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION (CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE)
+ constitutive_dotState(g,i,e)%p = a(1,1) * constitutive_RKCK45dotState(1,g,i,e)%p
+ crystallite_dotTemperature(g,i,e) = a(1,1) * RKCK45dotTemperature(1,g,i,e)
+ elseif (n == 2) then
+ constitutive_dotState(g,i,e)%p = a(1,2) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ + a(2,2) * constitutive_RKCK45dotState(2,g,i,e)%p
+ crystallite_dotTemperature(g,i,e) = a(1,2) * RKCK45dotTemperature(1,g,i,e) &
+ + a(2,2) * RKCK45dotTemperature(2,g,i,e)
+ elseif (n == 3) then
+ constitutive_dotState(g,i,e)%p = a(1,3) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ + a(2,3) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ + a(3,3) * constitutive_RKCK45dotState(3,g,i,e)%p
+ crystallite_dotTemperature(g,i,e) = a(1,3) * RKCK45dotTemperature(1,g,i,e) &
+ + a(2,3) * RKCK45dotTemperature(2,g,i,e) &
+ + a(3,3) * RKCK45dotTemperature(3,g,i,e)
+ elseif (n == 4) then
+ constitutive_dotState(g,i,e)%p = a(1,4) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ + a(2,4) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ + a(3,4) * constitutive_RKCK45dotState(3,g,i,e)%p &
+ + a(4,4) * constitutive_RKCK45dotState(4,g,i,e)%p
+ crystallite_dotTemperature(g,i,e) = a(1,4) * RKCK45dotTemperature(1,g,i,e) &
+ + a(2,4) * RKCK45dotTemperature(2,g,i,e) &
+ + a(3,4) * RKCK45dotTemperature(3,g,i,e) &
+ + a(4,4) * RKCK45dotTemperature(4,g,i,e)
+ elseif (n == 5) then
+ constitutive_dotState(g,i,e)%p = a(1,5) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ + a(2,5) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ + a(3,5) * constitutive_RKCK45dotState(3,g,i,e)%p &
+ + a(4,5) * constitutive_RKCK45dotState(4,g,i,e)%p &
+ + a(5,5) * constitutive_RKCK45dotState(5,g,i,e)%p
+ crystallite_dotTemperature(g,i,e) = a(1,5) * RKCK45dotTemperature(1,g,i,e) &
+ + a(2,5) * RKCK45dotTemperature(2,g,i,e) &
+ + a(3,5) * RKCK45dotTemperature(3,g,i,e) &
+ + a(4,5) * RKCK45dotTemperature(4,g,i,e) &
+ + a(5,5) * RKCK45dotTemperature(5,g,i,e)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
+ crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
+ + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- state jump ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- update dependent states ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
+ crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- stress integration ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e,c(n)) ! fraction of original time step
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- dot state and RK dot state---
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
- write(6,'(a,1x,i1)') '<< CRYST >> Runge--Kutta step',n+1_pInt
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
+ write(6,'(a,1x,i1)') '<< CRYST >> Runge--Kutta step',n+1_pInt
+ endif
#endif
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
- crystallite_Temperature(g,i,e), c(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
- crystallite_subFrac, g,i,e)
- crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
- crystallite_Temperature(g,i,e),g,i,e)
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
- .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
- if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- else ! if broken local...
- crystallite_todo(g,i,e) = .false. ! ... skip this one next time
- endif
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP END PARALLEL
-
-enddo
-
-
-! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE AND TEMPERATURE ---
-
-relStateResiduum = 0.0_pReal
-relTemperatureResiduum = 0.0_pReal
-!$OMP PARALLEL PRIVATE(mySizeDotState)
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- constitutive_RKCK45dotState(6,g,i,e)%p = constitutive_dotState(g,i,e)%p ! store Runge-Kutta dotState
- RKCK45dotTemperature(6,g,i,e) = crystallite_dotTemperature(g,i,e) ! store Runge-Kutta dotTemperature
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
-
- ! --- absolute residuum in state and temperature ---
- ! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION
- ! CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE
-
- stateResiduum(1:mySizeDotState,g,i,e) = &
- ( db(1) * constitutive_RKCK45dotState(1,g,i,e)%p(1:mySizeDotState) &
- + db(2) * constitutive_RKCK45dotState(2,g,i,e)%p(1:mySizeDotState) &
- + db(3) * constitutive_RKCK45dotState(3,g,i,e)%p(1:mySizeDotState) &
- + db(4) * constitutive_RKCK45dotState(4,g,i,e)%p(1:mySizeDotState) &
- + db(5) * constitutive_RKCK45dotState(5,g,i,e)%p(1:mySizeDotState) &
- + db(6) * constitutive_RKCK45dotState(6,g,i,e)%p(1:mySizeDotState)) &
- * crystallite_subdt(g,i,e)
- temperatureResiduum(g,i,e) = ( db(1) * RKCK45dotTemperature(1,g,i,e) &
- + db(2) * RKCK45dotTemperature(2,g,i,e) &
- + db(3) * RKCK45dotTemperature(3,g,i,e) &
- + db(4) * RKCK45dotTemperature(4,g,i,e) &
- + db(5) * RKCK45dotTemperature(5,g,i,e) &
- + db(6) * RKCK45dotTemperature(6,g,i,e)) &
- * crystallite_subdt(g,i,e)
-
- ! --- dot state and dot temperature ---
-
- constitutive_dotState(g,i,e)%p = b(1) * constitutive_RKCK45dotState(1,g,i,e)%p &
- + b(2) * constitutive_RKCK45dotState(2,g,i,e)%p &
- + b(3) * constitutive_RKCK45dotState(3,g,i,e)%p &
- + b(4) * constitutive_RKCK45dotState(4,g,i,e)%p &
- + b(5) * constitutive_RKCK45dotState(5,g,i,e)%p &
- + b(6) * constitutive_RKCK45dotState(6,g,i,e)%p
- crystallite_dotTemperature(g,i,e) = b(1) * RKCK45dotTemperature(1,g,i,e) &
- + b(2) * RKCK45dotTemperature(2,g,i,e) &
- + b(3) * RKCK45dotTemperature(3,g,i,e) &
- + b(4) * RKCK45dotTemperature(4,g,i,e) &
- + b(5) * RKCK45dotTemperature(5,g,i,e) &
- + b(6) * RKCK45dotTemperature(6,g,i,e)
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-! --- state and temperature update ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
- + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
- crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
- + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-! --- relative residui and state convergence ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) &
- relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s)
- if (crystallite_Temperature(g,i,e) > 0) then
- relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e)
- endif
- !$OMP FLUSH(relStateResiduum,relTemperatureResiduum)
-
- crystallite_todo(g,i,e) = &
- ( all( abs(relStateResiduum(:,g,i,e)) < rTol_crystalliteState &
- .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) ) &
- .and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature )
-
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
+ crystallite_Temperature(g,i,e), c(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
+ crystallite_subFrac, g,i,e)
+ crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
+ crystallite_Temperature(g,i,e),g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
+ .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
+ if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ else ! if broken local...
+ crystallite_todo(g,i,e) = .false. ! ... skip this one next time
+ endif
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP END PARALLEL
+
+ enddo
+
+
+ ! --- STATE UPDATE WITH ERROR ESTIMATE FOR STATE AND TEMPERATURE ---
+
+ relStateResiduum = 0.0_pReal
+ relTemperatureResiduum = 0.0_pReal
+ !$OMP PARALLEL PRIVATE(mySizeDotState)
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ constitutive_RKCK45dotState(6,g,i,e)%p = constitutive_dotState(g,i,e)%p ! store Runge-Kutta dotState
+ RKCK45dotTemperature(6,g,i,e) = crystallite_dotTemperature(g,i,e) ! store Runge-Kutta dotTemperature
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+
+ ! --- absolute residuum in state and temperature ---
+ ! NEED TO DO THE ADDITION IN THIS LENGTHY WAY BECAUSE OF PARALLELIZATION
+ ! CAN'T USE A REDUCTION CLAUSE ON A POINTER OR USER DEFINED TYPE
+
+ stateResiduum(1:mySizeDotState,g,i,e) = &
+ ( db(1) * constitutive_RKCK45dotState(1,g,i,e)%p(1:mySizeDotState) &
+ + db(2) * constitutive_RKCK45dotState(2,g,i,e)%p(1:mySizeDotState) &
+ + db(3) * constitutive_RKCK45dotState(3,g,i,e)%p(1:mySizeDotState) &
+ + db(4) * constitutive_RKCK45dotState(4,g,i,e)%p(1:mySizeDotState) &
+ + db(5) * constitutive_RKCK45dotState(5,g,i,e)%p(1:mySizeDotState) &
+ + db(6) * constitutive_RKCK45dotState(6,g,i,e)%p(1:mySizeDotState)) &
+ * crystallite_subdt(g,i,e)
+ temperatureResiduum(g,i,e) = ( db(1) * RKCK45dotTemperature(1,g,i,e) &
+ + db(2) * RKCK45dotTemperature(2,g,i,e) &
+ + db(3) * RKCK45dotTemperature(3,g,i,e) &
+ + db(4) * RKCK45dotTemperature(4,g,i,e) &
+ + db(5) * RKCK45dotTemperature(5,g,i,e) &
+ + db(6) * RKCK45dotTemperature(6,g,i,e)) &
+ * crystallite_subdt(g,i,e)
+
+ ! --- dot state and dot temperature ---
+
+ constitutive_dotState(g,i,e)%p = b(1) * constitutive_RKCK45dotState(1,g,i,e)%p &
+ + b(2) * constitutive_RKCK45dotState(2,g,i,e)%p &
+ + b(3) * constitutive_RKCK45dotState(3,g,i,e)%p &
+ + b(4) * constitutive_RKCK45dotState(4,g,i,e)%p &
+ + b(5) * constitutive_RKCK45dotState(5,g,i,e)%p &
+ + b(6) * constitutive_RKCK45dotState(6,g,i,e)%p
+ crystallite_dotTemperature(g,i,e) = b(1) * RKCK45dotTemperature(1,g,i,e) &
+ + b(2) * RKCK45dotTemperature(2,g,i,e) &
+ + b(3) * RKCK45dotTemperature(3,g,i,e) &
+ + b(4) * RKCK45dotTemperature(4,g,i,e) &
+ + b(5) * RKCK45dotTemperature(5,g,i,e) &
+ + b(6) * RKCK45dotTemperature(6,g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ ! --- state and temperature update ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
+ crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
+ + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ ! --- relative residui and state convergence ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) &
+ relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s)
+ if (crystallite_Temperature(g,i,e) > 0) then
+ relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e)
+ endif
+ !$OMP FLUSH(relStateResiduum,relTemperatureResiduum)
+
+ crystallite_todo(g,i,e) = &
+ ( all( abs(relStateResiduum(:,g,i,e)) < rTol_crystalliteState &
+ .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) ) &
+ .and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature )
+
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt&
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,'(a,i8,1x,i3,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
- write(6,*)
- write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> absolute residuum tolerance', &
- stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> relative residuum tolerance', &
- relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
- write(6,*)
- write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt&
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,'(a,i8,1x,i3,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
+ write(6,*)
+ write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> absolute residuum tolerance', &
+ stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> relative residuum tolerance', &
+ relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ endif
#endif
- endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- STATE JUMP ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- UPDATE DEPENDENT STATES IF RESIDUUM BELOW TOLERANCE ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
+ crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
+ endif
enddo; enddo; enddo
-!$OMP ENDDO
-
-
-! --- STATE JUMP ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-
-! --- UPDATE DEPENDENT STATES IF RESIDUUM BELOW TOLERANCE ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
- crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-
-! --- FINAL STRESS INTEGRATION STEP IF RESIDUUM BELOW TOLERANCE ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-
-! --- SET CONVERGENCE FLAG ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- crystallite_converged(g,i,e) = .true. ! if still "to do" then converged per definitionem
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- !$OMP CRITICAL (distributionState)
- debug_StateLoopDistribution(6,numerics_integrationMode) = &
- debug_StateLoopDistribution(6,numerics_integrationMode) + 1_pInt
- !$OMP END CRITICAL (distributionState)
- endif
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-!$OMP END PARALLEL
-
-
-! --- nonlocal convergence check ---
-
-if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
- !$OMP CRITICAL (write2out)
- write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged'
- write(6,*)
- !$OMP END CRITICAL (write2out)
-endif
-if (.not. singleRun) then ! if not requesting Integration of just a single IP
- if ( any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) then ! any non-local not yet converged (or broken)...
- crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
- endif
-
-endif
-
+ !$OMP ENDDO
+
+
+ ! --- FINAL STRESS INTEGRATION STEP IF RESIDUUM BELOW TOLERANCE ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- SET CONVERGENCE FLAG ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ crystallite_converged(g,i,e) = .true. ! if still "to do" then converged per definitionem
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ !$OMP CRITICAL (distributionState)
+ debug_StateLoopDistribution(6,numerics_integrationMode) = &
+ debug_StateLoopDistribution(6,numerics_integrationMode) + 1_pInt
+ !$OMP END CRITICAL (distributionState)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ !$OMP END PARALLEL
+
+
+ ! --- nonlocal convergence check ---
+
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
+ !$OMP CRITICAL (write2out)
+ write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged'
+ write(6,*)
+ !$OMP END CRITICAL (write2out)
+ endif
+ if (.not. singleRun) then ! if not requesting Integration of just a single IP
+ if ( any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) then ! any non-local not yet converged (or broken)...
+ crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
+ endif
+
+ endif
+
end subroutine crystallite_integrateStateRKCK45
-
-!********************************************************************
-! integrate stress, state and Temperature with
-! 1nd order Euler method with adaptive step size
-!********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief integrate stress, state and Temperature with 1st order Euler method with adaptive step size
+!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateAdaptiveEuler(gg,ii,ee)
-
-!*** variables and functions from other modules ***!
-use debug, only: debug_level, &
- debug_crystallite, &
- debug_levelBasic, &
- debug_levelExtensive, &
- debug_levelSelective, &
- debug_e, &
- debug_i, &
- debug_g, &
- debug_StateLoopDistribution
-use numerics, only: rTol_crystalliteState, &
- rTol_crystalliteTemperature, &
- numerics_integrationMode
-use FEsolving, only: FEsolving_execElem, &
- FEsolving_execIP
-use mesh, only: mesh_element, &
- mesh_NcpElems, &
- mesh_maxNips
-use material, only: homogenization_Ngrains, &
- homogenization_maxNgrains
-use constitutive, only: constitutive_sizeDotState, &
- constitutive_maxSizeDotState, &
- constitutive_state, &
- constitutive_aTolState, &
- constitutive_subState0, &
- constitutive_dotState, &
- constitutive_collectDotState, &
- constitutive_dotTemperature, &
- constitutive_microstructure
-
-implicit none
-
-
-!*** input variables ***!
-integer(pInt), optional, intent(in):: ee, & ! element index
- ii, & ! integration point index
- gg ! grain index
-!*** local variables ***!
-integer(pInt) e, & ! element index in element loop
- i, & ! integration point index in ip loop
- g, & ! grain index in grain loop
- mySizeDotState, & ! size of dot State
- s ! state index
-integer(pInt), dimension(2) :: eIter ! bounds for element iteration
-integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
- gIter ! bounds for grain iteration
-real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
- stateResiduum, & ! residuum from evolution in micrstructure
- relStateResiduum ! relative residuum from evolution in microstructure
-real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
- temperatureResiduum, & ! residuum from evolution in temperature
- relTemperatureResiduum ! relative residuum from evolution in temperature
-logical singleRun ! flag indicating computation for single (g,i,e) triple
-
-
-! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
-
-if (present(ee) .and. present(ii) .and. present(gg)) then
- eIter = ee
- iIter(1:2,ee) = ii
- gIter(1:2,ee) = gg
- singleRun = .true.
-else
- eIter = FEsolving_execElem(1:2)
- do e = eIter(1),eIter(2)
- iIter(1:2,e) = FEsolving_execIP(1:2,e)
- gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
- enddo
- singleRun = .false.
-endif
-
-
-stateResiduum = 0.0_pReal
-
-!$OMP PARALLEL PRIVATE(mySizeDotState)
-
-if (numerics_integrationMode == 1_pInt) then
-
- ! --- DOT STATE AND TEMPERATURE (EULER INTEGRATION) ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
- crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
- crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
- crystallite_Temperature(g,i,e),g,i,e)
- endif
+ use debug, only: debug_level, &
+ debug_crystallite, &
+ debug_levelBasic, &
+ debug_levelExtensive, &
+ debug_levelSelective, &
+ debug_e, &
+ debug_i, &
+ debug_g, &
+ debug_StateLoopDistribution
+ use numerics, only: rTol_crystalliteState, &
+ rTol_crystalliteTemperature, &
+ numerics_integrationMode
+ use FEsolving, only: FEsolving_execElem, &
+ FEsolving_execIP
+ use mesh, only: mesh_element, &
+ mesh_NcpElems, &
+ mesh_maxNips
+ use material, only: homogenization_Ngrains, &
+ homogenization_maxNgrains
+ use constitutive, only: constitutive_sizeDotState, &
+ constitutive_maxSizeDotState, &
+ constitutive_state, &
+ constitutive_aTolState, &
+ constitutive_subState0, &
+ constitutive_dotState, &
+ constitutive_collectDotState, &
+ constitutive_dotTemperature, &
+ constitutive_microstructure
+
+ implicit none
+
+ !*** input variables ***!
+ integer(pInt), optional, intent(in):: ee, & ! element index
+ ii, & ! integration point index
+ gg ! grain index
+ !*** local variables ***!
+ integer(pInt) e, & ! element index in element loop
+ i, & ! integration point index in ip loop
+ g, & ! grain index in grain loop
+ mySizeDotState, & ! size of dot State
+ s ! state index
+ integer(pInt), dimension(2) :: eIter ! bounds for element iteration
+ integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
+ gIter ! bounds for grain iteration
+ real(pReal), dimension(constitutive_maxSizeDotState,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+ stateResiduum, & ! residuum from evolution in micrstructure
+ relStateResiduum ! relative residuum from evolution in microstructure
+ real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
+ temperatureResiduum, & ! residuum from evolution in temperature
+ relTemperatureResiduum ! relative residuum from evolution in temperature
+ logical singleRun ! flag indicating computation for single (g,i,e) triple
+
+
+ ! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP ---
+
+ if (present(ee) .and. present(ii) .and. present(gg)) then
+ eIter = ee
+ iIter(1:2,ee) = ii
+ gIter(1:2,ee) = gg
+ singleRun = .true.
+ else
+ eIter = FEsolving_execElem(1:2)
+ do e = eIter(1),eIter(2)
+ iIter(1:2,e) = FEsolving_execIP(1:2,e)
+ gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
+ enddo
+ singleRun = .false.
+ endif
+
+
+ stateResiduum = 0.0_pReal
+
+ !$OMP PARALLEL PRIVATE(mySizeDotState)
+
+ if (numerics_integrationMode == 1_pInt) then
+
+ ! --- DOT STATE AND TEMPERATURE (EULER INTEGRATION) ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
+ crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
+ crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
+ crystallite_Temperature(g,i,e),g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
+ .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
+ if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ else ! if broken local...
+ crystallite_todo(g,i,e) = .false. ! ... skip this one next time
+ endif
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- STATE UPDATE (EULER INTEGRATION) ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ stateResiduum(1:mySizeDotState,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature
+ temperatureResiduum(g,i,e) = - 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
+ constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
+ + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
+ crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
+ + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- STATE JUMP ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- UPDATE DEPENDENT STATES (EULER INTEGRATION) ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
+ crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ endif
+
+ ! --- STRESS INTEGRATION (EULER INTEGRATION) ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
- .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
- if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- else ! if broken local...
- crystallite_todo(g,i,e) = .false. ! ... skip this one next time
- endif
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- STATE UPDATE (EULER INTEGRATION) ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- stateResiduum(1:mySizeDotState,g,i,e) = - 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature
- temperatureResiduum(g,i,e) = - 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
- constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
- + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
- crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
- + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- STATE JUMP ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- UPDATE DEPENDENT STATES (EULER INTEGRATION) ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
- crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-endif
-
-! --- STRESS INTEGRATION (EULER INTEGRATION) ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-
-if (numerics_integrationMode == 1_pInt) then
-
- ! --- DOT STATE AND TEMPERATURE (HEUN METHOD) ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
- crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
- crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
- crystallite_Temperature(g,i,e),g,i,e)
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
- .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
- if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- else ! if broken local...
- crystallite_todo(g,i,e) = .false. ! ... skip this one next time
- endif
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- ERROR ESTIMATE FOR STATE AND TEMPERATURE (HEUN METHOD) ---
-
- !$OMP SINGLE
- relStateResiduum = 0.0_pReal
- relTemperatureResiduum = 0.0_pReal
- !$OMP END SINGLE
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
-
-
- ! --- contribution of heun step to absolute residui ---
-
- stateResiduum(1:mySizeDotState,g,i,e) = stateResiduum(1:mySizeDotState,g,i,e) &
- + 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature
- temperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) &
- + 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
- !$OMP FLUSH(stateResiduum,temperatureResiduum)
-
- ! --- relative residui ---
-
- forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) &
- relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s)
- if (crystallite_Temperature(g,i,e) > 0_pInt) then
- relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e)
- endif
- !$OMP FLUSH(relStateResiduum,relTemperatureResiduum)
-
+ !$OMP ENDDO
+
+
+ if (numerics_integrationMode == 1_pInt) then
+
+ ! --- DOT STATE AND TEMPERATURE (HEUN METHOD) ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
+ crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
+ crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
+ crystallite_Temperature(g,i,e),g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ if ( any(constitutive_dotState(g,i,e)%p /= constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
+ .or. crystallite_dotTemperature(g,i,e) /= crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
+ if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ else ! if broken local...
+ crystallite_todo(g,i,e) = .false. ! ... skip this one next time
+ endif
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- ERROR ESTIMATE FOR STATE AND TEMPERATURE (HEUN METHOD) ---
+
+ !$OMP SINGLE
+ relStateResiduum = 0.0_pReal
+ relTemperatureResiduum = 0.0_pReal
+ !$OMP END SINGLE
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+
+
+ ! --- contribution of heun step to absolute residui ---
+
+ stateResiduum(1:mySizeDotState,g,i,e) = stateResiduum(1:mySizeDotState,g,i,e) &
+ + 0.5_pReal * constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state and temperature
+ temperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) &
+ + 0.5_pReal * crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
+ !$OMP FLUSH(stateResiduum,temperatureResiduum)
+
+ ! --- relative residui ---
+
+ forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) &
+ relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s)
+ if (crystallite_Temperature(g,i,e) > 0_pInt) then
+ relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e)
+ endif
+ !$OMP FLUSH(relStateResiduum,relTemperatureResiduum)
+
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
- write(6,*)
- write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> absolute residuum tolerance', &
- stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> relative residuum tolerance', &
- relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
- write(6,*)
- write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
- - 2.0_pReal * stateResiduum(1:mySizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum
- write(6,*)
- write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
+ write(6,*)
+ write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> absolute residuum tolerance', &
+ stateResiduum(1:mySizeDotState,g,i,e) / constitutive_aTolState(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ write(6,'(a,/,(12x,12(f12.1,1x)))') '<< CRYST >> relative residuum tolerance', &
+ relStateResiduum(1:mySizeDotState,g,i,e) / rTol_crystalliteState
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState) &
+ - 2.0_pReal * stateResiduum(1:mySizeDotState,g,i,e) / crystallite_subdt(g,i,e) ! calculate former dotstate from higher order solution and state residuum
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.5,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ endif
#endif
-
- ! --- converged ? ---
-
- if ( all( abs(relStateResiduum(:,g,i,e)) < rTol_crystalliteState &
- .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState)) &
- .and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature ) then
- crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- !$OMP CRITICAL (distributionState)
- debug_StateLoopDistribution(2,numerics_integrationMode) = &
- debug_StateLoopDistribution(2,numerics_integrationMode) + 1_pInt
- !$OMP END CRITICAL (distributionState)
- endif
- endif
-
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-elseif (numerics_integrationMode > 1) then ! stiffness calculation
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- !$OMP CRITICAL (distributionState)
- debug_StateLoopDistribution(2,numerics_integrationMode) = &
- debug_StateLoopDistribution(2,numerics_integrationMode) + 1_pInt
- !$OMP END CRITICAL (distributionState)
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-endif
-
-!$OMP END PARALLEL
-
-
-! --- NONLOCAL CONVERGENCE CHECK ---
-
-if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
- !$OMP CRITICAL (write2out)
- write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged'
- write(6,*)
- !$OMP END CRITICAL (write2out)
-endif
-if (.not. singleRun) then ! if not requesting Integration of just a single IP
- if ( any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) then ! any non-local not yet converged (or broken)...
- crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
- endif
-endif
-
+
+ ! --- converged ? ---
+
+ if ( all( abs(relStateResiduum(:,g,i,e)) < rTol_crystalliteState &
+ .or. abs(stateResiduum(1:mySizeDotState,g,i,e)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState)) &
+ .and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature ) then
+ crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ !$OMP CRITICAL (distributionState)
+ debug_StateLoopDistribution(2,numerics_integrationMode) = &
+ debug_StateLoopDistribution(2,numerics_integrationMode) + 1_pInt
+ !$OMP END CRITICAL (distributionState)
+ endif
+ endif
+
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ elseif (numerics_integrationMode > 1) then ! stiffness calculation
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ !$OMP CRITICAL (distributionState)
+ debug_StateLoopDistribution(2,numerics_integrationMode) = &
+ debug_StateLoopDistribution(2,numerics_integrationMode) + 1_pInt
+ !$OMP END CRITICAL (distributionState)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ endif
+
+ !$OMP END PARALLEL
+
+
+ ! --- NONLOCAL CONVERGENCE CHECK ---
+
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
+ !$OMP CRITICAL (write2out)
+ write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged'
+ write(6,*)
+ !$OMP END CRITICAL (write2out)
+ endif
+ if (.not. singleRun) then ! if not requesting Integration of just a single IP
+ if ( any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) then ! any non-local not yet converged (or broken)...
+ crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
+ endif
+ endif
+
end subroutine crystallite_integrateStateAdaptiveEuler
-
-!********************************************************************
-! integrate stress, state and Temperature with
-! 1st order explicit Euler method
-!********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief integrate stress, state and Temperature with 1st order explicit Euler method
+!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateEuler(gg,ii,ee)
-
-!*** variables and functions from other modules ***!
-use numerics, only: numerics_integrationMode, &
- numerics_timeSyncing
-use debug, only: debug_level, &
- debug_crystallite, &
- debug_levelBasic, &
- debug_levelExtensive, &
- debug_levelSelective, &
- debug_e, &
- debug_i, &
- debug_g, &
- debug_StateLoopDistribution
-use FEsolving, only: FEsolving_execElem, &
- FEsolving_execIP
-use mesh, only: mesh_element, &
- mesh_NcpElems
-use material, only: homogenization_Ngrains
-use constitutive, only: constitutive_sizeDotState, &
- constitutive_state, &
- constitutive_subState0, &
- constitutive_dotState, &
- constitutive_collectDotState, &
- constitutive_dotTemperature, &
- constitutive_microstructure
-
-implicit none
-!*** input variables ***!
-integer(pInt), optional, intent(in):: ee, & ! element index
- ii, & ! integration point index
- gg ! grain index
-!*** local variables ***!
-integer(pInt) e, & ! element index in element loop
- i, & ! integration point index in ip loop
- g, & ! grain index in grain loop
- mySizeDotState
-integer(pInt), dimension(2) :: eIter ! bounds for element iteration
-integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
- gIter ! bounds for grain iteration
-logical singleRun ! flag indicating computation for single (g,i,e) triple
-
-
-if (present(ee) .and. present(ii) .and. present(gg)) then
- eIter = ee
- iIter(1:2,ee) = ii
- gIter(1:2,ee) = gg
- singleRun = .true.
-else
- eIter = FEsolving_execElem(1:2)
- do e = eIter(1),eIter(2)
- iIter(1:2,e) = FEsolving_execIP(1:2,e)
- gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
- enddo
- singleRun = .false.
-endif
-
-
-!$OMP PARALLEL
-
-if (numerics_integrationMode == 1_pInt) then
-
- ! --- DOT STATE AND TEMPERATURE ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
- crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
- crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
- crystallite_Temperature(g,i,e),g,i,e)
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
- .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
- if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- else ! if broken local...
- crystallite_todo(g,i,e) = .false. ! ... skip this one next time
- endif
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- UPDATE STATE AND TEMPERATURE ---
-
- !$OMP DO PRIVATE(mySizeDotState)
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
- + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
- crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
- + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
+ use numerics, only: numerics_integrationMode, &
+ numerics_timeSyncing
+ use debug, only: debug_level, &
+ debug_crystallite, &
+ debug_levelBasic, &
+ debug_levelExtensive, &
+ debug_levelSelective, &
+ debug_e, &
+ debug_i, &
+ debug_g, &
+ debug_StateLoopDistribution
+ use FEsolving, only: FEsolving_execElem, &
+ FEsolving_execIP
+ use mesh, only: mesh_element, &
+ mesh_NcpElems
+ use material, only: homogenization_Ngrains
+ use constitutive, only: constitutive_sizeDotState, &
+ constitutive_state, &
+ constitutive_subState0, &
+ constitutive_dotState, &
+ constitutive_collectDotState, &
+ constitutive_dotTemperature, &
+ constitutive_microstructure
+
+ implicit none
+ !*** input variables ***!
+ integer(pInt), optional, intent(in):: ee, & ! element index
+ ii, & ! integration point index
+ gg ! grain index
+ !*** local variables ***!
+ integer(pInt) e, & ! element index in element loop
+ i, & ! integration point index in ip loop
+ g, & ! grain index in grain loop
+ mySizeDotState
+ integer(pInt), dimension(2) :: eIter ! bounds for element iteration
+ integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
+ gIter ! bounds for grain iteration
+ logical singleRun ! flag indicating computation for single (g,i,e) triple
+
+
+ if (present(ee) .and. present(ii) .and. present(gg)) then
+ eIter = ee
+ iIter(1:2,ee) = ii
+ gIter(1:2,ee) = gg
+ singleRun = .true.
+ else
+ eIter = FEsolving_execElem(1:2)
+ do e = eIter(1),eIter(2)
+ iIter(1:2,e) = FEsolving_execIP(1:2,e)
+ gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
+ enddo
+ singleRun = .false.
+ endif
+
+
+ !$OMP PARALLEL
+
+ if (numerics_integrationMode == 1_pInt) then
+
+ ! --- DOT STATE AND TEMPERATURE ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
+ crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
+ crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
+ crystallite_Temperature(g,i,e),g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
+ .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
+ if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ else ! if broken local...
+ crystallite_todo(g,i,e) = .false. ! ... skip this one next time
+ endif
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- UPDATE STATE AND TEMPERATURE ---
+
+ !$OMP DO PRIVATE(mySizeDotState)
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
+ + constitutive_dotState(g,i,e)%p(1:mySizeDotState) * crystallite_subdt(g,i,e)
+ crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
+ + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> update state at el ip g ',e,i,g
- write(6,*)
- write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> update state at el ip g ',e,i,g
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> dotState', constitutive_dotState(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ endif
#endif
- endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- STATE JUMP ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local...
+ .and. .not. numerics_timeSyncing) then
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- UPDATE DEPENDENT STATES ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
+ crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
+ endif
enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- STATE JUMP ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local...
- .and. .not. numerics_timeSyncing) then
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- UPDATE DEPENDENT STATES ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
- crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
- endif
+ !$OMP ENDDO
+
+ endif
+
+
+ ! --- STRESS INTEGRATION ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local...
+ .and. .not. numerics_timeSyncing) then
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
enddo; enddo; enddo
- !$OMP ENDDO
-
-endif
-
-
-! --- STRESS INTEGRATION ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local...
- .and. .not. numerics_timeSyncing) then
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-
-! --- SET CONVERGENCE FLAG ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- crystallite_converged(g,i,e) = .true. ! if still "to do" then converged per definitionem
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- !$OMP CRITICAL (distributionState)
- debug_StateLoopDistribution(1,numerics_integrationMode) = &
- debug_StateLoopDistribution(1,numerics_integrationMode) + 1_pInt
- !$OMP END CRITICAL (distributionState)
- endif
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-!$OMP END PARALLEL
-
-
-! --- CHECK NON-LOCAL CONVERGENCE ---
-
-if (.not. singleRun) then ! if not requesting Integration of just a single IP
- if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity) & ! any non-local not yet converged (or broken)...
- .and. .not. numerics_timeSyncing) then
- crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
- endif
-endif
-
+ !$OMP ENDDO
+
+
+ ! --- SET CONVERGENCE FLAG ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ crystallite_converged(g,i,e) = .true. ! if still "to do" then converged per definitionem
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ !$OMP CRITICAL (distributionState)
+ debug_StateLoopDistribution(1,numerics_integrationMode) = &
+ debug_StateLoopDistribution(1,numerics_integrationMode) + 1_pInt
+ !$OMP END CRITICAL (distributionState)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ !$OMP END PARALLEL
+
+
+ ! --- CHECK NON-LOCAL CONVERGENCE ---
+
+ if (.not. singleRun) then ! if not requesting Integration of just a single IP
+ if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity) & ! any non-local not yet converged (or broken)...
+ .and. .not. numerics_timeSyncing) then
+ crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
+ endif
+ endif
+
end subroutine crystallite_integrateStateEuler
-
-!********************************************************************
-! integrate stress, state and Temperature with
-! adaptive 1st order explicit Euler method
-! using Fixed Point Iteration to adapt the stepsize
-!********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief integrate stress, state and Temperature with adaptive 1st order explicit Euler method
+!> using Fixed Point Iteration to adapt the stepsize
+!--------------------------------------------------------------------------------------------------
subroutine crystallite_integrateStateFPI(gg,ii,ee)
-
-!*** variables and functions from other modules ***!
-use debug, only: debug_e, &
- debug_i, &
- debug_g, &
- debug_level,&
- debug_crystallite, &
- debug_levelBasic, &
- debug_levelExtensive, &
- debug_levelSelective, &
- debug_StateLoopDistribution
-use numerics, only: nState, &
- numerics_integrationMode, &
- rTol_crystalliteState, &
- rTol_crystalliteTemperature
-use FEsolving, only: FEsolving_execElem, &
- FEsolving_execIP
-use mesh, only: mesh_element, &
- mesh_NcpElems
-use material, only: homogenization_Ngrains
-use constitutive, only: constitutive_subState0, &
- constitutive_state, &
- constitutive_sizeDotState, &
- constitutive_maxSizeDotState, &
- constitutive_dotState, &
- constitutive_collectDotState, &
- constitutive_dotTemperature, &
- constitutive_microstructure, &
- constitutive_previousDotState, &
- constitutive_previousDotState2, &
- constitutive_aTolState
-
-implicit none
-
-!*** input variables ***!
-integer(pInt), optional, intent(in):: ee, & ! element index
- ii, & ! integration point index
- gg ! grain index
-
-!*** output variables ***!
-
-!*** local variables ***!
-integer(pInt) NiterationState, & ! number of iterations in state loop
- e, & ! element index in element loop
- i, & ! integration point index in ip loop
- g, & ! grain index in grain loop
- mySizeDotState
-integer(pInt), dimension(2) :: eIter ! bounds for element iteration
-integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
- gIter ! bounds for grain iteration
-real(pReal) dot_prod12, &
- dot_prod22, &
- stateDamper, & ! damper for integration of state
- temperatureResiduum
-real(pReal), dimension(constitutive_maxSizeDotState) :: &
- stateResiduum, &
- tempState
-logical singleRun ! flag indicating computation for single (g,i,e) triple
-
-singleRun = present(ee) .and. present(ii) .and. present(gg)
-if (singleRun) then
- eIter = ee
- iIter(1:2,ee) = ii
- gIter(1:2,ee) = gg
-else
- eIter = FEsolving_execElem(1:2)
- do e = eIter(1),eIter(2)
- iIter(1:2,e) = FEsolving_execIP(1:2,e)
- gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
- enddo
-endif
-
-
-! --+>> PREGUESS FOR STATE <<+--
-
-!$OMP PARALLEL
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- constitutive_previousDotState(g,i,e)%p = 0.0_pReal
- constitutive_previousDotState2(g,i,e)%p = 0.0_pReal
- enddo; enddo; enddo
-!$OMP ENDDO
-
-
-! --- DOT STATES ---
-
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
- crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
- crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
- crystallite_Temperature(g,i,e),g,i,e)
- endif
+ use debug, only: debug_e, &
+ debug_i, &
+ debug_g, &
+ debug_level,&
+ debug_crystallite, &
+ debug_levelBasic, &
+ debug_levelExtensive, &
+ debug_levelSelective, &
+ debug_StateLoopDistribution
+ use numerics, only: nState, &
+ numerics_integrationMode, &
+ rTol_crystalliteState, &
+ rTol_crystalliteTemperature
+ use FEsolving, only: FEsolving_execElem, &
+ FEsolving_execIP
+ use mesh, only: mesh_element, &
+ mesh_NcpElems
+ use material, only: homogenization_Ngrains
+ use constitutive, only: constitutive_subState0, &
+ constitutive_state, &
+ constitutive_sizeDotState, &
+ constitutive_maxSizeDotState, &
+ constitutive_dotState, &
+ constitutive_collectDotState, &
+ constitutive_dotTemperature, &
+ constitutive_microstructure, &
+ constitutive_previousDotState, &
+ constitutive_previousDotState2, &
+ constitutive_aTolState
+
+ implicit none
+ integer(pInt), optional, intent(in):: ee, & ! element index
+ ii, & ! integration point index
+ gg ! grain index
+
+ !*** local variables ***!
+ integer(pInt) NiterationState, & ! number of iterations in state loop
+ e, & ! element index in element loop
+ i, & ! integration point index in ip loop
+ g, & ! grain index in grain loop
+ mySizeDotState
+ integer(pInt), dimension(2) :: eIter ! bounds for element iteration
+ integer(pInt), dimension(2,mesh_NcpElems) :: iIter, & ! bounds for ip iteration
+ gIter ! bounds for grain iteration
+ real(pReal) dot_prod12, &
+ dot_prod22, &
+ stateDamper, & ! damper for integration of state
+ temperatureResiduum
+ real(pReal), dimension(constitutive_maxSizeDotState) :: &
+ stateResiduum, &
+ tempState
+ logical singleRun ! flag indicating computation for single (g,i,e) triple
+
+ singleRun = present(ee) .and. present(ii) .and. present(gg)
+ if (singleRun) then
+ eIter = ee
+ iIter(1:2,ee) = ii
+ gIter(1:2,ee) = gg
+ else
+ eIter = FEsolving_execElem(1:2)
+ do e = eIter(1),eIter(2)
+ iIter(1:2,e) = FEsolving_execIP(1:2,e)
+ gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
+ enddo
+ endif
+
+
+ ! --+>> PREGUESS FOR STATE <<+--
+
+ !$OMP PARALLEL
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ constitutive_previousDotState(g,i,e)%p = 0.0_pReal
+ constitutive_previousDotState2(g,i,e)%p = 0.0_pReal
enddo; enddo; enddo
-!$OMP ENDDO
-!$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e)) then
- if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
- .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
- if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
- !$OMP END CRITICAL (checkTodo)
- else ! broken one was local...
- crystallite_todo(g,i,e) = .false. ! ... done (and broken)
- endif
- endif
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-
-! --- UPDATE STATE AND TEMPERATURE ---
-
-!$OMP DO PRIVATE(mySizeDotState)
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e)) then
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
- + constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e)
- crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
- + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
- endif
- enddo; enddo; enddo
-!$OMP ENDDO
-
-!$OMP END PARALLEL
-
-
-! --+>> STATE LOOP <<+--
-
-NiterationState = 0_pInt
-do while (any(crystallite_todo .and. .not. crystallite_converged) .and. NiterationState < nState ) ! convergence loop for crystallite
- NiterationState = NiterationState + 1_pInt
-
- !$OMP PARALLEL
-
- ! --- UPDATE DEPENDENT STATES ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
- crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
- endif
- constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! remember previous dotState
- constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p ! remember current dotState
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- STRESS INTEGRATION ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ... then all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- !$OMP SINGLE
- !$OMP CRITICAL (write2out)
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
- write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration'
- endif
- !$OMP END CRITICAL (write2out)
- !$OMP END SINGLE
-
-
- ! --- DOT STATE AND TEMPERATURE ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
- crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
- crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
- crystallite_Temperature(g,i,e),g,i,e)
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
- if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
- .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
- crystallite_todo(g,i,e) = .false. ! ... skip me next time
- if (.not. crystallite_localPlasticity(g,i,e)) then ! if me is non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- UPDATE STATE AND TEMPERATURE ---
-
- !$OMP DO PRIVATE(dot_prod12,dot_prod22,statedamper,mySizeDotState,stateResiduum,temperatureResiduum,tempState)
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
-
- ! --- state damper ---
-
- dot_prod12 = dot_product( constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p, &
- constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p )
- dot_prod22 = dot_product( constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p, &
- constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p )
- if ( dot_prod22 > 0.0_pReal &
- .and. ( dot_prod12 < 0.0_pReal &
- .or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) then
- statedamper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
- else
- statedamper = 1.0_pReal
- endif
-
- ! --- get residui ---
-
- mySizeDotState = constitutive_sizeDotState(g,i,e)
- stateResiduum(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
- - constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
- - (constitutive_dotState(g,i,e)%p * statedamper &
- + constitutive_previousDotState(g,i,e)%p * (1.0_pReal - statedamper)) * crystallite_subdt(g,i,e)
- temperatureResiduum = crystallite_Temperature(g,i,e) &
- - crystallite_subTemperature0(g,i,e) &
- - crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
-
- ! --- correct state and temperature with residuum ---
-
- tempState(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) - stateResiduum(1:mySizeDotState) ! need to copy to local variable, since we cant flush a pointer in openmp
- crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - temperatureResiduum
- !$OMP FLUSH(crystallite_Temperature)
+ !$OMP ENDDO
+
+
+ ! --- DOT STATES ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
+ crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
+ crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
+ crystallite_Temperature(g,i,e),g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e)) then
+ if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
+ .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
+ if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken)
+ !$OMP END CRITICAL (checkTodo)
+ else ! broken one was local...
+ crystallite_todo(g,i,e) = .false. ! ... done (and broken)
+ endif
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- UPDATE STATE AND TEMPERATURE ---
+
+ !$OMP DO PRIVATE(mySizeDotState)
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e)) then
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ + constitutive_dotState(g,i,e)%p * crystallite_subdt(g,i,e)
+ crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) &
+ + crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ !$OMP END PARALLEL
+
+
+ ! --+>> STATE LOOP <<+--
+
+ NiterationState = 0_pInt
+ do while (any(crystallite_todo .and. .not. crystallite_converged) .and. NiterationState < nState ) ! convergence loop for crystallite
+ NiterationState = NiterationState + 1_pInt
+
+ !$OMP PARALLEL
+
+ ! --- UPDATE DEPENDENT STATES ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
+ crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
+ endif
+ constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! remember previous dotState
+ constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p ! remember current dotState
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- STRESS INTEGRATION ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ... then all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ !$OMP SINGLE
+ !$OMP CRITICAL (write2out)
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
+ write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration'
+ endif
+ !$OMP END CRITICAL (write2out)
+ !$OMP END SINGLE
+
+
+ ! --- DOT STATE AND TEMPERATURE ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Fe, crystallite_Fp, &
+ crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
+ crystallite_dotTemperature(g,i,e) = constitutive_dotTemperature(crystallite_Tstar_v(1:6,g,i,e), &
+ crystallite_Temperature(g,i,e),g,i,e)
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+ if ( any(constitutive_dotState(g,i,e)%p/=constitutive_dotState(g,i,e)%p) & ! NaN occured in dotState
+ .or. crystallite_dotTemperature(g,i,e)/=crystallite_dotTemperature(g,i,e) ) then ! NaN occured in dotTemperature
+ crystallite_todo(g,i,e) = .false. ! ... skip me next time
+ if (.not. crystallite_localPlasticity(g,i,e)) then ! if me is non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- UPDATE STATE AND TEMPERATURE ---
+
+ !$OMP DO PRIVATE(dot_prod12,dot_prod22,statedamper,mySizeDotState,stateResiduum,temperatureResiduum,tempState)
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then
+
+ ! --- state damper ---
+
+ dot_prod12 = dot_product( constitutive_dotState(g,i,e)%p - constitutive_previousDotState(g,i,e)%p, &
+ constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p )
+ dot_prod22 = dot_product( constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p, &
+ constitutive_previousDotState(g,i,e)%p - constitutive_previousDotState2(g,i,e)%p )
+ if ( dot_prod22 > 0.0_pReal &
+ .and. ( dot_prod12 < 0.0_pReal &
+ .or. dot_product(constitutive_dotState(g,i,e)%p, constitutive_previousDotState(g,i,e)%p) < 0.0_pReal) ) then
+ statedamper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
+ else
+ statedamper = 1.0_pReal
+ endif
+
+ ! --- get residui ---
+
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ stateResiduum(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
+ - constitutive_subState0(g,i,e)%p(1:mySizeDotState) &
+ - (constitutive_dotState(g,i,e)%p * statedamper &
+ + constitutive_previousDotState(g,i,e)%p * (1.0_pReal - statedamper)) * crystallite_subdt(g,i,e)
+ temperatureResiduum = crystallite_Temperature(g,i,e) &
+ - crystallite_subTemperature0(g,i,e) &
+ - crystallite_dotTemperature(g,i,e) * crystallite_subdt(g,i,e)
+
+ ! --- correct state and temperature with residuum ---
+
+ tempState(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) - stateResiduum(1:mySizeDotState) ! need to copy to local variable, since we cant flush a pointer in openmp
+ crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - temperatureResiduum
+ !$OMP FLUSH(crystallite_Temperature)
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> update state at el ip g ',e,i,g
- write(6,*)
- write(6,'(a,f6.1)') '<< CRYST >> statedamper ',statedamper
- write(6,*)
- write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> state residuum',stateResiduum(1:mySizeDotState)
- write(6,*)
- write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state',tempState(1:mySizeDotState)
- write(6,*)
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> update state at el ip g ',e,i,g
+ write(6,*)
+ write(6,'(a,f6.1)') '<< CRYST >> statedamper ',statedamper
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> state residuum',stateResiduum(1:mySizeDotState)
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state',tempState(1:mySizeDotState)
+ write(6,*)
+ endif
#endif
-
- ! --- store corrected dotState --- (cannot do this before state update, because not sure how to flush pointers in openmp)
-
- constitutive_dotState(g,i,e)%p = constitutive_dotState(g,i,e)%p * statedamper &
- + constitutive_previousDotState(g,i,e)%p * (1.0_pReal - statedamper)
-
-
- ! --- converged ? ---
-
- if ( all( abs(stateResiduum(1:mySizeDotState)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) &
- .or. abs(stateResiduum(1:mySizeDotState)) < rTol_crystalliteState &
- * abs(tempState(1:mySizeDotState)) ) &
- .and. (abs(temperatureResiduum) < rTol_crystalliteTemperature * crystallite_Temperature(g,i,e) &
- .or. crystallite_Temperature(g,i,e) == 0.0_pReal) ) then
- crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- !$OMP CRITICAL (distributionState)
- debug_StateLoopDistribution(NiterationState,numerics_integrationMode) = &
- debug_StateLoopDistribution(NiterationState,numerics_integrationMode) + 1_pInt
- !$OMP END CRITICAL (distributionState)
- endif
- endif
- constitutive_state(g,i,e)%p(1:mySizeDotState) = tempState(1:mySizeDotState) ! copy local backup to global pointer
-
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
-
- ! --- STATE JUMP ---
-
- !$OMP DO
- do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
- !$OMP FLUSH(crystallite_todo)
- if (crystallite_todo(g,i,e) .and. crystallite_converged(g,i,e)) then ! converged and still alive...
- crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
- !$OMP FLUSH(crystallite_todo)
- if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken
- crystallite_converged(g,i,e) = .false.
- if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
- !$OMP CRITICAL (checkTodo)
- crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
- !$OMP END CRITICAL (checkTodo)
- endif
- endif
- endif
- enddo; enddo; enddo
- !$OMP ENDDO
-
- !$OMP END PARALLEL
-
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
- !$OMP CRITICAL(write2out)
- write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
- ' grains converged after state integration no. ', NiterationState
- write(6,*)
- !$OMP END CRITICAL(write2out)
- endif
-
-
- ! --- NON-LOCAL CONVERGENCE CHECK ---
-
- if (.not. singleRun) then ! if not requesting Integration of just a single IP
- if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) then ! any non-local not yet converged (or broken)...
- crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
- endif
- endif
-
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
- !$OMP CRITICAL(write2out)
- write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)),' grains converged after non-local check'
- write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after state integration no. ',&
- NiterationState
- write(6,*)
- !$OMP END CRITICAL(write2out)
- endif
-
-enddo ! crystallite convergence loop
-
-
+
+ ! --- store corrected dotState --- (cannot do this before state update, because not sure how to flush pointers in openmp)
+
+ constitutive_dotState(g,i,e)%p = constitutive_dotState(g,i,e)%p * statedamper &
+ + constitutive_previousDotState(g,i,e)%p * (1.0_pReal - statedamper)
+
+
+ ! --- converged ? ---
+
+ if ( all( abs(stateResiduum(1:mySizeDotState)) < constitutive_aTolState(g,i,e)%p(1:mySizeDotState) &
+ .or. abs(stateResiduum(1:mySizeDotState)) < rTol_crystalliteState &
+ * abs(tempState(1:mySizeDotState)) ) &
+ .and. (abs(temperatureResiduum) < rTol_crystalliteTemperature * crystallite_Temperature(g,i,e) &
+ .or. crystallite_Temperature(g,i,e) == 0.0_pReal) ) then
+ crystallite_converged(g,i,e) = .true. ! ... converged per definitionem
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ !$OMP CRITICAL (distributionState)
+ debug_StateLoopDistribution(NiterationState,numerics_integrationMode) = &
+ debug_StateLoopDistribution(NiterationState,numerics_integrationMode) + 1_pInt
+ !$OMP END CRITICAL (distributionState)
+ endif
+ endif
+ constitutive_state(g,i,e)%p(1:mySizeDotState) = tempState(1:mySizeDotState) ! copy local backup to global pointer
+
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+
+ ! --- STATE JUMP ---
+
+ !$OMP DO
+ do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
+ !$OMP FLUSH(crystallite_todo)
+ if (crystallite_todo(g,i,e) .and. crystallite_converged(g,i,e)) then ! converged and still alive...
+ crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e)
+ !$OMP FLUSH(crystallite_todo)
+ if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken
+ crystallite_converged(g,i,e) = .false.
+ if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local...
+ !$OMP CRITICAL (checkTodo)
+ crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals skipped
+ !$OMP END CRITICAL (checkTodo)
+ endif
+ endif
+ endif
+ enddo; enddo; enddo
+ !$OMP ENDDO
+
+ !$OMP END PARALLEL
+
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
+ !$OMP CRITICAL(write2out)
+ write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), &
+ ' grains converged after state integration no. ', NiterationState
+ write(6,*)
+ !$OMP END CRITICAL(write2out)
+ endif
+
+
+ ! --- NON-LOCAL CONVERGENCE CHECK ---
+
+ if (.not. singleRun) then ! if not requesting Integration of just a single IP
+ if (any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) then ! any non-local not yet converged (or broken)...
+ crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged
+ endif
+ endif
+
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then
+ !$OMP CRITICAL(write2out)
+ write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)),' grains converged after non-local check'
+ write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after state integration no. ',&
+ NiterationState
+ write(6,*)
+ !$OMP END CRITICAL(write2out)
+ endif
+
+ enddo ! crystallite convergence loop
+
end subroutine crystallite_integrateStateFPI
-
-!***********************************************************
-!* calculates a jump in the state according to the current *
-!* state and the current stress *
-!***********************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates a jump in the state according to the current state and the current stress
+!--------------------------------------------------------------------------------------------------
function crystallite_stateJump(g,i,e)
-
-!*** variables and functions from other modules ***!
-use debug, only: debug_level, &
- debug_crystallite, &
- debug_levelExtensive, &
- debug_levelSelective, &
- debug_e, &
- debug_i, &
- debug_g
-use FEsolving, only: FEsolving_execElem, &
- FEsolving_execIP
-use mesh, only: mesh_element, &
- mesh_NcpElems
-use material, only: homogenization_Ngrains
-use constitutive, only: constitutive_sizeDotState, &
- constitutive_state, &
- constitutive_deltaState, &
- constitutive_collectDeltaState, &
- constitutive_microstructure
-
-implicit none
-
-!*** input variables ***!
-integer(pInt), intent(in):: e, & ! element index
- i, & ! integration point index
- g ! grain index
-
-!*** output variables ***!
-logical crystallite_stateJump
-
-!*** local variables ***!
-integer(pInt) mySizeDotState
-
-
-crystallite_stateJump = .false.
-
-call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Temperature(g,i,e), g,i,e)
-
-mySizeDotState = constitutive_sizeDotState(g,i,e)
-if (any(constitutive_deltaState(g,i,e)%p(1:mySizeDotState) /= constitutive_deltaState(g,i,e)%p(1:mySizeDotState))) then
- return
-endif
-
-constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
- + constitutive_deltaState(g,i,e)%p(1:mySizeDotState)
-
+ use debug, only: debug_level, &
+ debug_crystallite, &
+ debug_levelExtensive, &
+ debug_levelSelective, &
+ debug_e, &
+ debug_i, &
+ debug_g
+ use FEsolving, only: FEsolving_execElem, &
+ FEsolving_execIP
+ use mesh, only: mesh_element, &
+ mesh_NcpElems
+ use material, only: homogenization_Ngrains
+ use constitutive, only: constitutive_sizeDotState, &
+ constitutive_state, &
+ constitutive_deltaState, &
+ constitutive_collectDeltaState, &
+ constitutive_microstructure
+
+ implicit none
+ integer(pInt), intent(in):: e, & ! element index
+ i, & ! integration point index
+ g ! grain index
+
+ !*** output variables ***!
+ logical crystallite_stateJump
+
+ !*** local variables ***!
+ integer(pInt) mySizeDotState
+
+
+ crystallite_stateJump = .false.
+
+ call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,g,i,e), crystallite_Temperature(g,i,e), g,i,e)
+
+ mySizeDotState = constitutive_sizeDotState(g,i,e)
+ if (any(constitutive_deltaState(g,i,e)%p(1:mySizeDotState) /= constitutive_deltaState(g,i,e)%p(1:mySizeDotState))) then
+ return
+ endif
+
+ constitutive_state(g,i,e)%p(1:mySizeDotState) = constitutive_state(g,i,e)%p(1:mySizeDotState) &
+ + constitutive_deltaState(g,i,e)%p(1:mySizeDotState)
+
#ifndef _OPENMP
-if (any(constitutive_deltaState(g,i,e)%p(1:mySizeDotState) /= 0.0_pReal) &
- .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> update state at el ip g ',e,i,g
- write(6,*)
- write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> deltaState', constitutive_deltaState(g,i,e)%p(1:mySizeDotState)
- write(6,*)
- write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
- write(6,*)
-endif
+ if (any(constitutive_deltaState(g,i,e)%p(1:mySizeDotState) /= 0.0_pReal) &
+ .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> update state at el ip g ',e,i,g
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> deltaState', constitutive_deltaState(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ write(6,'(a,/,(12x,12(e12.6,1x)))') '<< CRYST >> new state', constitutive_state(g,i,e)%p(1:mySizeDotState)
+ write(6,*)
+ endif
#endif
-
-crystallite_stateJump = .true.
-
+
+ crystallite_stateJump = .true.
+
end function crystallite_stateJump
-!***********************************************************************
-!*** calculation of stress (P) with time integration ***
-!*** based on a residuum in Lp and intermediate ***
-!*** acceleration of the Newton-Raphson correction ***
-!***********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief calculation of stress (P) with time integration based on a residuum in Lp and
+!> intermediate acceleration of the Newton-Raphson correction
+!--------------------------------------------------------------------------------------------------
function crystallite_integrateStress(&
- g,& ! grain number
- i,& ! integration point number
- e,& ! element number
- timeFraction &
- )
-
-
-use prec, only: pLongInt
-use numerics, only: nStress, &
- aTol_crystalliteStress, &
- rTol_crystalliteStress, &
- iJacoLpresiduum, &
- numerics_integrationMode
-use debug, only: debug_level, &
- debug_crystallite, &
- debug_levelBasic, &
- debug_levelExtensive, &
- debug_levelSelective, &
- debug_e, &
- debug_i, &
- debug_g, &
- debug_cumLpCalls, &
- debug_cumLpTicks, &
- debug_StressLoopDistribution
-use constitutive, only: constitutive_LpAndItsTangent, &
- constitutive_TandItsTangent, &
- constitutive_homogenizedC
-use math, only: math_mul33x33, &
- math_mul33xx33, &
- math_mul66x6, &
- math_mul99x99, &
- math_transpose33, &
- math_inv33, &
- math_invert33, &
- math_invert, &
- math_det33, &
- math_norm33, &
- math_I3, &
- math_identity2nd, &
- math_Mandel66to3333, &
- math_Mandel6to33, &
- math_Mandel33to6, &
- math_Plain3333to99, &
- math_Plain33to9, &
- math_Plain9to33
-
-implicit none
-
-!*** input variables ***!
-integer(pInt), intent(in):: e, & ! element index
- i, & ! integration point index
- g ! grain index
-real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep
-
-!*** output variables ***!
-logical crystallite_integrateStress ! flag indicating if integration suceeded
-
-!*** local variables ***!
-real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep
- Fp_current, & ! plastic deformation gradient at start of timestep
- Fp_new, & ! plastic deformation gradient at end of timestep
- Fe_new, & ! elastic deformation gradient at end of timestep
- invFp_new, & ! inverse of Fp_new
- invFp_current, & ! inverse of Fp_current
- Lpguess, & ! current guess for plastic velocity gradient
- Lpguess_old, & ! known last good guess for plastic velocity gradient
- Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law
- residuum, & ! current residuum of plastic velocity gradient
- residuum_old, & ! last residuum of plastic velocity gradient
- deltaLp, & ! direction of next guess
- Tstar,& ! 2nd Piola-Kirchhoff Stress
- A,&
- B, &
- Fe ! elastic deformation gradient
-real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
-real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK
-integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK
-real(pReal), dimension(9,9) :: dLp_dT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law
- dT_dFe_constitutive, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law
- dFe_dLp, & ! partial derivative of elastic deformation gradient
- dR_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme)
- dR_dLp2 ! working copy of dRdLp
-real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress
- dFe_dLp3333 ! partial derivative of elastic deformation gradient
-real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress
- det, & ! determinant
- steplength0, &
- steplength, &
- dt, & ! time increment
- aTol
-logical error ! flag indicating an error
-integer(pInt) NiterationStress, & ! number of stress integrations
- ierr, & ! error indicator for LAPACK
- n, &
- o, &
- p, &
- jacoCounter ! counter to check for Jacobian update
-integer(pLongInt) tick, &
- tock, &
- tickrate, &
- maxticks
-
+ g,& ! grain number
+ i,& ! integration point number
+ e,& ! element number
+ timeFraction &
+ )
+ use prec, only: pLongInt
+ use numerics, only: nStress, &
+ aTol_crystalliteStress, &
+ rTol_crystalliteStress, &
+ iJacoLpresiduum, &
+ numerics_integrationMode
+ use debug, only: debug_level, &
+ debug_crystallite, &
+ debug_levelBasic, &
+ debug_levelExtensive, &
+ debug_levelSelective, &
+ debug_e, &
+ debug_i, &
+ debug_g, &
+ debug_cumLpCalls, &
+ debug_cumLpTicks, &
+ debug_StressLoopDistribution
+ use constitutive, only: constitutive_LpAndItsTangent, &
+ constitutive_TandItsTangent, &
+ constitutive_homogenizedC
+ use math, only: math_mul33x33, &
+ math_mul33xx33, &
+ math_mul66x6, &
+ math_mul99x99, &
+ math_transpose33, &
+ math_inv33, &
+ math_invert33, &
+ math_invert, &
+ math_det33, &
+ math_norm33, &
+ math_I3, &
+ math_identity2nd, &
+ math_Mandel66to3333, &
+ math_Mandel6to33, &
+ math_Mandel33to6, &
+ math_Plain3333to99, &
+ math_Plain33to9, &
+ math_Plain9to33
-!* be pessimistic
-
-crystallite_integrateStress = .false.
-#ifndef _OPENMP
-if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress at el ip g ',e,i,g
-endif
-#endif
-
-
-!* only integrate over fraction of timestep?
-
-if (present(timeFraction)) then
- dt = crystallite_subdt(g,i,e) * timeFraction
- Fg_new = crystallite_subF0(1:3,1:3,g,i,e) + (crystallite_subF(1:3,1:3,g,i,e) - crystallite_subF0(1:3,1:3,g,i,e)) * timeFraction
-else
- dt = crystallite_subdt(g,i,e)
- Fg_new = crystallite_subF(1:3,1:3,g,i,e)
-endif
-
+ implicit none
+ integer(pInt), intent(in):: e, & ! element index
+ i, & ! integration point index
+ g ! grain index
+ real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep
-!* feed local variables
-
-Fp_current = crystallite_subFp0(1:3,1:3,g,i,e) ! "Fp_current" is only used as temp var here...
-Lpguess_old = crystallite_Lp(1:3,1:3,g,i,e) ! consider present Lp good (i.e. worth remembering) ...
-Lpguess = crystallite_Lp(1:3,1:3,g,i,e) ! ... and take it as first guess
-
-
-!* inversion of Fp_current...
-
-invFp_current = math_inv33(Fp_current)
-if (all(invFp_current == 0.0_pReal)) then ! ... failed?
-#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el ip g ',e,i,g
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0_pInt) then
- write(6,*)
- write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',math_transpose33(Fp_current(1:3,1:3))
- endif
- endif
-#endif
- return
-endif
-A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
+ !*** output variables ***!
+ logical crystallite_integrateStress ! flag indicating if integration suceeded
-
-!* start LpLoop with normal step length
-
-NiterationStress = 0_pInt
-jacoCounter = 0_pInt
-steplength0 = 1.0_pReal
-steplength = steplength0
-residuum_old = 0.0_pReal
-
-LpLoop: do
- NiterationStress = NiterationStress + 1_pInt
-
+ !*** local variables ***!
+ real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep
+ Fp_current, & ! plastic deformation gradient at start of timestep
+ Fp_new, & ! plastic deformation gradient at end of timestep
+ Fe_new, & ! elastic deformation gradient at end of timestep
+ invFp_new, & ! inverse of Fp_new
+ invFp_current, & ! inverse of Fp_current
+ Lpguess, & ! current guess for plastic velocity gradient
+ Lpguess_old, & ! known last good guess for plastic velocity gradient
+ Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law
+ residuum, & ! current residuum of plastic velocity gradient
+ residuum_old, & ! last residuum of plastic velocity gradient
+ deltaLp, & ! direction of next guess
+ Tstar,& ! 2nd Piola-Kirchhoff Stress
+ A,&
+ B, &
+ Fe ! elastic deformation gradient
+ real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
+ real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK
+ integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK
+ real(pReal), dimension(9,9) :: dLp_dT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law
+ dT_dFe_constitutive, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law
+ dFe_dLp, & ! partial derivative of elastic deformation gradient
+ dR_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme)
+ dR_dLp2 ! working copy of dRdLp
+ real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress
+ dFe_dLp3333 ! partial derivative of elastic deformation gradient
+ real(pReal) p_hydro, & ! volumetric part of 2nd Piola-Kirchhoff Stress
+ det, & ! determinant
+ steplength0, &
+ steplength, &
+ dt, & ! time increment
+ aTol
+ logical error ! flag indicating an error
+ integer(pInt) NiterationStress, & ! number of stress integrations
+ ierr, & ! error indicator for LAPACK
+ n, &
+ o, &
+ p, &
+ jacoCounter ! counter to check for Jacobian update
+ integer(pLongInt) tick, &
+ tock, &
+ tickrate, &
+ maxticks
+ external :: dgesv
- !* too many loops required ?
-
- if (NiterationStress > nStress) then
+ !* be pessimistic
+
+ crystallite_integrateStress = .false.
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- write(6,'(a,i3,a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress reached loop limit',nStress,' at el ip g ',e,i,g
- write(6,*)
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress at el ip g ',e,i,g
+ endif
#endif
- return
- endif
-
-
- !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
-
- B = math_I3 - dt*Lpguess
- Fe = math_mul33x33(A,B) ! current elastic deformation tensor
- call constitutive_TandItsTangent(Tstar, dT_dFe3333, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative
- Tstar_v = math_Mandel33to6(Tstar)
- p_hydro = sum(Tstar_v(1:3)) / 3.0_pReal
- forall(n=1_pInt:3_pInt) Tstar_v(n) = Tstar_v(n) - p_hydro ! get deviatoric stress tensor
-
+
+
+ !* only integrate over fraction of timestep?
+
+ if (present(timeFraction)) then
+ dt = crystallite_subdt(g,i,e) * timeFraction
+ Fg_new = crystallite_subF0(1:3,1:3,g,i,e) + (crystallite_subF(1:3,1:3,g,i,e) - crystallite_subF0(1:3,1:3,g,i,e)) * timeFraction
+ else
+ dt = crystallite_subdt(g,i,e)
+ Fg_new = crystallite_subF(1:3,1:3,g,i,e)
+ endif
+
- !* calculate plastic velocity gradient and its tangent from constitutive law
-
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
- call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
-
- call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e)
-
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
- !$OMP CRITICAL (debugTimingLpTangent)
- debug_cumLpCalls = debug_cumLpCalls + 1_pInt
- debug_cumLpTicks = debug_cumLpTicks + tock-tick
- !$OMP FLUSH (debug_cumLpTicks)
- if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
- !$OMP END CRITICAL (debugTimingLpTangent)
- endif
-
+ !* feed local variables
+
+ Fp_current = crystallite_subFp0(1:3,1:3,g,i,e) ! "Fp_current" is only used as temp var here...
+ Lpguess_old = crystallite_Lp(1:3,1:3,g,i,e) ! consider present Lp good (i.e. worth remembering) ...
+ Lpguess = crystallite_Lp(1:3,1:3,g,i,e) ! ... and take it as first guess
+
+
+ !* inversion of Fp_current...
+
+ invFp_current = math_inv33(Fp_current)
+ if (all(invFp_current == 0.0_pReal)) then ! ... failed?
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,'(a,i3)') '<< CRYST >> iteration ', NiterationStress
- write(6,*)
- write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
- write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el ip g ',e,i,g
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0_pInt) then
+ write(6,*)
+ write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',math_transpose33(Fp_current(1:3,1:3))
+ endif
+ endif
#endif
-
-
- !* update current residuum and check for convergence of loop
+ return
+ endif
+ A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
- aTol = max(rTol_crystalliteStress * max(math_norm33(Lpguess),math_norm33(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
- aTol_crystalliteStress) ! minimum lower cutoff
- residuum = Lpguess - Lp_constitutive
-
- if (any(residuum /= residuum)) then ! NaN in residuum...
+
+ !* start LpLoop with normal step length
+
+ NiterationStress = 0_pInt
+ jacoCounter = 0_pInt
+ steplength0 = 1.0_pReal
+ steplength = steplength0
+ residuum_old = 0.0_pReal
+
+ LpLoop: do
+ NiterationStress = NiterationStress + 1_pInt
+
+
+ !* too many loops required ?
+
+ if (NiterationStress > nStress) then
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- write(6,'(a,i8,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el ip g ',e,i,g,&
- ' ; iteration ', NiterationStress,&
- ' >> returning..!'
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ write(6,'(a,i3,a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress reached loop limit',nStress,' at el ip g ',e,i,g
+ write(6,*)
+ endif
#endif
- return ! ...me = .false. to inform integrator about problem
- elseif (math_norm33(residuum) < aTol) then ! converged if below absolute tolerance
- exit LpLoop ! ...leave iteration loop
- elseif (math_norm33(residuum) < math_norm33(residuum_old) .or. NiterationStress == 1_pInt ) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
- residuum_old = residuum ! ...remember old values and...
- Lpguess_old = Lpguess
- steplength = steplength0 ! ...proceed with normal step length (calculate new search direction)
- else ! not converged and residuum not improved...
- steplength = 0.5_pReal * steplength ! ...try with smaller step length in same direction
- Lpguess = Lpguess_old + steplength * deltaLp
- cycle LpLoop
- endif
-
-
- !* calculate Jacobian for correction term
-
- if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then
- dFe_dLp3333 = 0.0_pReal
- do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
- dFe_dLp3333(p,o,1:3,p) = A(o,1:3) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) delta(j,l)
- enddo; enddo
- dFe_dLp3333 = -dt * dFe_dLp3333
- dFe_dLp = math_Plain3333to99(dFe_dLp3333)
- dT_dFe_constitutive = math_Plain3333to99(dT_dFe3333)
- dR_dLp = math_identity2nd(9_pInt) - &
- math_mul99x99(dLp_dT_constitutive, math_mul99x99(dT_dFe_constitutive , dFe_dLp))
- dR_dLp2 = dR_dLp ! will be overwritten in first call to LAPACK routine
- work = math_plain33to9(residuum)
+ return
+ endif
+
+
+ !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
+
+ B = math_I3 - dt*Lpguess
+ Fe = math_mul33x33(A,B) ! current elastic deformation tensor
+ call constitutive_TandItsTangent(Tstar, dT_dFe3333, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative
+ Tstar_v = math_Mandel33to6(Tstar)
+ p_hydro = sum(Tstar_v(1:3)) / 3.0_pReal
+ forall(n=1_pInt:3_pInt) Tstar_v(n) = Tstar_v(n) - p_hydro ! get deviatoric stress tensor
+
+
+ !* calculate plastic velocity gradient and its tangent from constitutive law
+
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) &
+ call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
+
+ call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e)
+
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
+ !$OMP CRITICAL (debugTimingLpTangent)
+ debug_cumLpCalls = debug_cumLpCalls + 1_pInt
+ debug_cumLpTicks = debug_cumLpTicks + tock-tick
+ !$OMP FLUSH (debug_cumLpTicks)
+ if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
+ !$OMP END CRITICAL (debugTimingLpTangent)
+ endif
+
+#ifndef _OPENMP
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,'(a,i3)') '<< CRYST >> iteration ', NiterationStress
+ write(6,*)
+ write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
+ write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
+ endif
+#endif
+
+
+ !* update current residuum and check for convergence of loop
+
+ aTol = max(rTol_crystalliteStress * max(math_norm33(Lpguess),math_norm33(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error
+ aTol_crystalliteStress) ! minimum lower cutoff
+ residuum = Lpguess - Lp_constitutive
+
+ if (any(residuum /= residuum)) then ! NaN in residuum...
+#ifndef _OPENMP
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ write(6,'(a,i8,1x,i2,1x,i3,a,i3,a)') '<< CRYST >> integrateStress encountered NaN at el ip g ',e,i,g,&
+ ' ; iteration ', NiterationStress,&
+ ' >> returning..!'
+ endif
+#endif
+ return ! ...me = .false. to inform integrator about problem
+ elseif (math_norm33(residuum) < aTol) then ! converged if below absolute tolerance
+ exit LpLoop ! ...leave iteration loop
+ elseif (math_norm33(residuum) < math_norm33(residuum_old) .or. NiterationStress == 1_pInt ) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
+ residuum_old = residuum ! ...remember old values and...
+ Lpguess_old = Lpguess
+ steplength = steplength0 ! ...proceed with normal step length (calculate new search direction)
+ else ! not converged and residuum not improved...
+ steplength = 0.5_pReal * steplength ! ...try with smaller step length in same direction
+ Lpguess = Lpguess_old + steplength * deltaLp
+ cycle LpLoop
+ endif
+
+
+ !* calculate Jacobian for correction term
+
+ if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then
+ dFe_dLp3333 = 0.0_pReal
+ do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
+ dFe_dLp3333(p,o,1:3,p) = A(o,1:3) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) delta(j,l)
+ enddo; enddo
+ dFe_dLp3333 = -dt * dFe_dLp3333
+ dFe_dLp = math_Plain3333to99(dFe_dLp3333)
+ dT_dFe_constitutive = math_Plain3333to99(dT_dFe3333)
+ dR_dLp = math_identity2nd(9_pInt) - &
+ math_mul99x99(dLp_dT_constitutive, math_mul99x99(dT_dFe_constitutive , dFe_dLp))
+ dR_dLp2 = dR_dLp ! will be overwritten in first call to LAPACK routine
+ work = math_plain33to9(residuum)
#if(FLOAT==8)
- call dgesv(9,1,dR_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
+ call dgesv(9,1,dR_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
#elif(FLOAT==4)
- call sgesv(9,1,dR_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
+ call sgesv(9,1,dR_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp
#endif
- if (ierr /= 0_pInt) then
+ if (ierr /= 0_pInt) then
#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,*)
- write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dR_dLp)
- write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(dFe_dLp)
- write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFe_constitutive',transpose(dT_dFe_constitutive)
- write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(dLp_dT_constitutive)
- write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',math_transpose33(A)
- write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',math_transpose33(B)
- write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive',math_transpose33(Lp_constitutive)
- write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess',math_transpose33(Lpguess)
- endif
- endif
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)&
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,*)
+ write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dR_dLp)
+ write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(dFe_dLp)
+ write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFe_constitutive',transpose(dT_dFe_constitutive)
+ write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(dLp_dT_constitutive)
+ write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',math_transpose33(A)
+ write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',math_transpose33(B)
+ write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive',math_transpose33(Lp_constitutive)
+ write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lpguess',math_transpose33(Lpguess)
+ endif
+ endif
#endif
- return
- endif
- deltaLp = - math_plain9to33(work)
- endif
- jacoCounter = jacoCounter + 1_pInt ! increase counter for jaco update
-
- Lpguess = Lpguess + steplength * deltaLp
-
-enddo LpLoop
-
-
-!* calculate new plastic and elastic deformation gradient
-
-invFp_new = math_mul33x33(invFp_current,B)
-invFp_new = invFp_new/math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
-call math_invert33(invFp_new,Fp_new,det,error)
-if (error .or. any(Fp_new/=Fp_new)) then
-#ifndef _OPENMP
- if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',&
- e,i,g, ' ; iteration ', NiterationStress
- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,*)
- write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new)
- endif
- endif
-#endif
- return
-endif
-Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
-
-
-!* add volumetric component to 2nd Piola-Kirchhoff stress and calculate 1st Piola-Kirchhoff stress
-
-forall (n=1_pInt:3_pInt) Tstar_v(n) = Tstar_v(n) + p_hydro
-crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), math_transpose33(invFp_new)))
+ return
+ endif
+ deltaLp = - math_plain9to33(work)
+ endif
+ jacoCounter = jacoCounter + 1_pInt ! increase counter for jaco update
-
-!* store local values in global variables
-
-crystallite_Lp(1:3,1:3,g,i,e) = Lpguess
-crystallite_Tstar_v(1:6,g,i,e) = Tstar_v
-crystallite_Fp(1:3,1:3,g,i,e) = Fp_new
-crystallite_Fe(1:3,1:3,g,i,e) = Fe_new
-crystallite_invFp(1:3,1:3,g,i,e) = invFp_new
-
-
-!* set return flag to true
-
-crystallite_integrateStress = .true.
+ Lpguess = Lpguess + steplength * deltaLp
+
+ enddo LpLoop
+
+
+ !* calculate new plastic and elastic deformation gradient
+
+ invFp_new = math_mul33x33(invFp_current,B)
+ invFp_new = invFp_new/math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
+ call math_invert33(invFp_new,Fp_new,det,error)
+ if (error .or. any(Fp_new/=Fp_new)) then
#ifndef _OPENMP
-if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt &
- .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
- .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
- write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal
- write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Cauchy / MPa', &
- math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose33(Fg_new)) / 1.0e6_pReal / math_det33(Fg_new)
- write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fe Lp Fe^-1', &
- math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order
- write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
-endif
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',&
+ e,i,g, ' ; iteration ', NiterationStress
+ if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,*)
+ write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new)
+ endif
+ endif
#endif
-
-if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
- !$OMP CRITICAL (distributionStress)
- debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) = &
- debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1_pInt
- !$OMP END CRITICAL (distributionStress)
-endif
-
+ return
+ endif
+ Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
+
+
+ !* add volumetric component to 2nd Piola-Kirchhoff stress and calculate 1st Piola-Kirchhoff stress
+
+ forall (n=1_pInt:3_pInt) Tstar_v(n) = Tstar_v(n) + p_hydro
+ crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), math_transpose33(invFp_new)))
+
+
+ !* store local values in global variables
+
+ crystallite_Lp(1:3,1:3,g,i,e) = Lpguess
+ crystallite_Tstar_v(1:6,g,i,e) = Tstar_v
+ crystallite_Fp(1:3,1:3,g,i,e) = Fp_new
+ crystallite_Fe(1:3,1:3,g,i,e) = Fe_new
+ crystallite_invFp(1:3,1:3,g,i,e) = invFp_new
+
+
+ !* set return flag to true
+
+ crystallite_integrateStress = .true.
+#ifndef _OPENMP
+ if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt &
+ .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) &
+ .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
+ write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal
+ write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Cauchy / MPa', &
+ math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose33(Fg_new)) / 1.0e6_pReal / math_det33(Fg_new)
+ write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fe Lp Fe^-1', &
+ math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order
+ write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
+ endif
+#endif
+
+ if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
+ !$OMP CRITICAL (distributionStress)
+ debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) = &
+ debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1_pInt
+ !$OMP END CRITICAL (distributionStress)
+ endif
+
end function crystallite_integrateStress
-
-!********************************************************************
-! calculates orientations and disorientations (in case of single grain ips)
-!********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates orientations and disorientations (in case of single grain ips)
+!--------------------------------------------------------------------------------------------------
subroutine crystallite_orientations
-
-!*** variables and functions from other modules ***!
-
-use math, only: math_pDecomposition, &
- math_RtoQ, &
- math_qDisorientation, &
- math_qConj
-use FEsolving, only: FEsolving_execElem, &
- FEsolving_execIP
-use IO, only: IO_warning
-use material, only: material_phase, &
- homogenization_Ngrains, &
- phase_localPlasticity, &
- phase_plasticityInstance
-use mesh, only: mesh_element, &
- mesh_ipNeighborhood, &
- FE_NipNeighbors, &
- FE_geomtype
-use constitutive_nonlocal, only: constitutive_nonlocal_structure, &
- constitutive_nonlocal_updateCompatibility
-
-implicit none
-
-!*** input variables ***!
-
-!*** output variables ***!
-
-!*** local variables ***!
-integer(pInt) e, & ! element index
- i, & ! integration point index
- g, & ! grain index
- n, & ! neighbor index
- neighboring_e, & ! element index of my neighbor
- neighboring_i, & ! integration point index of my neighbor
- myPhase, & ! phase
- neighboringPhase, &
- myInstance, & ! instance of plasticity
- neighboringInstance, &
- myStructure, & ! lattice structure
- neighboringStructure
-real(pReal), dimension(3,3) :: U, R
-real(pReal), dimension(4) :: orientation
-logical error
-
-! --- CALCULATE ORIENTATION AND LATTICE ROTATION ---
-
-!$OMP PARALLEL DO PRIVATE(error,U,R,orientation)
- do e = FEsolving_execElem(1),FEsolving_execElem(2)
- do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
- do g = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
-
- call math_pDecomposition(crystallite_Fe(1:3,1:3,g,i,e), U, R, error) ! polar decomposition of Fe
- if (error) then
- call IO_warning(650_pInt, e, i, g)
- orientation = [1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! fake orientation
- else
- orientation = math_RtoQ(transpose(R))
- endif
- crystallite_rotation(1:4,g,i,e) = math_qDisorientation(crystallite_orientation0(1:4,g,i,e), & ! active rotation from ori0
- orientation, & ! to current orientation
- 0_pInt ) ! we don't want symmetry here
- crystallite_orientation(1:4,g,i,e) = orientation
- enddo
- enddo
- enddo
-!$OMP END PARALLEL DO
-
-
-! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
-! --- we use crystallite_orientation from above, so need a seperate loop
-
-!$OMP PARALLEL DO PRIVATE(myPhase,myInstance,myStructure,neighboring_e,neighboring_i,neighboringPhase,&
-!$OMP neighboringInstance,neighboringStructure)
- do e = FEsolving_execElem(1),FEsolving_execElem(2)
- do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
- myPhase = material_phase(1,i,e) ! get my phase
- if (.not. phase_localPlasticity(myPhase)) then ! if nonlocal model
- myInstance = phase_plasticityInstance(myPhase)
- myStructure = constitutive_nonlocal_structure(myInstance) ! get my crystal structure
-
-
- ! --- calculate disorientation between me and my neighbor ---
-
- do n = 1_pInt,FE_NipNeighbors(FE_geomtype(mesh_element(2,e))) ! loop through my neighbors
- neighboring_e = mesh_ipNeighborhood(1,n,i,e)
- neighboring_i = mesh_ipNeighborhood(2,n,i,e)
- if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists
- neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase
- if (.not. phase_localPlasticity(neighboringPhase)) then ! neighbor got also nonlocal plasticity
- neighboringInstance = phase_plasticityInstance(neighboringPhase)
- neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure
- if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me
- crystallite_disorientation(:,n,1,i,e) = &
- math_qDisorientation( crystallite_orientation(1:4,1,i,e), &
- crystallite_orientation(1:4,1,neighboring_i,neighboring_e), &
- crystallite_symmetryID(1,i,e)) ! calculate disorientation
- else ! for neighbor with different phase
- crystallite_disorientation(:,n,1,i,e) = (/0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal/) ! 180 degree rotation about 100 axis
- endif
- else ! for neighbor with local plasticity
- crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity
- endif
- else ! no existing neighbor
- crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity
- endif
- enddo
-
-
- ! --- calculate compatibility and transmissivity between me and my neighbor ---
-
- call constitutive_nonlocal_updateCompatibility(crystallite_orientation,i,e)
-
- endif
- enddo
- enddo
-!$OMP END PARALLEL DO
-
+ use math, only: math_pDecomposition, &
+ math_RtoQ, &
+ math_qDisorientation, &
+ math_qConj
+ use FEsolving, only: FEsolving_execElem, &
+ FEsolving_execIP
+ use IO, only: IO_warning
+ use material, only: material_phase, &
+ homogenization_Ngrains, &
+ phase_localPlasticity, &
+ phase_plasticityInstance
+ use mesh, only: mesh_element, &
+ mesh_ipNeighborhood, &
+ FE_NipNeighbors, &
+ FE_geomtype
+ use constitutive_nonlocal, only: constitutive_nonlocal_structure, &
+ constitutive_nonlocal_updateCompatibility
+
+ implicit none
+ integer(pInt) e, & ! element index
+ i, & ! integration point index
+ g, & ! grain index
+ n, & ! neighbor index
+ neighboring_e, & ! element index of my neighbor
+ neighboring_i, & ! integration point index of my neighbor
+ myPhase, & ! phase
+ neighboringPhase, &
+ myInstance, & ! instance of plasticity
+ neighboringInstance, &
+ myStructure, & ! lattice structure
+ neighboringStructure
+ real(pReal), dimension(3,3) :: U, R
+ real(pReal), dimension(4) :: orientation
+ logical error
+
+ ! --- CALCULATE ORIENTATION AND LATTICE ROTATION ---
+
+ !$OMP PARALLEL DO PRIVATE(error,U,R,orientation)
+ do e = FEsolving_execElem(1),FEsolving_execElem(2)
+ do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
+ do g = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
+
+ call math_pDecomposition(crystallite_Fe(1:3,1:3,g,i,e), U, R, error) ! polar decomposition of Fe
+ if (error) then
+ call IO_warning(650_pInt, e, i, g)
+ orientation = [1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! fake orientation
+ else
+ orientation = math_RtoQ(transpose(R))
+ endif
+ crystallite_rotation(1:4,g,i,e) = math_qDisorientation(crystallite_orientation0(1:4,g,i,e), & ! active rotation from ori0
+ orientation, & ! to current orientation
+ 0_pInt ) ! we don't want symmetry here
+ crystallite_orientation(1:4,g,i,e) = orientation
+ enddo
+ enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+
+ ! --- UPDATE SOME ADDITIONAL VARIABLES THAT ARE NEEDED FOR NONLOCAL MATERIAL ---
+ ! --- we use crystallite_orientation from above, so need a seperate loop
+
+ !$OMP PARALLEL DO PRIVATE(myPhase,myInstance,myStructure,neighboring_e,neighboring_i,neighboringPhase,&
+ !$OMP neighboringInstance,neighboringStructure)
+ do e = FEsolving_execElem(1),FEsolving_execElem(2)
+ do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
+ myPhase = material_phase(1,i,e) ! get my phase
+ if (.not. phase_localPlasticity(myPhase)) then ! if nonlocal model
+ myInstance = phase_plasticityInstance(myPhase)
+ myStructure = constitutive_nonlocal_structure(myInstance) ! get my crystal structure
+
+
+ ! --- calculate disorientation between me and my neighbor ---
+
+ do n = 1_pInt,FE_NipNeighbors(FE_geomtype(mesh_element(2,e))) ! loop through my neighbors
+ neighboring_e = mesh_ipNeighborhood(1,n,i,e)
+ neighboring_i = mesh_ipNeighborhood(2,n,i,e)
+ if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists
+ neighboringPhase = material_phase(1,neighboring_i,neighboring_e) ! get my neighbor's phase
+ if (.not. phase_localPlasticity(neighboringPhase)) then ! neighbor got also nonlocal plasticity
+ neighboringInstance = phase_plasticityInstance(neighboringPhase)
+ neighboringStructure = constitutive_nonlocal_structure(neighboringInstance) ! get my neighbor's crystal structure
+ if (myStructure == neighboringStructure) then ! if my neighbor has same crystal structure like me
+ crystallite_disorientation(:,n,1,i,e) = &
+ math_qDisorientation( crystallite_orientation(1:4,1,i,e), &
+ crystallite_orientation(1:4,1,neighboring_i,neighboring_e), &
+ crystallite_symmetryID(1,i,e)) ! calculate disorientation
+ else ! for neighbor with different phase
+ crystallite_disorientation(:,n,1,i,e) = (/0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal/) ! 180 degree rotation about 100 axis
+ endif
+ else ! for neighbor with local plasticity
+ crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity
+ endif
+ else ! no existing neighbor
+ crystallite_disorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! homomorphic identity
+ endif
+ enddo
+
+
+ ! --- calculate compatibility and transmissivity between me and my neighbor ---
+
+ call constitutive_nonlocal_updateCompatibility(crystallite_orientation,i,e)
+
+ endif
+ enddo
+ enddo
+ !$OMP END PARALLEL DO
+
end subroutine crystallite_orientations
-!********************************************************************
-! return results of particular grain
-!********************************************************************
+!--------------------------------------------------------------------------------------------------
+!> @brief return results of particular grain
+!--------------------------------------------------------------------------------------------------
function crystallite_postResults(&
dt,& ! time increment
g,& ! grain number
i,& ! integration point number
e & ! element number
)
-
- !*** variables and functions from other modules ***!
use math, only: math_qToEuler, &
math_qToAxisAngle, &
math_mul33x33, &
@@ -3508,8 +3451,6 @@ function crystallite_postResults(&
constitutive_homogenizedC
implicit none
-
- !*** input variables ***!
integer(pInt), intent(in):: e, & ! element index
i, & ! integration point index
g ! grain index
@@ -3622,7 +3563,4 @@ function crystallite_postResults(&
end function crystallite_postResults
-
-END MODULE
-!##############################################################
-
+end module crystallite