From 42b96354dbf79c592b5b96ca76976bf56e673bed Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 21 Feb 2013 23:08:36 +0000 Subject: [PATCH] doxygen comments --- code/crystallite.f90 | 4580 +++++++++++++++++++++--------------------- 1 file changed, 2259 insertions(+), 2321 deletions(-) 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