diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 2ab66d970..4fa85173d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -9,117 +9,117 @@ !-------------------------------------------------------------------------------------------------- module crystallite - use prec, only: & - pReal - use rotations, only: & - rotation - use FEsolving, only: & - FEsolving_execElem, & - FEsolving_execIP - use material, only: & - homogenization_Ngrains - - implicit none - - private - character(len=64), dimension(:,:), allocatable, private :: & - crystallite_output !< name of each post result output - integer, public, protected :: & - crystallite_maxSizePostResults !< description not available - integer, dimension(:), allocatable, public, protected :: & - crystallite_sizePostResults !< description not available - integer, dimension(:,:), allocatable, private :: & - crystallite_sizePostResult !< description not available - - real(pReal), dimension(:,:,:), allocatable, public :: & - crystallite_dt !< requested time increment of each grain - 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 - type(rotation), dimension(:,:,:), allocatable, private :: & - crystallite_orientation, & !< orientation - crystallite_orientation0 !< initial orientation - real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & - crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - crystallite_P !< 1st Piola-Kirchhoff stress per grain - real(pReal), dimension(:,:,:,:,:), allocatable, public :: & - crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) - crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc - crystallite_partionedS0, & !< 2nd Piola-Kirchhoff stress vector at start of homog inc - 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_Fi, & !< current intermediate def grad (end of converged time step) - crystallite_Fi0, & !< intermediate def grad at start of FE inc - crystallite_partionedFi0,& !< intermediate 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_Li, & !< current intermediate velocitiy grad (end of converged time step) - crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc - crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc - real(pReal), dimension(:,:,:,:,:), allocatable, private :: & - crystallite_subS0, & !< 2nd Piola-Kirchhoff stress vector 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_invFi, & !< inverse of current intermediate def grad (end of converged time step) - crystallite_subFi0,& !< intermediate 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_subLi0 !< intermediate velocity grad at start of crystallite inc - real(pReal), dimension(:,:,:,:,:,:,:), allocatable, public :: & - crystallite_dPdF !< current individual dPdF per grain (end of converged time step) - logical, dimension(:,:,:), allocatable, public :: & - crystallite_requested !< used by upper level (homogenization) to request crystallite calculation - logical, dimension(:,:,:), allocatable, private :: & - crystallite_converged, & !< convergence flag - crystallite_todo, & !< flag to indicate need for further computation - crystallite_localPlasticity !< indicates this grain to have purely local constitutive law - - enum, bind(c) - enumerator :: undefined_ID, & - phase_ID, & - texture_ID, & - volume_ID, & - orientation_ID, & - grainrotation_ID, & - defgrad_ID, & - fe_ID, & - fp_ID, & - fi_ID, & - lp_ID, & - li_ID, & - p_ID, & - s_ID, & - elasmatrix_ID, & - neighboringip_ID, & - neighboringelement_ID - end enum - integer(kind(undefined_ID)),dimension(:,:), allocatable, private :: & - crystallite_outputID !< ID of each post result output - procedure(), pointer :: integrateState - - public :: & - crystallite_init, & - crystallite_stress, & - crystallite_stressTangent, & - crystallite_orientations, & - crystallite_push33ToRef, & - crystallite_postResults - private :: & - integrateStress, & - integrateState, & - integrateStateFPI, & - integrateStateEuler, & - integrateStateAdaptiveEuler, & - integrateStateRK4, & - integrateStateRKCK45, & - stateJump + use prec, only: & + pReal + use rotations, only: & + rotation + use FEsolving, only: & + FEsolving_execElem, & + FEsolving_execIP + use material, only: & + homogenization_Ngrains + + implicit none + + private + character(len=64), dimension(:,:), allocatable, private :: & + crystallite_output !< name of each post result output + integer, public, protected :: & + crystallite_maxSizePostResults !< description not available + integer, dimension(:), allocatable, public, protected :: & + crystallite_sizePostResults !< description not available + integer, dimension(:,:), allocatable, private :: & + crystallite_sizePostResult !< description not available + + real(pReal), dimension(:,:,:), allocatable, public :: & + crystallite_dt !< requested time increment of each grain + 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 + type(rotation), dimension(:,:,:), allocatable, private :: & + crystallite_orientation, & !< orientation + crystallite_orientation0 !< initial orientation + real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & + crystallite_Fe, & !< current "elastic" def grad (end of converged time step) + crystallite_P !< 1st Piola-Kirchhoff stress per grain + real(pReal), dimension(:,:,:,:,:), allocatable, public :: & + crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) + crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc + crystallite_partionedS0, & !< 2nd Piola-Kirchhoff stress vector at start of homog inc + 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_Fi, & !< current intermediate def grad (end of converged time step) + crystallite_Fi0, & !< intermediate def grad at start of FE inc + crystallite_partionedFi0,& !< intermediate 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_Li, & !< current intermediate velocitiy grad (end of converged time step) + crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc + crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc + real(pReal), dimension(:,:,:,:,:), allocatable, private :: & + crystallite_subS0, & !< 2nd Piola-Kirchhoff stress vector 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_invFi, & !< inverse of current intermediate def grad (end of converged time step) + crystallite_subFi0,& !< intermediate 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_subLi0 !< intermediate velocity grad at start of crystallite inc + real(pReal), dimension(:,:,:,:,:,:,:), allocatable, public :: & + crystallite_dPdF !< current individual dPdF per grain (end of converged time step) + logical, dimension(:,:,:), allocatable, public :: & + crystallite_requested !< used by upper level (homogenization) to request crystallite calculation + logical, dimension(:,:,:), allocatable, private :: & + crystallite_converged, & !< convergence flag + crystallite_todo, & !< flag to indicate need for further computation + crystallite_localPlasticity !< indicates this grain to have purely local constitutive law + + enum, bind(c) + enumerator :: undefined_ID, & + phase_ID, & + texture_ID, & + volume_ID, & + orientation_ID, & + grainrotation_ID, & + defgrad_ID, & + fe_ID, & + fp_ID, & + fi_ID, & + lp_ID, & + li_ID, & + p_ID, & + s_ID, & + elasmatrix_ID, & + neighboringip_ID, & + neighboringelement_ID + end enum + integer(kind(undefined_ID)),dimension(:,:), allocatable, private :: & + crystallite_outputID !< ID of each post result output + procedure(), pointer :: integrateState + + public :: & + crystallite_init, & + crystallite_stress, & + crystallite_stressTangent, & + crystallite_orientations, & + crystallite_push33ToRef, & + crystallite_postResults + private :: & + integrateStress, & + integrateState, & + integrateStateFPI, & + integrateStateEuler, & + integrateStateAdaptiveEuler, & + integrateStateRK4, & + integrateStateRKCK45, & + stateJump contains @@ -129,274 +129,274 @@ contains !-------------------------------------------------------------------------------------------------- subroutine crystallite_init #ifdef DEBUG - use debug, only: & - debug_info, & - debug_reset, & - debug_level, & - debug_crystallite, & - debug_levelBasic + use debug, only: & + debug_info, & + debug_reset, & + debug_level, & + debug_crystallite, & + debug_levelBasic #endif - use numerics, only: & - numerics_integrator, & - worldrank, & - usePingPong - use math, only: & - math_I3, & - math_EulerToR, & - math_inv33 - use mesh, only: & - theMesh, & - mesh_element - use IO, only: & - IO_stringValue, & - IO_write_jobFile, & - IO_error - use material - use config, only: & - config_deallocate, & - config_crystallite, & - crystallite_name - use constitutive, only: & - constitutive_initialFi, & - constitutive_microstructure ! derived (shortcut) quantities of given state - - implicit none - - integer, parameter :: FILEUNIT=434 - logical, dimension(:,:), allocatable :: devNull - integer :: & - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e, & !< counter in element loop - o = 0, & !< counter in output loop - r, & - cMax, & !< maximum number of integration point components - iMax, & !< maximum number of integration points - eMax, & !< maximum number of elements - myNcomponents, & !< number of components at current IP - mySize - - character(len=65536), dimension(:), allocatable :: str - - write(6,'(/,a)') ' <<<+- crystallite init -+>>>' - - cMax = homogenization_maxNgrains - iMax = theMesh%elem%nIPs - eMax = theMesh%nElems - - allocate(crystallite_S0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_partionedS0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_S(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subS0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_P(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_F0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_partionedF0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_partionedF(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subF0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subF(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_Fp0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_partionedFp0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subFp0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_Fp(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_invFp(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_Fi0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_partionedFi0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subFi0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_Fi(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_invFi(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_Fe(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_Lp0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_partionedLp0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subLp0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_Lp(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_Li0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_partionedLi0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subLi0(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_Li(3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_dt(cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subdt(cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subFrac(cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_subStep(cMax,iMax,eMax), source=0.0_pReal) - allocate(crystallite_orientation(cMax,iMax,eMax)) - allocate(crystallite_orientation0(cMax,iMax,eMax)) - allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) - allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) - allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) - allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) - allocate(crystallite_output(maxval(crystallite_Noutput), & - size(config_crystallite))) ; crystallite_output = '' - allocate(crystallite_outputID(maxval(crystallite_Noutput), & - size(config_crystallite)), source=undefined_ID) - allocate(crystallite_sizePostResults(size(config_crystallite)),source=0) - allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & - size(config_crystallite)), source=0) - - select case(numerics_integrator) - case(1) - integrateState => integrateStateFPI - case(2) - integrateState => integrateStateEuler - case(3) - integrateState => integrateStateAdaptiveEuler - case(4) - integrateState => integrateStateRK4 - case(5) - integrateState => integrateStateRKCK45 - end select - - - - do c = 1, size(config_crystallite) + use numerics, only: & + numerics_integrator, & + worldrank, & + usePingPong + use math, only: & + math_I3, & + math_EulerToR, & + math_inv33 + use mesh, only: & + theMesh, & + mesh_element + use IO, only: & + IO_stringValue, & + IO_write_jobFile, & + IO_error + use material + use config, only: & + config_deallocate, & + config_crystallite, & + crystallite_name + use constitutive, only: & + constitutive_initialFi, & + constitutive_microstructure ! derived (shortcut) quantities of given state + + implicit none + + integer, parameter :: FILEUNIT=434 + logical, dimension(:,:), allocatable :: devNull + integer :: & + c, & !< counter in integration point component loop + i, & !< counter in integration point loop + e, & !< counter in element loop + o = 0, & !< counter in output loop + r, & + cMax, & !< maximum number of integration point components + iMax, & !< maximum number of integration points + eMax, & !< maximum number of elements + myNcomponents, & !< number of components at current IP + mySize + + character(len=65536), dimension(:), allocatable :: str + + write(6,'(/,a)') ' <<<+- crystallite init -+>>>' + + cMax = homogenization_maxNgrains + iMax = theMesh%elem%nIPs + eMax = theMesh%nElems + + allocate(crystallite_S0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedS0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_S(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subS0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_P(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_F0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedF0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedF(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subF0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subF(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Fp0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedFp0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subFp0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Fp(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_invFp(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Fi0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedFi0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subFi0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Fi(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_invFi(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Fe(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Lp0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedLp0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subLp0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Lp(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Li0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_partionedLi0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subLi0(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_Li(3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_dt(cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subdt(cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subFrac(cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_subStep(cMax,iMax,eMax), source=0.0_pReal) + allocate(crystallite_orientation(cMax,iMax,eMax)) + allocate(crystallite_orientation0(cMax,iMax,eMax)) + allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) + allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) + allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) + allocate(crystallite_converged(cMax,iMax,eMax), source=.true.) + allocate(crystallite_output(maxval(crystallite_Noutput), & + size(config_crystallite))) ; crystallite_output = '' + allocate(crystallite_outputID(maxval(crystallite_Noutput), & + size(config_crystallite)), source=undefined_ID) + allocate(crystallite_sizePostResults(size(config_crystallite)),source=0) + allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), & + size(config_crystallite)), source=0) + + select case(numerics_integrator) + case(1) + integrateState => integrateStateFPI + case(2) + integrateState => integrateStateEuler + case(3) + integrateState => integrateStateAdaptiveEuler + case(4) + integrateState => integrateStateRK4 + case(5) + integrateState => integrateStateRKCK45 + end select + + + + do c = 1, size(config_crystallite) #if defined(__GFORTRAN__) - str = ['GfortranBug86277'] - str = config_crystallite(c)%getStrings('(output)',defaultVal=str) - if (str(1) == 'GfortranBug86277') str = [character(len=65536)::] + str = ['GfortranBug86277'] + str = config_crystallite(c)%getStrings('(output)',defaultVal=str) + if (str(1) == 'GfortranBug86277') str = [character(len=65536)::] #else - str = config_crystallite(c)%getStrings('(output)',defaultVal=[character(len=65536)::]) + str = config_crystallite(c)%getStrings('(output)',defaultVal=[character(len=65536)::]) #endif - do o = 1, size(str) - crystallite_output(o,c) = str(o) - outputName: select case(str(o)) - case ('phase') outputName - crystallite_outputID(o,c) = phase_ID - case ('texture') outputName - crystallite_outputID(o,c) = texture_ID - case ('volume') outputName - crystallite_outputID(o,c) = volume_ID - case ('orientation') outputName - crystallite_outputID(o,c) = orientation_ID - case ('grainrotation') outputName - crystallite_outputID(o,c) = grainrotation_ID - case ('defgrad','f') outputName ! ToDo: no alias (f only) - crystallite_outputID(o,c) = defgrad_ID - case ('fe') outputName - crystallite_outputID(o,c) = fe_ID - case ('fp') outputName - crystallite_outputID(o,c) = fp_ID - case ('fi') outputName - crystallite_outputID(o,c) = fi_ID - case ('lp') outputName - crystallite_outputID(o,c) = lp_ID - case ('li') outputName - crystallite_outputID(o,c) = li_ID - case ('p','firstpiola','1stpiola') outputName ! ToDo: no alias (p only) - crystallite_outputID(o,c) = p_ID - case ('s','tstar','secondpiola','2ndpiola') outputName ! ToDo: no alias (s only) - crystallite_outputID(o,c) = s_ID - case ('elasmatrix') outputName - crystallite_outputID(o,c) = elasmatrix_ID - case ('neighboringip') outputName ! ToDo: this is not a result, it is static. Should be written out by mesh - crystallite_outputID(o,c) = neighboringip_ID - case ('neighboringelement') outputName ! ToDo: this is not a result, it is static. Should be written out by mesh - crystallite_outputID(o,c) = neighboringelement_ID - case default outputName - call IO_error(105,ext_msg=trim(str(o))//' (Crystallite)') - end select outputName - enddo - enddo + do o = 1, size(str) + crystallite_output(o,c) = str(o) + outputName: select case(str(o)) + case ('phase') outputName + crystallite_outputID(o,c) = phase_ID + case ('texture') outputName + crystallite_outputID(o,c) = texture_ID + case ('volume') outputName + crystallite_outputID(o,c) = volume_ID + case ('orientation') outputName + crystallite_outputID(o,c) = orientation_ID + case ('grainrotation') outputName + crystallite_outputID(o,c) = grainrotation_ID + case ('defgrad','f') outputName ! ToDo: no alias (f only) + crystallite_outputID(o,c) = defgrad_ID + case ('fe') outputName + crystallite_outputID(o,c) = fe_ID + case ('fp') outputName + crystallite_outputID(o,c) = fp_ID + case ('fi') outputName + crystallite_outputID(o,c) = fi_ID + case ('lp') outputName + crystallite_outputID(o,c) = lp_ID + case ('li') outputName + crystallite_outputID(o,c) = li_ID + case ('p','firstpiola','1stpiola') outputName ! ToDo: no alias (p only) + crystallite_outputID(o,c) = p_ID + case ('s','tstar','secondpiola','2ndpiola') outputName ! ToDo: no alias (s only) + crystallite_outputID(o,c) = s_ID + case ('elasmatrix') outputName + crystallite_outputID(o,c) = elasmatrix_ID + case ('neighboringip') outputName ! ToDo: this is not a result, it is static. Should be written out by mesh + crystallite_outputID(o,c) = neighboringip_ID + case ('neighboringelement') outputName ! ToDo: this is not a result, it is static. Should be written out by mesh + crystallite_outputID(o,c) = neighboringelement_ID + case default outputName + call IO_error(105,ext_msg=trim(str(o))//' (Crystallite)') + end select outputName + enddo + enddo - do r = 1,size(config_crystallite) - do o = 1,crystallite_Noutput(r) - select case(crystallite_outputID(o,r)) - case(phase_ID,texture_ID,volume_ID) - mySize = 1 - case(orientation_ID,grainrotation_ID) - mySize = 4 - case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,p_ID,s_ID) - mySize = 9 - case(elasmatrix_ID) - mySize = 36 - case(neighboringip_ID,neighboringelement_ID) - mySize = theMesh%elem%nIPneighbors - case default - mySize = 0 - end select - crystallite_sizePostResult(o,r) = mySize - crystallite_sizePostResults(r) = crystallite_sizePostResults(r) + mySize - enddo - enddo - - crystallite_maxSizePostResults = & - maxval(crystallite_sizePostResults(microstructure_crystallite),microstructure_active) + do r = 1,size(config_crystallite) + do o = 1,crystallite_Noutput(r) + select case(crystallite_outputID(o,r)) + case(phase_ID,texture_ID,volume_ID) + mySize = 1 + case(orientation_ID,grainrotation_ID) + mySize = 4 + case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,p_ID,s_ID) + mySize = 9 + case(elasmatrix_ID) + mySize = 36 + case(neighboringip_ID,neighboringelement_ID) + mySize = theMesh%elem%nIPneighbors + case default + mySize = 0 + end select + crystallite_sizePostResult(o,r) = mySize + crystallite_sizePostResults(r) = crystallite_sizePostResults(r) + mySize + enddo + enddo + + crystallite_maxSizePostResults = & + maxval(crystallite_sizePostResults(microstructure_crystallite),microstructure_active) !-------------------------------------------------------------------------------------------------- ! write description file for crystallite output - if (worldrank == 0) then - call IO_write_jobFile(FILEUNIT,'outputCrystallite') - - do r = 1,size(config_crystallite) - if (any(microstructure_crystallite(mesh_element(4,:)) == r)) then - write(FILEUNIT,'(/,a,/)') '['//trim(crystallite_name(r))//']' - do o = 1,crystallite_Noutput(r) - write(FILEUNIT,'(a,i4)') trim(crystallite_output(o,r))//char(9),crystallite_sizePostResult(o,r) - enddo - endif - enddo - - close(FILEUNIT) - endif - - call config_deallocate('material.config/crystallite') + if (worldrank == 0) then + call IO_write_jobFile(FILEUNIT,'outputCrystallite') + + do r = 1,size(config_crystallite) + if (any(microstructure_crystallite(mesh_element(4,:)) == r)) then + write(FILEUNIT,'(/,a,/)') '['//trim(crystallite_name(r))//']' + do o = 1,crystallite_Noutput(r) + write(FILEUNIT,'(a,i4)') trim(crystallite_output(o,r))//char(9),crystallite_sizePostResult(o,r) + enddo + endif + enddo + + close(FILEUNIT) + endif + + call config_deallocate('material.config/crystallite') !-------------------------------------------------------------------------------------------------- ! initialize !$OMP PARALLEL DO PRIVATE(myNcomponents,i,c) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents) - crystallite_Fp0(1:3,1:3,c,i,e) = math_EulerToR(material_EulerAngles(1:3,c,i,e)) ! plastic def gradient reflects init orientation - crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e) - crystallite_F0(1:3,1:3,c,i,e) = math_I3 - crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phase(c,i,e)) - crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), & - crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration - crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) - crystallite_Fi(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e) - crystallite_requested(c,i,e) = .true. - endforall - enddo - !$OMP END PARALLEL DO - - if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? - - crystallite_partionedFp0 = crystallite_Fp0 - crystallite_partionedFi0 = crystallite_Fi0 - crystallite_partionedF0 = crystallite_F0 - crystallite_partionedF = crystallite_F0 - - call crystallite_orientations() - crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations - - !$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1,homogenization_Ngrains(mesh_element(3,e)) - call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), & - crystallite_Fp(1:3,1:3,c,i,e), & - c,i,e) ! update dependent state variables to be consistent with basic states + do e = FEsolving_execElem(1),FEsolving_execElem(2) + myNcomponents = homogenization_Ngrains(mesh_element(3,e)) + forall (i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), c = 1:myNcomponents) + crystallite_Fp0(1:3,1:3,c,i,e) = math_EulerToR(material_EulerAngles(1:3,c,i,e)) ! plastic def gradient reflects init orientation + crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e) + crystallite_F0(1:3,1:3,c,i,e) = math_I3 + crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phase(c,i,e)) + crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), & + crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration + crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) + crystallite_Fi(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e) + crystallite_requested(c,i,e) = .true. + endforall + enddo + !$OMP END PARALLEL DO + + if(any(.not. crystallite_localPlasticity) .and. .not. usePingPong) call IO_error(601) ! exit if nonlocal but no ping-pong ToDo: Why not check earlier? or in nonlocal? + + crystallite_partionedFp0 = crystallite_Fp0 + crystallite_partionedFi0 = crystallite_Fi0 + crystallite_partionedF0 = crystallite_F0 + crystallite_partionedF = crystallite_F0 + + call crystallite_orientations() + crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations + + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + do c = 1,homogenization_Ngrains(mesh_element(3,e)) + call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), & + crystallite_Fp(1:3,1:3,c,i,e), & + c,i,e) ! update dependent state variables to be consistent with basic states + enddo enddo - enddo - enddo - !$OMP END PARALLEL DO - - devNull = crystallite_stress() - call crystallite_stressTangent + enddo + !$OMP END PARALLEL DO + + devNull = crystallite_stress() + call crystallite_stressTangent #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then - write(6,'(a42,1x,i10)') ' # of elements: ', eMax - write(6,'(a42,1x,i10)') 'max # of integration points/element: ', iMax - write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax - write(6,'(a42,1x,i10)') 'max # of neigbours/integration point: ', theMesh%elem%nIPneighbors - write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) - flush(6) - endif - - call debug_info - call debug_reset + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then + write(6,'(a42,1x,i10)') ' # of elements: ', eMax + write(6,'(a42,1x,i10)') 'max # of integration points/element: ', iMax + write(6,'(a42,1x,i10)') 'max # of constituents/integration point: ', cMax + write(6,'(a42,1x,i10)') 'max # of neigbours/integration point: ', theMesh%elem%nIPneighbors + write(6,'(a42,1x,i10)') ' # of nonlocal constituents: ',count(.not. crystallite_localPlasticity) + flush(6) + endif + + call debug_info + call debug_reset #endif end subroutine crystallite_init @@ -406,270 +406,270 @@ end subroutine crystallite_init !> @brief calculate stress (P) !-------------------------------------------------------------------------------------------------- function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC) - use prec, only: & - tol_math_check, & - dNeq0 - use numerics, only: & - subStepMinCryst, & - subStepSizeCryst, & - stepIncreaseCryst + use prec, only: & + tol_math_check, & + dNeq0 + use numerics, only: & + subStepMinCryst, & + subStepSizeCryst, & + stepIncreaseCryst #ifdef DEBUG - use debug, only: & - debug_level, & - debug_crystallite, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_e, & - debug_i, & - debug_g + use debug, only: & + debug_level, & + debug_crystallite, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective, & + debug_e, & + debug_i, & + debug_g #endif - use IO, only: & - IO_warning, & - IO_error - use math, only: & - math_inv33 - use mesh, only: & - theMesh, & - mesh_element - use material, only: & - homogenization_Ngrains, & - plasticState, & - sourceState, & - phase_Nsources, & - phaseAt, phasememberAt - - implicit none - logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress - real(pReal), intent(in), optional :: & - dummyArgumentToPreventInternalCompilerErrorWithGCC - real(pReal) :: & - formerSubStep - integer :: & - NiterationCrystallite, & ! number of iterations in crystallite loop - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e, & !< counter in element loop - startIP, endIP, & - s - + use IO, only: & + IO_warning, & + IO_error + use math, only: & + math_inv33 + use mesh, only: & + theMesh, & + mesh_element + use material, only: & + homogenization_Ngrains, & + plasticState, & + sourceState, & + phase_Nsources, & + phaseAt, phasememberAt + + implicit none + logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress + real(pReal), intent(in), optional :: & + dummyArgumentToPreventInternalCompilerErrorWithGCC + real(pReal) :: & + formerSubStep + integer :: & + NiterationCrystallite, & ! number of iterations in crystallite loop + c, & !< counter in integration point component loop + i, & !< counter in integration point loop + e, & !< counter in element loop + startIP, endIP, & + s + #ifdef DEBUG - if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & - .and. FEsolving_execElem(1) <= debug_e & - .and. debug_e <= FEsolving_execElem(2)) then - write(6,'(/,a,i8,1x,i2,1x,i3)') '<< CRYST stress >> boundary and initial values at el ip ipc ', & - debug_e,debug_i, debug_g - write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> F ', & - transpose(crystallite_partionedF(1:3,1:3,debug_g,debug_i,debug_e)) - write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> F0 ', & - transpose(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e)) - write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Fp0', & - transpose(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e)) - write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Fi0', & - transpose(crystallite_partionedFi0(1:3,1:3,debug_g,debug_i,debug_e)) - write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Lp0', & - transpose(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e)) - write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Li0', & - transpose(crystallite_partionedLi0(1:3,1:3,debug_g,debug_i,debug_e)) - endif + if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0 & + .and. FEsolving_execElem(1) <= debug_e & + .and. debug_e <= FEsolving_execElem(2)) then + write(6,'(/,a,i8,1x,i2,1x,i3)') '<< CRYST stress >> boundary and initial values at el ip ipc ', & + debug_e,debug_i, debug_g + write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> F ', & + transpose(crystallite_partionedF(1:3,1:3,debug_g,debug_i,debug_e)) + write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> F0 ', & + transpose(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e)) + write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Fp0', & + transpose(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e)) + write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Fi0', & + transpose(crystallite_partionedFi0(1:3,1:3,debug_g,debug_i,debug_e)) + write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Lp0', & + transpose(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e)) + write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Li0', & + transpose(crystallite_partionedLi0(1:3,1:3,debug_g,debug_i,debug_e)) + endif #endif !-------------------------------------------------------------------------------------------------- ! initialize to starting condition - crystallite_subStep = 0.0_pReal - !$OMP PARALLEL DO - elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,homogenization_Ngrains(mesh_element(3,e)) - homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then - plasticState (phaseAt(c,i,e))%subState0( :,phasememberAt(c,i,e)) = & - plasticState (phaseAt(c,i,e))%partionedState0(:,phasememberAt(c,i,e)) + crystallite_subStep = 0.0_pReal + !$OMP PARALLEL DO + elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1,homogenization_Ngrains(mesh_element(3,e)) + homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then + plasticState (phaseAt(c,i,e))%subState0( :,phasememberAt(c,i,e)) = & + plasticState (phaseAt(c,i,e))%partionedState0(:,phasememberAt(c,i,e)) - do s = 1, phase_Nsources(phaseAt(c,i,e)) - sourceState(phaseAt(c,i,e))%p(s)%subState0( :,phasememberAt(c,i,e)) = & - sourceState(phaseAt(c,i,e))%p(s)%partionedState0(:,phasememberAt(c,i,e)) - enddo - crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e) - crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_partionedLp0(1:3,1:3,c,i,e) - crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partionedFi0(1:3,1:3,c,i,e) - crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_partionedLi0(1:3,1:3,c,i,e) - crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) - crystallite_subS0(1:3,1:3,c,i,e) = crystallite_partionedS0(1:3,1:3,c,i,e) - crystallite_subFrac(c,i,e) = 0.0_pReal - crystallite_subStep(c,i,e) = 1.0_pReal/subStepSizeCryst - crystallite_todo(c,i,e) = .true. - crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst - endif homogenizationRequestsCalculation - enddo; enddo - enddo elementLooping1 - !$OMP END PARALLEL DO + do s = 1, phase_Nsources(phaseAt(c,i,e)) + sourceState(phaseAt(c,i,e))%p(s)%subState0( :,phasememberAt(c,i,e)) = & + sourceState(phaseAt(c,i,e))%p(s)%partionedState0(:,phasememberAt(c,i,e)) + enddo + crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e) + crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_partionedLp0(1:3,1:3,c,i,e) + crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partionedFi0(1:3,1:3,c,i,e) + crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_partionedLi0(1:3,1:3,c,i,e) + crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) + crystallite_subS0(1:3,1:3,c,i,e) = crystallite_partionedS0(1:3,1:3,c,i,e) + crystallite_subFrac(c,i,e) = 0.0_pReal + crystallite_subStep(c,i,e) = 1.0_pReal/subStepSizeCryst + crystallite_todo(c,i,e) = .true. + crystallite_converged(c,i,e) = .false. ! pretend failed step of 1/subStepSizeCryst + endif homogenizationRequestsCalculation + enddo; enddo + enddo elementLooping1 + !$OMP END PARALLEL DO - singleRun: if (FEsolving_execELem(1) == FEsolving_execElem(2) .and. & - FEsolving_execIP(1,FEsolving_execELem(1))==FEsolving_execIP(2,FEsolving_execELem(1))) then - startIP = FEsolving_execIP(1,FEsolving_execELem(1)) - endIP = startIP - else singleRun - startIP = 1 - endIP = theMesh%elem%nIPs - endif singleRun + singleRun: if (FEsolving_execELem(1) == FEsolving_execElem(2) .and. & + FEsolving_execIP(1,FEsolving_execELem(1))==FEsolving_execIP(2,FEsolving_execELem(1))) then + startIP = FEsolving_execIP(1,FEsolving_execELem(1)) + endIP = startIP + else singleRun + startIP = 1 + endIP = theMesh%elem%nIPs + endif singleRun - NiterationCrystallite = 0 - cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) - NiterationCrystallite = NiterationCrystallite + 1 + NiterationCrystallite = 0 + cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) + NiterationCrystallite = NiterationCrystallite + 1 #ifdef DEBUG - if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0) & - write(6,'(a,i6)') '<< CRYST stress >> crystallite iteration ',NiterationCrystallite + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0) & + write(6,'(a,i6)') '<< CRYST stress >> crystallite iteration ',NiterationCrystallite #endif - !$OMP PARALLEL DO PRIVATE(formerSubStep) - elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1,homogenization_Ngrains(mesh_element(3,e)) + !$OMP PARALLEL DO PRIVATE(formerSubStep) + elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + do c = 1,homogenization_Ngrains(mesh_element(3,e)) !-------------------------------------------------------------------------------------------------- ! wind forward - if (crystallite_converged(c,i,e)) then - formerSubStep = crystallite_subStep(c,i,e) - crystallite_subFrac(c,i,e) = crystallite_subFrac(c,i,e) + crystallite_subStep(c,i,e) - crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & - stepIncreaseCryst * crystallite_subStep(c,i,e)) + if (crystallite_converged(c,i,e)) then + formerSubStep = crystallite_subStep(c,i,e) + crystallite_subFrac(c,i,e) = crystallite_subFrac(c,i,e) + crystallite_subStep(c,i,e) + crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & + stepIncreaseCryst * crystallite_subStep(c,i,e)) - crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? - if (crystallite_todo(c,i,e)) then - crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) - crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) - crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) - crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) - crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e) - crystallite_subS0 (1:3,1:3,c,i,e) = crystallite_S (1:3,1:3,c,i,e) - !if abbrevation, make c and p private in omp - plasticState( phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e)) & - = plasticState(phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) - do s = 1, phase_Nsources(phaseAt(c,i,e)) - sourceState( phaseAt(c,i,e))%p(s)%subState0(:,phasememberAt(c,i,e)) & - = sourceState(phaseAt(c,i,e))%p(s)%state( :,phasememberAt(c,i,e)) - enddo + crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? + if (crystallite_todo(c,i,e)) then + crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) + crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) + crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) + crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) + crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e) + crystallite_subS0 (1:3,1:3,c,i,e) = crystallite_S (1:3,1:3,c,i,e) + !if abbrevation, make c and p private in omp + plasticState( phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e)) & + = plasticState(phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) + do s = 1, phase_Nsources(phaseAt(c,i,e)) + sourceState( phaseAt(c,i,e))%p(s)%subState0(:,phasememberAt(c,i,e)) & + = sourceState(phaseAt(c,i,e))%p(s)%state( :,phasememberAt(c,i,e)) + enddo #ifdef DEBUG - if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0 & - .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) & - write(6,'(a,f12.8,a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST stress >> winding forward from ', & - crystallite_subFrac(c,i,e)-formerSubStep,' to current crystallite_subfrac ', & - crystallite_subFrac(c,i,e),' in crystallite_stress at el ip ipc ',e,i,c + if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0 & + .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) & + write(6,'(a,f12.8,a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST stress >> winding forward from ', & + crystallite_subFrac(c,i,e)-formerSubStep,' to current crystallite_subfrac ', & + crystallite_subFrac(c,i,e),' in crystallite_stress at el ip ipc ',e,i,c #endif - endif + endif !-------------------------------------------------------------------------------------------------- ! cut back (reduced time and restore) - else - crystallite_subStep(c,i,e) = subStepSizeCryst * crystallite_subStep(c,i,e) - crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) - crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp (1:3,1:3,c,i,e)) - crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) - crystallite_invFi(1:3,1:3,c,i,e) = math_inv33(crystallite_Fi (1:3,1:3,c,i,e)) - crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) - if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback - crystallite_Lp (1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) - crystallite_Li (1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e) - endif - plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) & - = plasticState(phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e)) - do s = 1, phase_Nsources(phaseAt(c,i,e)) - sourceState( phaseAt(c,i,e))%p(s)%state( :,phasememberAt(c,i,e)) & - = sourceState(phaseAt(c,i,e))%p(s)%subState0(:,phasememberAt(c,i,e)) - enddo + else + crystallite_subStep(c,i,e) = subStepSizeCryst * crystallite_subStep(c,i,e) + crystallite_Fp (1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) + crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp (1:3,1:3,c,i,e)) + crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) + crystallite_invFi(1:3,1:3,c,i,e) = math_inv33(crystallite_Fi (1:3,1:3,c,i,e)) + crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) + if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback + crystallite_Lp (1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) + crystallite_Li (1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e) + endif + plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) & + = plasticState(phaseAt(c,i,e))%subState0(:,phasememberAt(c,i,e)) + do s = 1, phase_Nsources(phaseAt(c,i,e)) + sourceState( phaseAt(c,i,e))%p(s)%state( :,phasememberAt(c,i,e)) & + = sourceState(phaseAt(c,i,e))%p(s)%subState0(:,phasememberAt(c,i,e)) + enddo - ! cant restore dotState here, since not yet calculated in first cutback after initialization - crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst ! still on track or already done (beyond repair) + ! cant restore dotState here, since not yet calculated in first cutback after initialization + crystallite_todo(c,i,e) = crystallite_subStep(c,i,e) > subStepMinCryst ! still on track or already done (beyond repair) #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & - .or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0)) then - if (crystallite_todo(c,i,e)) then - write(6,'(a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST stress >> cutback with new crystallite_subStep: ', & - crystallite_subStep(c,i,e),' at el ip ipc ',e,i,c - else - write(6,'(a,i8,1x,i2,1x,i3,/)') '<< CRYST stress >> reached minimum step size at el ip ipc ',e,i,c - endif - endif + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & + .or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0)) then + if (crystallite_todo(c,i,e)) then + write(6,'(a,f12.8,a,i8,1x,i2,1x,i3,/)') '<< CRYST stress >> cutback with new crystallite_subStep: ', & + crystallite_subStep(c,i,e),' at el ip ipc ',e,i,c + else + write(6,'(a,i8,1x,i2,1x,i3,/)') '<< CRYST stress >> reached minimum step size at el ip ipc ',e,i,c + endif + endif #endif - endif + endif !-------------------------------------------------------------------------------------------------- ! prepare for integration - if (crystallite_todo(c,i,e)) then - crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & - + crystallite_subStep(c,i,e) * (crystallite_partionedF (1:3,1:3,c,i,e) & - - crystallite_partionedF0(1:3,1:3,c,i,e)) - crystallite_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), & - crystallite_invFp(1:3,1:3,c,i,e)), & - crystallite_invFi(1:3,1:3,c,i,e)) - crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) - crystallite_converged(c,i,e) = .false. - endif + if (crystallite_todo(c,i,e)) then + crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & + + crystallite_subStep(c,i,e) * (crystallite_partionedF (1:3,1:3,c,i,e) & + - crystallite_partionedF0(1:3,1:3,c,i,e)) + crystallite_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), & + crystallite_invFp(1:3,1:3,c,i,e)), & + crystallite_invFi(1:3,1:3,c,i,e)) + crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) + crystallite_converged(c,i,e) = .false. + endif - enddo - enddo - enddo elementLooping3 - !$OMP END PARALLEL DO + enddo + enddo + enddo elementLooping3 + !$OMP END PARALLEL DO #ifdef DEBUG - if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0) then - write(6,'(/,a,f8.5,a,f8.5,/)') '<< CRYST stress >> ',minval(crystallite_subStep), & - ' ≤ subStep ≤ ',maxval(crystallite_subStep) - write(6,'(/,a,f8.5,a,f8.5,/)') '<< CRYST stress >> ',minval(crystallite_subFrac), & - ' ≤ subFrac ≤ ',maxval(crystallite_subFrac) - flush(6) - if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0) then - write(6,'(/,a,f8.5,1x,a,1x,f8.5,1x,a)') '<< CRYST stress >> subFrac + subStep = ',& - crystallite_subFrac(debug_g,debug_i,debug_e),'+',crystallite_subStep(debug_g,debug_i,debug_e),'@selective' - flush(6) - endif - endif + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0) then + write(6,'(/,a,f8.5,a,f8.5,/)') '<< CRYST stress >> ',minval(crystallite_subStep), & + ' ≤ subStep ≤ ',maxval(crystallite_subStep) + write(6,'(/,a,f8.5,a,f8.5,/)') '<< CRYST stress >> ',minval(crystallite_subFrac), & + ' ≤ subFrac ≤ ',maxval(crystallite_subFrac) + flush(6) + if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0) then + write(6,'(/,a,f8.5,1x,a,1x,f8.5,1x,a)') '<< CRYST stress >> subFrac + subStep = ',& + crystallite_subFrac(debug_g,debug_i,debug_e),'+',crystallite_subStep(debug_g,debug_i,debug_e),'@selective' + flush(6) + endif + endif #endif !-------------------------------------------------------------------------------------------------- ! integrate --- requires fully defined state array (basic + dependent state) - if (any(crystallite_todo)) call integrateState() ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation - where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & ! do not try non-converged but fully cutbacked any further - crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation + if (any(crystallite_todo)) call integrateState ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation + where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & ! do not try non-converged but fully cutbacked any further + crystallite_todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation - enddo cutbackLooping + enddo cutbackLooping ! return whether converged or not - crystallite_stress = .false. - elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) - enddo - enddo elementLooping5 + crystallite_stress = .false. + elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + crystallite_stress(i,e) = all(crystallite_converged(:,i,e)) + enddo + enddo elementLooping5 #ifdef DEBUG - elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1,homogenization_Ngrains(mesh_element(3,e)) - if (.not. crystallite_converged(c,i,e)) then - if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & - write(6,'(a,i8,1x,i2,1x,i3,/)') '<< CRYST stress >> no convergence at el ip ipc ', & - e,i,c - endif - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & - .or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0)) then - write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST stress >> solution at el ip ipc ',e,i,c - write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CRYST stress >> P / MPa', & - transpose(crystallite_P(1:3,1:3,c,i,e))*1.0e-6_pReal - write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Fp', & - transpose(crystallite_Fp(1:3,1:3,c,i,e)) - write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Fi', & - transpose(crystallite_Fi(1:3,1:3,c,i,e)) - write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST stress >> Lp', & - transpose(crystallite_Lp(1:3,1:3,c,i,e)) - write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST stress >> Li', & - transpose(crystallite_Li(1:3,1:3,c,i,e)) - flush(6) - endif - enddo - enddo - enddo elementLooping6 + elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + do c = 1,homogenization_Ngrains(mesh_element(3,e)) + if (.not. crystallite_converged(c,i,e)) then + if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & + write(6,'(a,i8,1x,i2,1x,i3,/)') '<< CRYST stress >> no convergence at el ip ipc ', & + e,i,c + endif + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & + .or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0)) then + write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST stress >> solution at el ip ipc ',e,i,c + write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CRYST stress >> P / MPa', & + transpose(crystallite_P(1:3,1:3,c,i,e))*1.0e-6_pReal + write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Fp', & + transpose(crystallite_Fp(1:3,1:3,c,i,e)) + write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST stress >> Fi', & + transpose(crystallite_Fi(1:3,1:3,c,i,e)) + write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST stress >> Lp', & + transpose(crystallite_Lp(1:3,1:3,c,i,e)) + write(6,'(a,/,3(12x,3(f14.9,1x)/),/)') '<< CRYST stress >> Li', & + transpose(crystallite_Li(1:3,1:3,c,i,e)) + flush(6) + endif + enddo + enddo + enddo elementLooping6 #endif end function crystallite_stress @@ -679,163 +679,163 @@ end function crystallite_stress !> @brief calculate tangent (dPdF) !-------------------------------------------------------------------------------------------------- subroutine crystallite_stressTangent() - use prec, only: & - tol_math_check, & - dNeq0 - use IO, only: & - IO_warning, & - IO_error - use math, only: & - math_inv33, & - math_identity2nd, & - math_3333to99, & - math_99to3333, & - math_I3, & - math_mul3333xx3333, & - math_mul33xx33, & - math_invert2, & - math_det33 - use mesh, only: & - mesh_element - use material, only: & - homogenization_Ngrains - use constitutive, only: & - constitutive_SandItsTangents, & - constitutive_LpAndItsTangents, & - constitutive_LiAndItsTangents + use prec, only: & + tol_math_check, & + dNeq0 + use IO, only: & + IO_warning, & + IO_error + use math, only: & + math_inv33, & + math_identity2nd, & + math_3333to99, & + math_99to3333, & + math_I3, & + math_mul3333xx3333, & + math_mul33xx33, & + math_invert2, & + math_det33 + use mesh, only: & + mesh_element + use material, only: & + homogenization_Ngrains + use constitutive, only: & + constitutive_SandItsTangents, & + constitutive_LpAndItsTangents, & + constitutive_LiAndItsTangents - implicit none - integer :: & - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e, & !< counter in element loop - o, & - p + implicit none + integer :: & + c, & !< counter in integration point component loop + i, & !< counter in integration point loop + e, & !< counter in element loop + o, & + p - real(pReal), dimension(3,3) :: temp_33_1, devNull,invSubFi0, temp_33_2, temp_33_3, temp_33_4 - real(pReal), dimension(3,3,3,3) :: dSdFe, & - dSdF, & - dSdFi, & - dLidS, & - dLidFi, & - dLpdS, & - dLpdFi, & - dFidS, & - dFpinvdF, & - rhs_3333, & - lhs_3333, & - temp_3333 - real(pReal), dimension(9,9):: temp_99 - logical :: error + real(pReal), dimension(3,3) :: temp_33_1, devNull,invSubFi0, temp_33_2, temp_33_3, temp_33_4 + real(pReal), dimension(3,3,3,3) :: dSdFe, & + dSdF, & + dSdFi, & + dLidS, & + dLidFi, & + dLpdS, & + dLpdFi, & + dFidS, & + dFpinvdF, & + rhs_3333, & + lhs_3333, & + temp_3333 + real(pReal), dimension(9,9):: temp_99 + logical :: error - !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,invSubFi0,o,p, & - !$OMP rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error) - elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1,homogenization_Ngrains(mesh_element(3,e)) + !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,invSubFi0,o,p, & + !$OMP rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error) + elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + do c = 1,homogenization_Ngrains(mesh_element(3,e)) - call constitutive_SandItsTangents(devNull,dSdFe,dSdFi, & - crystallite_Fe(1:3,1:3,c,i,e), & - crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent - call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & - crystallite_S (1:3,1:3,c,i,e), & - crystallite_Fi(1:3,1:3,c,i,e), & - c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration + call constitutive_SandItsTangents(devNull,dSdFe,dSdFi, & + crystallite_Fe(1:3,1:3,c,i,e), & + crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent + call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & + crystallite_S (1:3,1:3,c,i,e), & + crystallite_Fi(1:3,1:3,c,i,e), & + c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration - if (sum(abs(dLidS)) < tol_math_check) then - dFidS = 0.0_pReal - else - invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)) - lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal - do o=1,3; do p=1,3 - lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & - + crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) - lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & - + crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e) - rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & - - crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p)) - enddo;enddo - call math_invert2(temp_99,error,math_3333to99(lhs_3333)) - if (error) then - call IO_warning(warning_ID=600,el=e,ip=i,g=c, & - ext_msg='inversion error in analytic tangent calculation') - dFidS = 0.0_pReal - else - dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) - endif - dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS - endif + if (sum(abs(dLidS)) < tol_math_check) then + dFidS = 0.0_pReal + else + invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)) + lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal + do o=1,3; do p=1,3 + lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & + + crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) + lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & + + crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e) + rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & + - crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p)) + enddo;enddo + call math_invert2(temp_99,error,math_3333to99(lhs_3333)) + if (error) then + call IO_warning(warning_ID=600,el=e,ip=i,g=c, & + ext_msg='inversion error in analytic tangent calculation') + dFidS = 0.0_pReal + else + dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) + endif + dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS + endif - call constitutive_LpAndItsTangents(devNull,dLpdS,dLpdFi, & - crystallite_S (1:3,1:3,c,i,e), & - crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration - dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS + call constitutive_LpAndItsTangents(devNull,dLpdS,dLpdFi, & + crystallite_S (1:3,1:3,c,i,e), & + crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration + dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS !-------------------------------------------------------------------------------------------------- ! calculate dSdF - temp_33_1 = transpose(matmul(crystallite_invFp(1:3,1:3,c,i,e), & - crystallite_invFi(1:3,1:3,c,i,e))) - temp_33_2 = matmul( crystallite_subF (1:3,1:3,c,i,e), & - math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))) - temp_33_3 = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), & - crystallite_invFp (1:3,1:3,c,i,e)), & - math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))) + temp_33_1 = transpose(matmul(crystallite_invFp(1:3,1:3,c,i,e), & + crystallite_invFi(1:3,1:3,c,i,e))) + temp_33_2 = matmul( crystallite_subF (1:3,1:3,c,i,e), & + math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))) + temp_33_3 = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), & + crystallite_invFp (1:3,1:3,c,i,e)), & + math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))) - forall(p=1:3, o=1:3) - rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) - temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), & - crystallite_invFi(1:3,1:3,c,i,e)) & - + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) - end forall - lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) & - + math_mul3333xx3333(dSdFi,dFidS) + forall(p=1:3, o=1:3) + rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1) + temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), & + crystallite_invFi(1:3,1:3,c,i,e)) & + + matmul(temp_33_3,dLidS(1:3,1:3,p,o)) + end forall + lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) & + + math_mul3333xx3333(dSdFi,dFidS) - call math_invert2(temp_99,error,math_identity2nd(9)+math_3333to99(lhs_3333)) - if (error) then - call IO_warning(warning_ID=600,el=e,ip=i,g=c, & - ext_msg='inversion error in analytic tangent calculation') - dSdF = rhs_3333 - else - dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) - endif + call math_invert2(temp_99,error,math_identity2nd(9)+math_3333to99(lhs_3333)) + if (error) then + call IO_warning(warning_ID=600,el=e,ip=i,g=c, & + ext_msg='inversion error in analytic tangent calculation') + dSdF = rhs_3333 + else + dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) + endif !-------------------------------------------------------------------------------------------------- ! calculate dFpinvdF - temp_3333 = math_mul3333xx3333(dLpdS,dSdF) - forall(p=1:3, o=1:3) - dFpinvdF(1:3,1:3,p,o) & - = -crystallite_subdt(c,i,e) & - * matmul(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), & - matmul(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e))) - end forall + temp_3333 = math_mul3333xx3333(dLpdS,dSdF) + forall(p=1:3, o=1:3) + dFpinvdF(1:3,1:3,p,o) & + = -crystallite_subdt(c,i,e) & + * matmul(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), & + matmul(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e))) + end forall !-------------------------------------------------------------------------------------------------- ! assemble dPdF - temp_33_1 = matmul(crystallite_invFp(1:3,1:3,c,i,e), & - matmul(crystallite_S(1:3,1:3,c,i,e), & - transpose(crystallite_invFp(1:3,1:3,c,i,e)))) - temp_33_2 = matmul(crystallite_S(1:3,1:3,c,i,e), & - transpose(crystallite_invFp(1:3,1:3,c,i,e))) - temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e), & - crystallite_invFp(1:3,1:3,c,i,e)) - temp_33_4 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), & - crystallite_invFp(1:3,1:3,c,i,e)), & - crystallite_S(1:3,1:3,c,i,e)) + temp_33_1 = matmul(crystallite_invFp(1:3,1:3,c,i,e), & + matmul(crystallite_S(1:3,1:3,c,i,e), & + transpose(crystallite_invFp(1:3,1:3,c,i,e)))) + temp_33_2 = matmul(crystallite_S(1:3,1:3,c,i,e), & + transpose(crystallite_invFp(1:3,1:3,c,i,e))) + temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e), & + crystallite_invFp(1:3,1:3,c,i,e)) + temp_33_4 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), & + crystallite_invFp(1:3,1:3,c,i,e)), & + crystallite_S(1:3,1:3,c,i,e)) - crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal - do p=1, 3 - crystallite_dPdF(p,1:3,p,1:3,c,i,e) = transpose(temp_33_1) - enddo - forall(p=1:3, o=1:3) - crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & - matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33_2) + & - matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)),transpose(crystallite_invFp(1:3,1:3,c,i,e))) + & - matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o))) - end forall + crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal + do p=1, 3 + crystallite_dPdF(p,1:3,p,1:3,c,i,e) = transpose(temp_33_1) + enddo + forall(p=1:3, o=1:3) + crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & + matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33_2) + & + matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)),transpose(crystallite_invFp(1:3,1:3,c,i,e))) + & + matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o))) + end forall - enddo; enddo - enddo elementLooping - !$OMP END PARALLEL DO + enddo; enddo + enddo elementLooping + !$OMP END PARALLEL DO end subroutine crystallite_stressTangent @@ -844,44 +844,43 @@ end subroutine crystallite_stressTangent !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- subroutine crystallite_orientations - use math, only: & - math_rotationalPart33, & - math_RtoQ - use material, only: & - plasticState, & - material_phase, & - homogenization_Ngrains - use mesh, only: & - mesh_element - use lattice, only: & - lattice_qDisorientation - use plastic_nonlocal, only: & - plastic_nonlocal_updateCompatibility - - implicit none - integer & - c, & !< counter in integration point component loop - i, & !< counter in integration point loop - e !< counter in element loop - -!$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1,homogenization_Ngrains(mesh_element(3,e)) - call crystallite_orientation(c,i,e)%fromRotationMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) - enddo; enddo; enddo -!$OMP END PARALLEL DO + use math, only: & + math_rotationalPart33, & + math_RtoQ + use material, only: & + plasticState, & + material_phase, & + homogenization_Ngrains + use mesh, only: & + mesh_element + use lattice, only: & + lattice_qDisorientation + use plastic_nonlocal, only: & + plastic_nonlocal_updateCompatibility - ! --- we use crystallite_orientation from above, so need a separate loop - nonlocalPresent: if (any(plasticState%nonLocal)) then -!$OMP PARALLEL DO - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - if (plasticState(material_phase(1,i,e))%nonLocal) & ! if nonlocal model - call plastic_nonlocal_updateCompatibility(crystallite_orientation,i,e) - enddo; enddo -!$OMP END PARALLEL DO - endif nonlocalPresent + implicit none + integer & + c, & !< counter in integration point component loop + i, & !< counter in integration point loop + e !< counter in element loop + + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + do c = 1,homogenization_Ngrains(mesh_element(3,e)) + call crystallite_orientation(c,i,e)%fromRotationMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) + enddo; enddo; enddo + !$OMP END PARALLEL DO + + nonlocalPresent: if (any(plasticState%nonLocal)) then + !$OMP PARALLEL DO + do e = FEsolving_execElem(1),FEsolving_execElem(2) + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) + if (plasticState(material_phase(1,i,e))%nonLocal) & ! if nonlocal model + call plastic_nonlocal_updateCompatibility(crystallite_orientation,i,e) + enddo; enddo + !$OMP END PARALLEL DO + endif nonlocalPresent end subroutine crystallite_orientations @@ -890,24 +889,24 @@ end subroutine crystallite_orientations !> @brief Map 2nd order tensor to reference config !-------------------------------------------------------------------------------------------------- function crystallite_push33ToRef(ipc,ip,el, tensor33) - use math, only: & - math_inv33, & - math_EulerToR - use material, only: & - material_EulerAngles ! ToDo: Why stored? We also have crystallite_orientation0 - - implicit none - real(pReal), dimension(3,3) :: crystallite_push33ToRef - real(pReal), dimension(3,3), intent(in) :: tensor33 - real(pReal), dimension(3,3) :: T - integer, intent(in):: & - el, & ! element index - ip, & ! integration point index - ipc ! grain index - - T = matmul(math_EulerToR(material_EulerAngles(1:3,ipc,ip,el)), & - transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el)))) - crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) + use math, only: & + math_inv33, & + math_EulerToR + use material, only: & + material_EulerAngles ! ToDo: Why stored? We also have crystallite_orientation0 + + implicit none + real(pReal), dimension(3,3) :: crystallite_push33ToRef + real(pReal), dimension(3,3), intent(in) :: tensor33 + real(pReal), dimension(3,3) :: T + integer, intent(in):: & + el, & + ip, & + ipc + + T = matmul(math_EulerToR(material_EulerAngles(1:3,ipc,ip,el)), & + transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el)))) + crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) end function crystallite_push33ToRef @@ -1057,434 +1056,428 @@ end function crystallite_postResults !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -logical function integrateStress(& - ipc,& ! grain number - ip,& ! integration point number - el,& ! element number - timeFraction & - ) - use, intrinsic :: & - IEEE_arithmetic - use prec, only: tol_math_check, & - dEq0 - use numerics, only: nStress, & - aTol_crystalliteStress, & - rTol_crystalliteStress, & - iJacoLpresiduum, & - subStepSizeLp, & - subStepSizeLi +logical function integrateStress(ipc,ip,el,timeFraction) + use, intrinsic :: & + IEEE_arithmetic + use prec, only: tol_math_check, & + dEq0 + use numerics, only: nStress, & + aTol_crystalliteStress, & + rTol_crystalliteStress, & + iJacoLpresiduum, & + subStepSizeLp, & + subStepSizeLi #ifdef DEBUG - use debug, only: debug_level, & - debug_e, & - debug_i, & - debug_g, & - debug_crystallite, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective + use debug, only: debug_level, & + debug_e, & + debug_i, & + debug_g, & + debug_crystallite, & + debug_levelBasic, & + debug_levelExtensive, & + debug_levelSelective #endif - use constitutive, only: constitutive_LpAndItsTangents, & - constitutive_LiAndItsTangents, & - constitutive_SandItsTangents - use math, only: math_mul33xx33, & - math_mul3333xx3333, & - math_inv33, & - math_det33, & - math_I3, & - math_identity2nd, & - math_3333to99, & - math_33to9, & - math_9to33 - - implicit none - integer, intent(in):: el, & ! element index - ip, & ! integration point index - ipc ! grain index - real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep - - real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end 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 - Fi_new, & ! gradient of intermediate deformation stages - invFi_new, & - invFp_current, & ! inverse of Fp_current - invFi_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 - residuumLp, & ! current residuum of plastic velocity gradient - residuumLp_old, & ! last residuum of plastic velocity gradient - deltaLp, & ! direction of next guess - Liguess, & ! current guess for intermediate velocity gradient - Liguess_old, & ! known last good guess for intermediate velocity gradient - Li_constitutive, & ! intermediate velocity gradient resulting from constitutive law - residuumLi, & ! current residuum of intermediate velocity gradient - residuumLi_old, & ! last residuum of intermediate velocity gradient - deltaLi, & ! direction of next guess - S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration - A, & - B, & - Fe, & ! elastic deformation gradient - temp_33 - real(pReal), dimension(9) :: work ! needed for matrix inversion by LAPACK - integer, dimension(9) :: devNull ! needed for matrix inversion by LAPACK - real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme) - dRLp_dLp2, & ! working copy of dRdLp - dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme) - real(pReal), dimension(3,3,3,3):: dS_dFe, & ! partial derivative of 2nd Piola-Kirchhoff stress - dS_dFi, & - dFe_dLp, & ! partial derivative of elastic deformation gradient - dFe_dLi, & - dFi_dLi, & - dLp_dFi, & - dLi_dFi, & - dLp_dS, & - dLi_dS - real(pReal) detInvFi, & ! determinant of InvFi - steplengthLp, & - steplengthLi, & - dt, & ! time increment - aTolLp, & - aTolLi - integer NiterationStressLp, & ! number of stress integrations - NiterationStressLi, & ! number of inner stress integrations - ierr, & ! error indicator for LAPACK - o, & - p, & - jacoCounterLp, & - jacoCounterLi ! counters to check for Jacobian update - external :: & - dgesv - - !* be pessimistic - integrateStress = .false. + use constitutive, only: constitutive_LpAndItsTangents, & + constitutive_LiAndItsTangents, & + constitutive_SandItsTangents + use math, only: math_mul33xx33, & + math_mul3333xx3333, & + math_inv33, & + math_det33, & + math_I3, & + math_identity2nd, & + math_3333to99, & + math_33to9, & + math_9to33 + + implicit none + integer, intent(in):: el, & ! element index + ip, & ! integration point index + ipc ! grain index + real(pReal), optional, intent(in) :: timeFraction ! fraction of timestep + + real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end 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 + Fi_new, & ! gradient of intermediate deformation stages + invFi_new, & + invFp_current, & ! inverse of Fp_current + invFi_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 + residuumLp, & ! current residuum of plastic velocity gradient + residuumLp_old, & ! last residuum of plastic velocity gradient + deltaLp, & ! direction of next guess + Liguess, & ! current guess for intermediate velocity gradient + Liguess_old, & ! known last good guess for intermediate velocity gradient + Li_constitutive, & ! intermediate velocity gradient resulting from constitutive law + residuumLi, & ! current residuum of intermediate velocity gradient + residuumLi_old, & ! last residuum of intermediate velocity gradient + deltaLi, & ! direction of next guess + S, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration + A, & + B, & + Fe, & ! elastic deformation gradient + temp_33 + real(pReal), dimension(9) :: work ! needed for matrix inversion by LAPACK + integer, dimension(9) :: devNull ! needed for matrix inversion by LAPACK + real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme) + dRLp_dLp2, & ! working copy of dRdLp + dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme) + real(pReal), dimension(3,3,3,3):: dS_dFe, & ! partial derivative of 2nd Piola-Kirchhoff stress + dS_dFi, & + dFe_dLp, & ! partial derivative of elastic deformation gradient + dFe_dLi, & + dFi_dLi, & + dLp_dFi, & + dLi_dFi, & + dLp_dS, & + dLi_dS + real(pReal) detInvFi, & ! determinant of InvFi + steplengthLp, & + steplengthLi, & + dt, & ! time increment + aTolLp, & + aTolLi + integer NiterationStressLp, & ! number of stress integrations + NiterationStressLi, & ! number of inner stress integrations + ierr, & ! error indicator for LAPACK + o, & + p, & + jacoCounterLp, & + jacoCounterLi ! counters to check for Jacobian update + external :: & + dgesv + + !* be pessimistic + integrateStress = .false. #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) & - write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> at el ip ipc ',el,ip,ipc + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) & + write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> at el ip ipc ',el,ip,ipc #endif - if (present(timeFraction)) then - dt = crystallite_subdt(ipc,ip,el) * timeFraction - Fg_new = crystallite_subF0(1:3,1:3,ipc,ip,el) & - + (crystallite_subF(1:3,1:3,ipc,ip,el) - crystallite_subF0(1:3,1:3,ipc,ip,el)) * timeFraction - else - dt = crystallite_subdt(ipc,ip,el) - Fg_new = crystallite_subF(1:3,1:3,ipc,ip,el) - endif + if (present(timeFraction)) then + dt = crystallite_subdt(ipc,ip,el) * timeFraction + Fg_new = crystallite_subF0(1:3,1:3,ipc,ip,el) & + + (crystallite_subF(1:3,1:3,ipc,ip,el) - crystallite_subF0(1:3,1:3,ipc,ip,el)) * timeFraction + else + dt = crystallite_subdt(ipc,ip,el) + Fg_new = crystallite_subF(1:3,1:3,ipc,ip,el) + endif !* feed local variables - Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! ... and take it as first guess - Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! ... and take it as first guess - Liguess_old = Liguess - - invFp_current = math_inv33(crystallite_subFp0(1:3,1:3,ipc,ip,el)) - failedInversionFp: if (all(dEq0(invFp_current))) then + Lpguess = crystallite_Lp(1:3,1:3,ipc,ip,el) ! ... and take it as first guess + Liguess = crystallite_Li(1:3,1:3,ipc,ip,el) ! ... and take it as first guess + Liguess_old = Liguess + + invFp_current = math_inv33(crystallite_subFp0(1:3,1:3,ipc,ip,el)) + failedInversionFp: if (all(dEq0(invFp_current))) then #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & - write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on inversion of current Fp at el ip ipc ',& - el,ip,ipc - if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0) & - write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> current Fp ',transpose(crystallite_subFp0(1:3,1:3,ipc,ip,el)) + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & + write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on inversion of current Fp at el ip ipc ',& + el,ip,ipc + if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0) & + write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> current Fp ',transpose(crystallite_subFp0(1:3,1:3,ipc,ip,el)) #endif - return - endif failedInversionFp - A = matmul(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp + return + endif failedInversionFp + A = matmul(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp - invFi_current = math_inv33(crystallite_subFi0(1:3,1:3,ipc,ip,el)) - failedInversionFi: if (all(dEq0(invFi_current))) then + invFi_current = math_inv33(crystallite_subFi0(1:3,1:3,ipc,ip,el)) + failedInversionFi: if (all(dEq0(invFi_current))) then #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & - write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on inversion of current Fi at el ip ipc ',& - el,ip,ipc - if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0) & - write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> current Fi ', & - transpose(crystallite_subFi0(1:3,1:3,ipc,ip,el)) + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & + write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on inversion of current Fi at el ip ipc ',& + el,ip,ipc + if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0) & + write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> current Fi ', & + transpose(crystallite_subFi0(1:3,1:3,ipc,ip,el)) #endif - return - endif failedInversionFi - - !* start Li loop with normal step length - NiterationStressLi = 0 - jacoCounterLi = 0 - steplengthLi = 1.0_pReal - residuumLi_old = 0.0_pReal - - LiLoop: do - NiterationStressLi = NiterationStressLi + 1 - LiLoopLimit: if (NiterationStressLi > nStress) then + return + endif failedInversionFi + + !* start Li loop with normal step length + NiterationStressLi = 0 + jacoCounterLi = 0 + steplengthLi = 1.0_pReal + residuumLi_old = 0.0_pReal + + LiLoop: do + NiterationStressLi = NiterationStressLi + 1 + LiLoopLimit: if (NiterationStressLi > nStress) then #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & - write(6,'(a,i3,a,i8,1x,i2,1x,i3,/)') '<< CRYST integrateStress >> reached Li loop limit',nStress, & + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & + write(6,'(a,i3,a,i8,1x,i2,1x,i3,/)') '<< CRYST integrateStress >> reached Li loop limit',nStress, & ' at el ip ipc ', el,ip,ipc #endif - return - endif LiLoopLimit - - invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) - Fi_new = math_inv33(invFi_new) - detInvFi = math_det33(invFi_new) - - !* start Lp loop with normal step length - NiterationStressLp = 0 - jacoCounterLp = 0 - steplengthLp = 1.0_pReal - residuumLp_old = 0.0_pReal - Lpguess_old = Lpguess - - LpLoop: do - NiterationStressLp = NiterationStressLp + 1 - LpLoopLimit: if (NiterationStressLp > nStress) then + return + endif LiLoopLimit + + invFi_new = matmul(invFi_current,math_I3 - dt*Liguess) + Fi_new = math_inv33(invFi_new) + detInvFi = math_det33(invFi_new) + + !* start Lp loop with normal step length + NiterationStressLp = 0 + jacoCounterLp = 0 + steplengthLp = 1.0_pReal + residuumLp_old = 0.0_pReal + Lpguess_old = Lpguess + + LpLoop: do + NiterationStressLp = NiterationStressLp + 1 + LpLoopLimit: if (NiterationStressLp > nStress) then #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & - write(6,'(a,i3,a,i8,1x,i2,1x,i3,/)') '<< CRYST integrateStress >> reached Lp loop limit',nStress, & - ' at el ip ipc ', el,ip,ipc + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & + write(6,'(a,i3,a,i8,1x,i2,1x,i3,/)') '<< CRYST integrateStress >> reached Lp loop limit',nStress, & + ' at el ip ipc ', el,ip,ipc #endif - return - endif LpLoopLimit - - !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law - - B = math_I3 - dt*Lpguess - Fe = matmul(matmul(A,B), invFi_new) - call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration - - !* calculate plastic velocity gradient and its tangent from constitutive law - call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & - S, Fi_new, ipc, ip, el) - + return + endif LpLoopLimit + + !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law + + B = math_I3 - dt*Lpguess + Fe = matmul(matmul(A,B), invFi_new) + call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & + Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration + + !* calculate plastic velocity gradient and its tangent from constitutive law + call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & + S, Fi_new, ipc, ip, el) + #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then - write(6,'(a,i3,/)') '<< CRYST integrateStress >> iteration ', NiterationStressLp - write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lpguess', transpose(Lpguess) - write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lp_constitutive', transpose(Lp_constitutive) - write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Fi', transpose(Fi_new) - write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Fe', transpose(Fe) - write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> S', transpose(S) - endif + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,'(a,i3,/)') '<< CRYST integrateStress >> iteration ', NiterationStressLp + write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lpguess', transpose(Lpguess) + write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Lp_constitutive', transpose(Lp_constitutive) + write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Fi', transpose(Fi_new) + write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> Fe', transpose(Fe) + write(6,'(a,/,3(12x,3(e20.10,1x)/))') '<< CRYST integrateStress >> S', transpose(S) + endif #endif - !* update current residuum and check for convergence of loop - aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error - aTol_crystalliteStress) ! minimum lower cutoff - residuumLp = Lpguess - Lp_constitutive - - if (any(IEEE_is_NaN(residuumLp))) then + !* update current residuum and check for convergence of loop + aTolLp = max(rTol_crystalliteStress * max(norm2(Lpguess),norm2(Lp_constitutive)), & ! absolute tolerance from largest acceptable relative error + aTol_crystalliteStress) ! minimum lower cutoff + residuumLp = Lpguess - Lp_constitutive + + if (any(IEEE_is_NaN(residuumLp))) then #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & - write(6,'(a,i8,1x,i2,1x,i3,a,i3,a)') '<< CRYST integrateStress >> encountered NaN for Lp-residuum at el ip ipc ', & - el,ip,ipc, & - ' ; iteration ', NiterationStressLp,& - ' >> returning..!' + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & + write(6,'(a,i8,1x,i2,1x,i3,a,i3,a)') '<< CRYST integrateStress >> encountered NaN for Lp-residuum at el ip ipc ', & + el,ip,ipc, & + ' ; iteration ', NiterationStressLp,& + ' >> returning..!' #endif - return ! ...me = .false. to inform integrator about problem - elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance - exit LpLoop ! ...leave iteration loop - elseif ( NiterationStressLp == 1 & - .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLp_old = residuumLp ! ...remember old values and... - Lpguess_old = Lpguess - steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction - Lpguess = Lpguess_old + steplengthLp * deltaLp + return ! ...me = .false. to inform integrator about problem + elseif (norm2(residuumLp) < aTolLp) then ! converged if below absolute tolerance + exit LpLoop ! ...leave iteration loop + elseif ( NiterationStressLp == 1 & + .or. norm2(residuumLp) < norm2(residuumLp_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLp_old = residuumLp ! ...remember old values and... + Lpguess_old = Lpguess + steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLp = subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction + Lpguess = Lpguess_old + steplengthLp * deltaLp #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then - write(6,'(a,1x,f7.4)') '<< CRYST integrateStress >> linear search for Lpguess with step', steplengthLp - endif + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,'(a,1x,f7.4)') '<< CRYST integrateStress >> linear search for Lpguess with step', steplengthLp + endif #endif - cycle LpLoop - endif + cycle LpLoop + endif !* calculate Jacobian for correction term - if (mod(jacoCounterLp, iJacoLpresiduum) == 0) then - forall(o=1:3,p=1:3) dFe_dLp(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) - dFe_dLp = - dt * dFe_dLp - dRLp_dLp = math_identity2nd(9) & - - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) + if (mod(jacoCounterLp, iJacoLpresiduum) == 0) then + forall(o=1:3,p=1:3) dFe_dLp(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFe_dLp = - dt * dFe_dLp + dRLp_dLp = math_identity2nd(9) & + - math_3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp)) #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then - write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST integrateStress >> dLp_dS', math_3333to99(dLp_dS) - write(6,'(a,1x,e20.10)') '<< CRYST integrateStress >> dLp_dS norm', norm2(math_3333to99(dLp_dS)) - write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST integrateStress >> dRLp_dLp', dRLp_dLp-math_identity2nd(9) - write(6,'(a,1x,e20.10)') '<< CRYST integrateStress >> dRLp_dLp norm', norm2(dRLp_dLp-math_identity2nd(9)) - endif + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST integrateStress >> dLp_dS', math_3333to99(dLp_dS) + write(6,'(a,1x,e20.10)') '<< CRYST integrateStress >> dLp_dS norm', norm2(math_3333to99(dLp_dS)) + write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST integrateStress >> dRLp_dLp', dRLp_dLp-math_identity2nd(9) + write(6,'(a,1x,e20.10)') '<< CRYST integrateStress >> dRLp_dLp norm', norm2(dRLp_dLp-math_identity2nd(9)) + endif #endif - dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine - work = math_33to9(residuumLp) - call dgesv(9,1,dRLp_dLp2,9,devNull,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp - if (ierr /= 0) then + dRLp_dLp2 = dRLp_dLp ! will be overwritten in first call to LAPACK routine + work = math_33to9(residuumLp) + call dgesv(9,1,dRLp_dLp2,9,devNull,work,9,ierr) ! solve dRLp/dLp * delta Lp = -res for delta Lp + if (ierr /= 0) then #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then - write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on dR/dLp inversion at el ip ipc ', & - el,ip,ipc - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then - write(6,*) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dR_dLp',transpose(dRLp_dLp) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dFe_dLp',transpose(math_3333to99(dFe_dLp)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dS_dFe (cnst)',transpose(math_3333to99(dS_dFe)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dLp_dS (cnst)',transpose(math_3333to99(dLp_dS)) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> A',transpose(A) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> B',transpose(B) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Lp_constitutive',transpose(Lp_constitutive) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Lpguess',transpose(Lpguess) - endif - endif + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then + write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on dR/dLp inversion at el ip ipc ', & + el,ip,ipc + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,*) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dR_dLp',transpose(dRLp_dLp) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dFe_dLp',transpose(math_3333to99(dFe_dLp)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dS_dFe (cnst)',transpose(math_3333to99(dS_dFe)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dLp_dS (cnst)',transpose(math_3333to99(dLp_dS)) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> A',transpose(A) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> B',transpose(B) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Lp_constitutive',transpose(Lp_constitutive) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Lpguess',transpose(Lpguess) + endif + endif #endif - return - endif - deltaLp = - math_9to33(work) - endif - jacoCounterLp = jacoCounterLp + 1 + return + endif + deltaLp = - math_9to33(work) + endif + jacoCounterLp = jacoCounterLp + 1 - Lpguess = Lpguess + steplengthLp * deltaLp + Lpguess = Lpguess + steplengthLp * deltaLp - enddo LpLoop + enddo LpLoop !* calculate intermediate velocity gradient and its tangent from constitutive law - call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & - S, Fi_new, ipc, ip, el) + call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & + S, Fi_new, ipc, ip, el) #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Li_constitutive', transpose(Li_constitutive) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Liguess', transpose(Liguess) - endif + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Li_constitutive', transpose(Li_constitutive) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Liguess', transpose(Liguess) + endif #endif - !* update current residuum and check for convergence of loop - aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error - aTol_crystalliteStress) ! minimum lower cutoff - residuumLi = Liguess - Li_constitutive - if (any(IEEE_is_NaN(residuumLi))) then ! NaN in residuum... + !* update current residuum and check for convergence of loop + aTolLi = max(rTol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error + aTol_crystalliteStress) ! minimum lower cutoff + residuumLi = Liguess - Li_constitutive + if (any(IEEE_is_NaN(residuumLi))) then ! NaN in residuum... #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & - write(6,'(a,i8,1x,i2,1x,i3,a,i3,a)') '<< CRYST integrateStress >> encountered NaN for Li-residuum at el ip ipc ', & - el,ip,ipc, & - ' ; iteration ', NiterationStressLi,& - ' >> returning..!' + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) & + write(6,'(a,i8,1x,i2,1x,i3,a,i3,a)') '<< CRYST integrateStress >> encountered NaN for Li-residuum at el ip ipc ', & + el,ip,ipc, & + ' ; iteration ', NiterationStressLi,& + ' >> returning..!' #endif - return ! ...me = .false. to inform integrator about problem - elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance - exit LiLoop ! ...leave iteration loop - elseif ( NiterationStressLi == 1 & - .or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... - residuumLi_old = residuumLi ! ...remember old values and... - Liguess_old = Liguess - steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) - else ! not converged and residuum not improved... - steplengthLi = subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction - Liguess = Liguess_old + steplengthLi * deltaLi - cycle LiLoop - endif - - !* calculate Jacobian for correction term - if (mod(jacoCounterLi, iJacoLpresiduum) == 0) then - temp_33 = matmul(matmul(A,B),invFi_current) - forall(o=1:3,p=1:3) - dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) - dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current - end forall - forall(o=1:3,p=1:3) & - dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) - - dRLi_dLi = math_identity2nd(9) & - - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) + & - math_mul3333xx3333(dS_dFi, dFi_dLi))) & - - math_3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi)) - work = math_33to9(residuumLi) - call dgesv(9,1,dRLi_dLi,9,devNull,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li - if (ierr /= 0) then + return ! ...me = .false. to inform integrator about problem + elseif (norm2(residuumLi) < aTolLi) then ! converged if below absolute tolerance + exit LiLoop ! ...leave iteration loop + elseif ( NiterationStressLi == 1 & + .or. norm2(residuumLi) < norm2(residuumLi_old)) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumLi_old = residuumLi ! ...remember old values and... + Liguess_old = Liguess + steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthLi = subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction + Liguess = Liguess_old + steplengthLi * deltaLi + cycle LiLoop + endif + + !* calculate Jacobian for correction term + if (mod(jacoCounterLi, iJacoLpresiduum) == 0) then + temp_33 = matmul(matmul(A,B),invFi_current) + forall(o=1:3,p=1:3) + dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) + dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current + end forall + forall(o=1:3,p=1:3) & + dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) + + dRLi_dLi = math_identity2nd(9) & + - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) + & + math_mul3333xx3333(dS_dFi, dFi_dLi))) & + - math_3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi)) + work = math_33to9(residuumLi) + call dgesv(9,1,dRLi_dLi,9,devNull,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li + if (ierr /= 0) then #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then - write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on dR/dLi inversion at el ip ipc ', & - el,ip,ipc - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then - write(6,*) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dR_dLi',transpose(dRLi_dLi) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dFe_dLi',transpose(math_3333to99(dFe_dLi)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dS_dFi (cnst)',transpose(math_3333to99(dS_dFi)) - write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dLi_dS (cnst)',transpose(math_3333to99(dLi_dS)) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Li_constitutive',transpose(Li_constitutive) - write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Liguess',transpose(Liguess) - endif - endif + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then + write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on dR/dLi inversion at el ip ipc ', & + el,ip,ipc + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g)& + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,*) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dR_dLi',transpose(dRLi_dLi) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dFe_dLi',transpose(math_3333to99(dFe_dLi)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dS_dFi (cnst)',transpose(math_3333to99(dS_dFi)) + write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST integrateStress >> dLi_dS (cnst)',transpose(math_3333to99(dLi_dS)) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Li_constitutive',transpose(Li_constitutive) + write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST integrateStress >> Liguess',transpose(Liguess) + endif + endif #endif - return - endif - - deltaLi = - math_9to33(work) - endif - jacoCounterLi = jacoCounterLi + 1 - - Liguess = Liguess + steplengthLi * deltaLi - enddo LiLoop - - !* calculate new plastic and elastic deformation gradient - invFp_new = matmul(invFp_current,B) - invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize - Fp_new = math_inv33(invFp_new) - failedInversionInvFp: if (all(dEq0(Fp_new))) then + return + endif + + deltaLi = - math_9to33(work) + endif + jacoCounterLi = jacoCounterLi + 1 + + Liguess = Liguess + steplengthLi * deltaLi + enddo LiLoop + + !* calculate new plastic and elastic deformation gradient + invFp_new = matmul(invFp_current,B) + invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize + Fp_new = math_inv33(invFp_new) + failedInversionInvFp: if (all(dEq0(Fp_new))) then #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then - write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on invFp_new inversion at el ip ipc ', & - el,ip,ipc - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) & - write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> invFp_new',transpose(invFp_new) - endif + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0) then + write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST integrateStress >> failed on invFp_new inversion at el ip ipc ', & + el,ip,ipc + if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) & + write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> invFp_new',transpose(invFp_new) + endif #endif - return - endif failedInversionInvFp - Fe_new = matmul(matmul(Fg_new,invFp_new),invFi_new) + return + endif failedInversionInvFp + Fe_new = matmul(matmul(Fg_new,invFp_new),invFi_new) !-------------------------------------------------------------------------------------------------- ! stress integration was successful - integrateStress = .true. - crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(Fg_new,invFp_new), & - matmul(S,transpose(invFp_new))) - crystallite_S (1:3,1:3,ipc,ip,el) = S - crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess - crystallite_Li (1:3,1:3,ipc,ip,el) = Liguess - crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new - crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new - crystallite_Fe (1:3,1:3,ipc,ip,el) = Fe_new - crystallite_invFp(1:3,1:3,ipc,ip,el) = invFp_new - crystallite_invFi(1:3,1:3,ipc,ip,el) = invFi_new + integrateStress = .true. + crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(Fg_new,invFp_new),matmul(S,transpose(invFp_new))) + crystallite_S (1:3,1:3,ipc,ip,el) = S + crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess + crystallite_Li (1:3,1:3,ipc,ip,el) = Liguess + crystallite_Fp (1:3,1:3,ipc,ip,el) = Fp_new + crystallite_Fi (1:3,1:3,ipc,ip,el) = Fi_new + crystallite_Fe (1:3,1:3,ipc,ip,el) = Fe_new + crystallite_invFp(1:3,1:3,ipc,ip,el) = invFp_new + crystallite_invFi(1:3,1:3,ipc,ip,el) = invFi_new #ifdef DEBUG - if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0 & - .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & - .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then - write(6,'(a,/)') '<< CRYST integrateStress >> successful integration' - write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> P / MPa', & - transpose(crystallite_P(1:3,1:3,ipc,ip,el))*1.0e-6_pReal - write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Cauchy / MPa', & - matmul(crystallite_P(1:3,1:3,ipc,ip,el), transpose(Fg_new)) * 1.0e-6_pReal / math_det33(Fg_new) - write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fe Lp Fe^-1', & - transpose(matmul(Fe_new, matmul(crystallite_Lp(1:3,1:3,ipc,ip,el), math_inv33(Fe_new)))) - write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fp',transpose(crystallite_Fp(1:3,1:3,ipc,ip,el)) - write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fi',transpose(crystallite_Fi(1:3,1:3,ipc,ip,el)) - endif + if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0 & + .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & + .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0)) then + write(6,'(a,/)') '<< CRYST integrateStress >> successful integration' + write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> P / MPa', & + transpose(crystallite_P(1:3,1:3,ipc,ip,el))*1.0e-6_pReal + write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Cauchy / MPa', & + matmul(crystallite_P(1:3,1:3,ipc,ip,el), transpose(Fg_new)) * 1.0e-6_pReal / math_det33(Fg_new) + write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fe Lp Fe^-1', & + transpose(matmul(Fe_new, matmul(crystallite_Lp(1:3,1:3,ipc,ip,el), math_inv33(Fe_new)))) + write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fp',transpose(crystallite_Fp(1:3,1:3,ipc,ip,el)) + write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fi',transpose(crystallite_Fi(1:3,1:3,ipc,ip,el)) + endif #endif end function integrateStress