diff --git a/PRIVATE b/PRIVATE index df3b567c9..02a391e38 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit df3b567c90dbcb291a5c2563d3984fce346cb55b +Subproject commit 02a391e38d1a44a8cfd2f169687a9d9f314a8146 diff --git a/VERSION b/VERSION index 429ed3632..26c337c96 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-558-gdfba3c9b +v2.0.2-591-ga00d15b8 diff --git a/examples/ConfigFiles/Crystallite_All.config b/examples/ConfigFiles/Crystallite_All.config index 761380fcd..d46c3e0e6 100644 --- a/examples/ConfigFiles/Crystallite_All.config +++ b/examples/ConfigFiles/Crystallite_All.config @@ -5,15 +5,10 @@ (output) orientation # quaternion (output) eulerangles # orientation as Bunge triple in degree (output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) in crystal reference coordinates -(output) grainrotationx # deviation from initial orientation as angle in degrees around sample reference x axis -(output) grainrotationy # deviation from initial orientation as angle in degrees around sample reference y axis -(output) grainrotationz # deviation from initial orientation as angle in degrees around sample reference z axis -(output) f # deformation gradient tensor; synonyms: "defgrad" +(output) f # deformation gradient tensor (output) fe # elastic deformation gradient tensor (output) fp # plastic deformation gradient tensor -(output) e # total strain as Green-Lagrange tensor -(output) ee # elastic strain as Green-Lagrange tensor -(output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola" -(output) s # second Piola-Kichhoff stress tensor; synonyms: "tstar", "secondpiola", "2ndpiola" +(output) p # first Piola-Kichhoff stress tensor +(output) s # second Piola-Kichhoff stress tensor (output) lp # plastic velocity gradient tensor (output) elasmatrix # elastic stiffness matrix diff --git a/examples/SpectralMethod/Polycrystal/material.config b/examples/SpectralMethod/Polycrystal/material.config index 5073f165e..93a5a6710 100644 --- a/examples/SpectralMethod/Polycrystal/material.config +++ b/examples/SpectralMethod/Polycrystal/material.config @@ -18,8 +18,6 @@ mech none (output) f # deformation gradient tensor; synonyms: "defgrad" (output) fe # elastic deformation gradient tensor (output) fp # plastic deformation gradient tensor -(output) e # total strain as Green-Lagrange tensor -(output) ee # elastic strain as Green-Lagrange tensor (output) p # first Piola-Kichhoff stress tensor; synonyms: "firstpiola", "1stpiola" (output) lp # plastic velocity gradient tensor @@ -56,7 +54,6 @@ h0_slipslip 75e6 interaction_slipslip 1 1 1.4 1.4 1.4 1.4 atol_resistance 1 - #-------------------# #-------------------# diff --git a/installation/patch/python2to3.sh b/installation/patch/python2to3.sh index 1d86b0ce7..07eed0e1e 100755 --- a/installation/patch/python2to3.sh +++ b/installation/patch/python2to3.sh @@ -1,8 +1,8 @@ #! /usr/bin/env bash if [ $1x != 3to2x ]; then echo 'python2.7 to python' - find . -name '*.py' | xargs sed -i 's/usr\/bin\/env python2.7/usr\/bin\/env python/g' + find . -name '*.py' -type f | xargs sed -i 's/usr\/bin\/env python2.7/usr\/bin\/env python/g' else echo 'python to python2.7' - find . -name '*.py' | xargs sed -i 's/usr\/bin\/env python/usr\/bin\/env python2.7/g' + find . -name '*.py' -type f | xargs sed -i 's/usr\/bin\/env python/usr\/bin\/env python2.7/g' fi diff --git a/lib/damask/asciitable.py b/lib/damask/asciitable.py index c47dbd8a7..839adb944 100644 --- a/lib/damask/asciitable.py +++ b/lib/damask/asciitable.py @@ -1,14 +1,14 @@ # -*- coding: UTF-8 no BOM -*- -import os,sys +import os, sys import numpy as np # ------------------------------------------------------------------ # python 3 has no unicode object, this ensures that the code works on Python 2&3 try: - test=isinstance('test', unicode) -except(NameError): - unicode=str + test = isinstance('test', unicode) +except NameError: + unicode = str # ------------------------------------------------------------------ class ASCIItable(): @@ -63,16 +63,17 @@ class ASCIItable(): x): try: return float(x) - except: + except ValueError: return 0.0 # ------------------------------------------------------------------ def _removeCRLF(self, - string): + string): + """Delete any carriage return and line feed from string""" try: return string.replace('\n','').replace('\r','') - except: - return string + except AttributeError: + return str(string) # ------------------------------------------------------------------ @@ -93,22 +94,22 @@ class ASCIItable(): # ------------------------------------------------------------------ def input_close(self): - try: +# try: if self.__IO__['in'] != sys.stdin: self.__IO__['in'].close() - except: - pass +# except: +# pass # ------------------------------------------------------------------ def output_write(self, what): """Aggregate a single row (string) or list of (possibly containing further lists of) rows into output""" - if not isinstance(what, (str, unicode)): + if isinstance(what, (str, unicode)): + self.__IO__['output'] += [what] + else: try: for item in what: self.output_write(item) - except: + except TypeError: self.__IO__['output'] += [str(what)] - else: - self.__IO__['output'] += [what] return self.__IO__['buffered'] or self.output_flush() @@ -129,10 +130,10 @@ class ASCIItable(): # ------------------------------------------------------------------ def output_close(self, dismiss = False): - try: - if self.__IO__['out'] != sys.stdout: self.__IO__['out'].close() - except: - pass +# try: + if self.__IO__['out'] != sys.stdout: self.__IO__['out'].close() +# except: +# pass if dismiss and os.path.isfile(self.__IO__['out'].name): os.remove(self.__IO__['out'].name) elif self.__IO__['inPlace']: @@ -150,7 +151,7 @@ class ASCIItable(): try: self.__IO__['in'].seek(0) - except: + except IOError: pass firstline = self.__IO__['in'].readline().strip() @@ -170,7 +171,7 @@ class ASCIItable(): else: # other table format try: self.__IO__['in'].seek(0) # try to rewind - except: + except IOError: self.__IO__['readBuffer'] = [firstline] # or at least save data in buffer while self.data_read(advance = False, respectLabels = False): @@ -197,7 +198,9 @@ class ASCIItable(): """Write current header information (info + labels)""" head = ['{}\theader'.format(len(self.info)+self.__IO__['labeled'])] if header else [] head.append(self.info) - if self.__IO__['labeled']: head.append('\t'.join(map(self._quote,self.tags))) + if self.__IO__['labeled']: + head.append('\t'.join(map(self._quote,self.tags))) + if len(self.tags) == 0: raise ValueError('no labels present.') return self.output_write(head) @@ -257,13 +260,13 @@ class ASCIItable(): what, reset = False): """Add item or list to existing set of labels (and switch on labeling)""" - if not isinstance(what, (str, unicode)): + if isinstance(what, (str, unicode)): + self.tags += [self._removeCRLF(what)] + else: try: for item in what: self.labels_append(item) - except: + except TypeError: self.tags += [self._removeCRLF(str(what))] - else: - self.tags += [self._removeCRLF(what)] self.__IO__['labeled'] = True # switch on processing (in particular writing) of tags if reset: self.__IO__['tags'] = list(self.tags) # subsequent data_read uses current tags as data size @@ -410,13 +413,13 @@ class ASCIItable(): def info_append(self, what): """Add item or list to existing set of infos""" - if not isinstance(what, (str, unicode)): + if isinstance(what, (str, unicode)): + self.info += [self._removeCRLF(what)] + else: try: for item in what: self.info_append(item) - except: + except TypeError: self.info += [self._removeCRLF(str(what))] - else: - self.info += [self._removeCRLF(what)] # ------------------------------------------------------------------ def info_clear(self): @@ -468,10 +471,8 @@ class ASCIItable(): """Read whole data of all (given) labels as numpy array""" from collections import Iterable - try: - self.data_rewind() # try to wind back to start of data - except: - pass # assume/hope we are at data start already... + try: self.data_rewind() # try to wind back to start of data + except: pass # assume/hope we are at data start already... if labels is None or labels == []: use = None # use all columns (and keep labels intact) @@ -530,13 +531,13 @@ class ASCIItable(): # ------------------------------------------------------------------ def data_append(self, what): - if not isinstance(what, (str, unicode)): + if isinstance(what, (str, unicode)): + self.data += [what] + else: try: for item in what: self.data_append(item) - except: + except TypeError: self.data += [str(what)] - else: - self.data += [what] # ------------------------------------------------------------------ def data_set(self, @@ -581,7 +582,7 @@ class ASCIItable(): if len(items) > 2: if items[1].lower() == 'of': items = np.ones(datatype(items[0]))*datatype(items[2]) - elif items[1].lower() == 'to': + elif items[1].lower() == 'to': items = np.linspace(datatype(items[0]),datatype(items[2]), abs(datatype(items[2])-datatype(items[0]))+1,dtype=int) else: items = list(map(datatype,items)) diff --git a/processing/pre/geom_fromOsteonGeometry.py b/processing/pre/geom_fromOsteonGeometry.py index 807e5200e..dc481ef97 100755 --- a/processing/pre/geom_fromOsteonGeometry.py +++ b/processing/pre/geom_fromOsteonGeometry.py @@ -65,7 +65,7 @@ if np.any(np.array(options.size) <= 0.0): if filename == []: filename = [None] table = damask.ASCIItable(outname = filename[0], - buffered = False) + buffered = False, labeled=False) damask.util.report(scriptName,filename[0]) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index d34eeb7fe..c44f90b41 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -440,7 +440,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e math_Mandel33to6, & math_Plain99to3333 use material, only: & + phasememberAt, & phase_plasticity, & + phase_plasticityInstance, & material_phase, & material_homog, & temperature, & @@ -490,7 +492,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e ho, & !< homogenization tme !< thermal member position integer(pInt) :: & - i, j + i, j, instance, of ho = material_homog(ip,el) tme = thermalMapping(ho)%p(ip,el) @@ -509,8 +511,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) - dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of) case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el) @@ -817,6 +820,7 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac debug_constitutive, & debug_levelBasic use math, only: & + math_mul33x33, & math_Mandel6to33, & math_Mandel33to6, & math_mul33x33 @@ -824,6 +828,8 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac mesh_NcpElems, & mesh_maxNips use material, only: & + phasememberAt, & + phase_plasticityInstance, & phase_plasticity, & phase_source, & phase_Nsources, & @@ -881,38 +887,41 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac real(pReal), intent(in), dimension(6) :: & S6 !< 2nd Piola Kirchhoff stress (vector notation) real(pReal), dimension(3,3) :: & - Mstar + Mp integer(pInt) :: & ho, & !< homogenization tme, & !< thermal member position - s !< counter in source loop + s, & !< counter in source loop + instance, of ho = material_homog( ip,el) tme = thermalMapping(ho)%p(ip,el) - Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) + Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_dotState (math_Mandel33to6(Mstar),ipc,ip,el) + call plastic_isotropic_dotState (math_Mandel33to6(Mp),ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_dotState(math_Mandel33to6(Mstar),ipc,ip,el) + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + call plastic_phenopowerlaw_dotState(Mp,instance,of) case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_dotState(math_Mandel33to6(Mstar),ipc,ip,el) + call plastic_kinehardening_dotState(math_Mandel33to6(Mp),ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dotState (Mstar,temperature(ho)%p(tme), & + call plastic_dislotwin_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), & ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), & + call plastic_disloucla_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), & ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dotState (math_Mandel33to6(Mstar),FeArray,FpArray,temperature(ho)%p(tme), & + call plastic_nonlocal_dotState (math_Mandel33to6(Mp),FeArray,FpArray,temperature(ho)%p(tme), & subdt,subfracArray,ip,el) end select plasticityType @@ -1026,13 +1035,18 @@ end subroutine constitutive_collectDeltaState !-------------------------------------------------------------------------------------------------- !> @brief returns array of constitutive results !-------------------------------------------------------------------------------------------------- -function constitutive_postResults(S6, FeArray, ipc, ip, el) +function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el) use prec, only: & pReal + use math, only: & + math_Mandel6to33, & + math_mul33x33 use mesh, only: & mesh_NcpElems, & mesh_maxNips use material, only: & + phasememberAt, & + phase_plasticityInstance, & plasticState, & sourceState, & phase_plasticity, & @@ -1083,19 +1097,25 @@ function constitutive_postResults(S6, FeArray, ipc, ip, el) real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults + & sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: & constitutive_postResults + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & FeArray !< elastic deformation gradient real(pReal), intent(in), dimension(6) :: & S6 !< 2nd Piola Kirchhoff stress (vector notation) + real(pReal), dimension(3,3) :: & + Mp !< Mandel stress integer(pInt) :: & startPos, endPos integer(pInt) :: & ho, & !< homogenization tme, & !< thermal member position - s !< counter in source loop + s, of, instance !< counter in source loop constitutive_postResults = 0.0_pReal + Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6)) + ho = material_homog( ip,el) tme = thermalMapping(ho)%p(ip,el) @@ -1104,10 +1124,13 @@ function constitutive_postResults(S6, FeArray, ipc, ip, el) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) case (PLASTICITY_ISOTROPIC_ID) plasticityType - constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(S6,ipc,ip,el) - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType constitutive_postResults(startPos:endPos) = & - plastic_phenopowerlaw_postResults(S6,ipc,ip,el) + plastic_isotropic_postResults(S6,ipc,ip,el) + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + of = phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phase(ipc,ip,el)) + constitutive_postResults(startPos:endPos) = & + plastic_phenopowerlaw_postResults(Mp,instance,of) case (PLASTICITY_KINEHARDENING_ID) plasticityType constitutive_postResults(startPos:endPos) = & plastic_kinehardening_postResults(S6,ipc,ip,el) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 94ba805e8..4082749b2 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -38,6 +38,9 @@ module crystallite crystallite_orientation, & !< orientation as quaternion crystallite_orientation0, & !< initial orientation as quaternion crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame + real(pReal), dimension(:,:,:,:,:), allocatable, public, 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_Fp, & !< current plastic def grad (end of converged time step) crystallite_Fp0, & !< plastic def grad at start of FE inc @@ -50,14 +53,11 @@ module crystallite 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_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 - crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - crystallite_P !< 1st Piola-Kirchhoff stress per grain + crystallite_partionedLi0 !< intermediate velocity grad at start of homog inc real(pReal), dimension(:,:,:,:,:), allocatable, private :: & - crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_subFp0,& !< plastic def grad at start of crystallite inc crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step) @@ -71,8 +71,6 @@ module crystallite crystallite_dPdF, & !< current individual dPdF per grain (end of converged time step) crystallite_dPdF0, & !< individual dPdF per grain at start of FE inc crystallite_partioneddPdF0 !< individual dPdF per grain at start of homog inc - real(pReal), dimension(:,:,:,:,:,:,:), allocatable, private :: & - crystallite_fallbackdPdF !< dPdF fallback for non-converged grains (elastic prediction) logical, dimension(:,:,:), allocatable, public :: & crystallite_requested !< flag to request crystallite calculation logical, dimension(:,:,:), allocatable, public, protected :: & @@ -92,9 +90,6 @@ module crystallite phase_ID, & texture_ID, & volume_ID, & - grainrotationx_ID, & - grainrotationy_ID, & - grainrotationz_ID, & orientation_ID, & grainrotation_ID, & eulerangles_ID, & @@ -104,8 +99,6 @@ module crystallite fi_ID, & lp_ID, & li_ID, & - e_ID, & - ee_ID, & p_ID, & s_ID, & elasmatrix_ID, & @@ -124,13 +117,13 @@ module crystallite crystallite_postResults private :: & integrateState, & - crystallite_integrateStateFPI, & - crystallite_integrateStateEuler, & - crystallite_integrateStateAdaptiveEuler, & - crystallite_integrateStateRK4, & - crystallite_integrateStateRKCK45, & - crystallite_integrateStress, & - crystallite_stateJump + integrateStateFPI, & + integrateStateEuler, & + integrateStateAdaptiveEuler, & + integrateStateRK4, & + integrateStateRKCK45, & + integrateStress, & + stateJump contains @@ -144,12 +137,14 @@ subroutine crystallite_init compiler_version, & compiler_options #endif +#ifdef DEBUG use debug, only: & debug_info, & debug_reset, & debug_level, & debug_crystallite, & debug_levelBasic +#endif use numerics, only: & numerics_integrator, & worldrank, & @@ -177,8 +172,7 @@ subroutine crystallite_init use config, only: & config_deallocate, & config_crystallite, & - crystallite_name, & - material_Nphase + crystallite_name use constitutive, only: & constitutive_initialFi, & constitutive_microstructure ! derived (shortcut) quantities of given state @@ -192,7 +186,6 @@ subroutine crystallite_init e, & !< counter in element loop o = 0_pInt, & !< counter in output loop r, & - ph, & !< counter in crystallite loop cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax, & !< maximum number of elements @@ -201,8 +194,6 @@ subroutine crystallite_init mySize character(len=65536), dimension(:), allocatable :: str - character(len=65536) :: & - tag = '' write(6,'(/,a)') ' <<<+- crystallite init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -235,7 +226,6 @@ subroutine crystallite_init 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_subFe0(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) @@ -247,7 +237,6 @@ subroutine crystallite_init allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_dPdF0(3,3,3,3,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_partioneddPdF0(3,3,3,3,cMax,iMax,eMax),source=0.0_pReal) - allocate(crystallite_fallbackdPdF(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) @@ -276,15 +265,15 @@ subroutine crystallite_init select case(numerics_integrator(1)) case(1_pInt) - integrateState => crystallite_integrateStateFPI + integrateState => integrateStateFPI case(2_pInt) - integrateState => crystallite_integrateStateEuler + integrateState => integrateStateEuler case(3_pInt) - integrateState => crystallite_integrateStateAdaptiveEuler + integrateState => integrateStateAdaptiveEuler case(4_pInt) - integrateState => crystallite_integrateStateRK4 + integrateState => integrateStateRK4 case(5_pInt) - integrateState => crystallite_integrateStateRKCK45 + integrateState => integrateStateRKCK45 end select @@ -306,12 +295,6 @@ subroutine crystallite_init crystallite_outputID(o,c) = texture_ID case ('volume') outputName crystallite_outputID(o,c) = volume_ID - case ('grainrotationx') outputName - crystallite_outputID(o,c) = grainrotationx_ID - case ('grainrotationy') outputName - crystallite_outputID(o,c) = grainrotationy_ID - case ('grainrotationz') outputName - crystallite_outputID(o,c) = grainrotationx_ID case ('orientation') outputName crystallite_outputID(o,c) = orientation_ID case ('grainrotation') outputName @@ -330,10 +313,6 @@ subroutine crystallite_init crystallite_outputID(o,c) = lp_ID case ('li') outputName crystallite_outputID(o,c) = li_ID - case ('e') outputName - crystallite_outputID(o,c) = e_ID - case ('ee') outputName - crystallite_outputID(o,c) = ee_ID case ('p','firstpiola','1stpiola') outputName crystallite_outputID(o,c) = p_ID case ('s','tstar','secondpiola','2ndpiola') outputName @@ -345,7 +324,7 @@ subroutine crystallite_init case ('neighboringelement') outputName crystallite_outputID(o,c) = neighboringelement_ID case default outputName - call IO_error(105_pInt,ext_msg=tag//' (Crystallite)') + call IO_error(105_pInt,ext_msg=trim(str(o))//' (Crystallite)') end select outputName enddo enddo @@ -354,13 +333,13 @@ subroutine crystallite_init do r = 1_pInt,size(config_crystallite) do o = 1_pInt,crystallite_Noutput(r) select case(crystallite_outputID(o,r)) - case(phase_ID,texture_ID,volume_ID,grainrotationx_ID,grainrotationy_ID,grainrotationz_ID) + case(phase_ID,texture_ID,volume_ID) mySize = 1_pInt case(orientation_ID,grainrotation_ID) mySize = 4_pInt case(eulerangles_ID) mySize = 3_pInt - case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,e_ID,ee_ID,p_ID,s_ID) + case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,p_ID,s_ID) mySize = 9_pInt case(elasmatrix_ID) mySize = 36_pInt @@ -426,11 +405,10 @@ subroutine crystallite_init call crystallite_orientations() crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations - !$OMP PARALLEL DO PRIVATE(myNcomponents) + !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1_pInt,myNcomponents + do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) call constitutive_microstructure(crystallite_orientation, & ! pass orientation to constitutive module crystallite_Fe(1:3,1:3,c,i,e), & crystallite_Fp(1:3,1:3,c,i,e), & @@ -440,77 +418,21 @@ subroutine crystallite_init enddo !$OMP END PARALLEL DO -!-------------------------------------------------------------------------------------------------- -! propagate dependent states to materialpoint and boundary value problem level -! do ph = 1_pInt,material_Nphase -! plasticState(ph)%partionedState0(plasticState(ph)%offsetDeltaState+plasticState(ph)%sizeDeltaState: & -! plasticState(ph)%sizeState,:) & -! = plasticState(ph)%state(plasticState(ph)%offsetDeltaState+plasticState(ph)%sizeDeltaState: & -! plasticState(ph)%sizeState,:) -! plasticState(ph)%state0 (plasticState(ph)%offsetDeltaState+plasticState(ph)%sizeDeltaState: & -! plasticState(ph)%sizeState,:) & -! = plasticState(ph)%state(plasticState(ph)%offsetDeltaState+plasticState(ph)%sizeDeltaState: & -! plasticState(ph)%sizeState,:) -! enddo - call crystallite_stressAndItsTangent(.true.) ! request elastic answers - crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback -!-------------------------------------------------------------------------------------------------- -! debug output +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fe: ', shape(crystallite_Fe) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp: ', shape(crystallite_Fp) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fi: ', shape(crystallite_Fi) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp: ', shape(crystallite_Lp) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Li: ', shape(crystallite_Li) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_F0: ', shape(crystallite_F0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fp0: ', shape(crystallite_Fp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Fi0: ', shape(crystallite_Fi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Lp0: ', shape(crystallite_Lp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Li0: ', shape(crystallite_Li0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF: ', shape(crystallite_partionedF) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFi0: ', shape(crystallite_partionedFi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLi0: ', shape(crystallite_partionedLi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF: ', shape(crystallite_subF) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF0: ', shape(crystallite_subF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFe0: ', shape(crystallite_subFe0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFi0: ', shape(crystallite_subFi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subLi0: ', shape(crystallite_subLi0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_P: ', shape(crystallite_P) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_Tstar0_v: ', shape(crystallite_Tstar0_v) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedTstar0_v: ', shape(crystallite_partionedTstar0_v) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_dPdF0: ', shape(crystallite_dPdF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_partioneddPdF0: ', shape(crystallite_partioneddPdF0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation: ', shape(crystallite_orientation) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_orientation0: ', shape(crystallite_orientation0) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_rotation: ', shape(crystallite_rotation) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_disorientation: ', shape(crystallite_disorientation) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_dt: ', shape(crystallite_dt) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subdt: ', shape(crystallite_subdt) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFrac: ', shape(crystallite_subFrac) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_subStep: ', shape(crystallite_subStep) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_localPlasticity: ', shape(crystallite_localPlasticity) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_requested: ', shape(crystallite_requested) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_todo: ', shape(crystallite_todo) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_converged: ', shape(crystallite_converged) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_sizePostResults: ', shape(crystallite_sizePostResults) - write(6,'(a35,1x,7(i8,1x))') 'crystallite_sizePostResult: ', shape(crystallite_sizePostResult) - write(6,'(/,a35,1x,i10)') 'Number of nonlocal grains: ',count(.not. crystallite_localPlasticity) + 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: ', nMax + 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 @@ -527,6 +449,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) subStepSizeCryst, & stepIncreaseCryst, & numerics_timeSyncing +#ifdef DEBUG use debug, only: & debug_level, & debug_crystallite, & @@ -536,6 +459,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) debug_e, & debug_i, & debug_g +#endif use IO, only: & IO_warning, & IO_error @@ -594,7 +518,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) neighboring_i, & o, & p, & - myNcomponents, & mySource ! local variables used for calculating analytic Jacobian real(pReal), dimension(3,3) :: temp_33 @@ -613,7 +536,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) real(pReal), dimension(9,9):: temp_99 logical :: error - +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt & .and. FEsolving_execElem(1) <= debug_e & .and. debug_e <= FEsolving_execElem(2)) then @@ -632,15 +555,15 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> 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 PRIVATE(myNcomponents) + !$OMP PARALLEL DO elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) - do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1_pInt,myNcomponents + do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) 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)) @@ -656,9 +579,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_dPdF0(1:3,1:3,1:3,1:3,c,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,c,i,e) ! ...stiffness crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) ! ...def grad crystallite_subTstar0_v(1:6,c,i,e) = crystallite_partionedTstar0_v(1:6,c,i,e) !...2nd PK stress - crystallite_subFe0(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_subF0(1:3,1:3,c,i,e), & - math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))), & - math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)))! only needed later on for stiffness calculation crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subStep(c,i,e) = 1.0_pReal/subStepSizeCryst crystallite_todo(c,i,e) = .true. @@ -679,10 +599,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) NiterationCrystallite = 0_pInt cutbackLooping: do while (any(crystallite_todo(:,startIP:endIP,FEsolving_execELem(1):FEsolving_execElem(2)))) - +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> crystallite iteration ',NiterationCrystallite - +#endif timeSyncing1: if (any(.not. crystallite_localPlasticity) .and. numerics_timeSyncing) then ! Time synchronization can only be used for nonlocal calculations, and only there it makes sense. @@ -701,6 +621,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! and its not clear how to fix this, so all nonlocals become terminally ill. if (any(crystallite_syncSubFrac .and. .not. crystallite_converged(1,:,:))) then +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) @@ -709,6 +630,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo endif +#endif crystallite_syncSubFrac = .false. where(.not. crystallite_localPlasticity) crystallite_substep = 0.0_pReal @@ -723,8 +645,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMP END PARALLEL DO +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> time synchronization: wind forward' +#endif endif elseif (any(crystallite_syncSubFracCompleted)) then @@ -740,8 +664,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_clearToCutback(i,e) = crystallite_localPlasticity(1,i,e) .or. .not. crystallite_converged(1,i,e) enddo enddo +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> time synchronization: done, proceed with cutback' +#endif else ! Normal calculation. @@ -767,8 +693,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (all(crystallite_localPlasticity .or. crystallite_converged)) then if (all(crystallite_localPlasticity .or. crystallite_subStep + crystallite_subFrac >= 1.0_pReal)) then crystallite_clearToWindForward = .true. ! final wind forward +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> final wind forward' +#endif else !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -777,8 +705,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMP END PARALLEL DO +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> wind forward' +#endif endif else subFracIntermediate = maxval(crystallite_subFrac, mask=.not.crystallite_localPlasticity) @@ -864,8 +794,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMP END PARALLEL DO +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> time synchronization: cutback' +#endif else !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -874,8 +806,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo !$OMP END PARALLEL DO +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) & write(6,'(a,i6)') '<< CRYST >> cutback' +#endif endif endif endif @@ -899,30 +833,24 @@ subroutine crystallite_stressAndItsTangent(updateJaco) endif timeSyncing1 - !$OMP PARALLEL DO PRIVATE(myNcomponents,formerSubStep) + !$OMP PARALLEL DO PRIVATE(formerSubStep) elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do c = 1,myNcomponents + do c = 1,homogenization_Ngrains(mesh_element(3,e)) ! --- wind forward --- if (crystallite_converged(c,i,e) .and. crystallite_clearToWindForward(i,e)) then formerSubStep = crystallite_subStep(c,i,e) crystallite_subFrac(c,i,e) = crystallite_subFrac(c,i,e) + crystallite_subStep(c,i,e) - !$OMP FLUSH(crystallite_subFrac) crystallite_subStep(c,i,e) = min(1.0_pReal - crystallite_subFrac(c,i,e), & stepIncreaseCryst * crystallite_subStep(c,i,e)) - !$OMP FLUSH(crystallite_subStep) + if (crystallite_subStep(c,i,e) > 0.0_pReal) then crystallite_subF0(1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) ! ...def grad - !$OMP FLUSH(crystallite_subF0) crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp(1:3,1:3,c,i,e) ! ...plastic velocity gradient crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li(1:3,1:3,c,i,e) ! ...intermediate velocity gradient crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp(1:3,1:3,c,i,e) ! ...plastic def grad crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi(1:3,1:3,c,i,e) ! ...intermediate def grad - crystallite_subFe0(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(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)) ! only needed later on for stiffness calculation !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)) @@ -938,7 +866,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) else crystallite_todo(c,i,e) = .true. endif - !$OMP FLUSH(crystallite_todo) #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelBasic) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & @@ -949,7 +876,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) #endif else ! this crystallite just converged for the entire timestep crystallite_todo(c,i,e) = .false. ! so done here - !$OMP FLUSH(crystallite_todo) endif ! --- cutback --- @@ -960,15 +886,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) else crystallite_subStep(c,i,e) = subStepSizeCryst * crystallite_subStep(c,i,e) ! cut step in half and restore... endif - !$OMP FLUSH(crystallite_subStep) crystallite_Fp(1:3,1:3,c,i,e) = crystallite_subFp0(1:3,1:3,c,i,e) ! ...plastic def grad - !$OMP FLUSH(crystallite_Fp) crystallite_invFp(1:3,1:3,c,i,e) = math_inv33(crystallite_Fp(1:3,1:3,c,i,e)) - !$OMP FLUSH(crystallite_invFp) crystallite_Fi(1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) ! ...intermediate def grad - !$OMP FLUSH(crystallite_Fi) crystallite_invFi(1:3,1:3,c,i,e) = math_inv33(crystallite_Fi(1:3,1:3,c,i,e)) - !$OMP FLUSH(crystallite_invFi) crystallite_Lp(1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) ! ...plastic velocity grad crystallite_Li(1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e) ! ...intermediate velocity grad plasticState (phaseAt(c,i,e))%state( :,phasememberAt(c,i,e)) = & @@ -981,7 +902,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! 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) - !$OMP FLUSH(crystallite_todo) #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & @@ -1002,10 +922,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if (crystallite_todo(c,i,e) .and. (crystallite_clearToWindForward(i,e) .or. crystallite_clearToCutback(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_subStep(c,i,e) * (crystallite_partionedF(1:3,1:3,c,i,e) & - crystallite_partionedF0(1:3,1:3,c,i,e)) - !$OMP FLUSH(crystallite_subF) crystallite_Fe(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(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)) @@ -1021,11 +939,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) timeSyncing2: if(numerics_timeSyncing) then if (any(.not. crystallite_localPlasticity .and. .not. crystallite_todo .and. .not. crystallite_converged & .and. crystallite_subStep <= subStepMinCryst)) then ! no way of rescuing a nonlocal ip that violated the lower time step limit, ... +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - do c = 1,myNcomponents + do c = 1,homogenization_Ngrains(mesh_element(3,e)) if (.not. crystallite_localPlasticity(c,i,e) .and. .not. crystallite_todo(c,i,e) & .and. .not. crystallite_converged(c,i,e) .and. crystallite_subStep(c,i,e) <= subStepMinCryst) & write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> nonlocal violated minimum subStep at el ip ipc ',e,i,c @@ -1033,6 +951,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo enddo elementLooping4 endif +#endif where(.not. crystallite_localPlasticity) crystallite_todo = .false. ! ... so let all nonlocal ips die peacefully crystallite_subStep = 0.0_pReal @@ -1040,6 +959,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) endif endif timeSyncing2 +#ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt) then write(6,'(/,a,f8.5)') '<< CRYST >> min(subStep) ',minval(crystallite_subStep) write(6,'(a,f8.5)') '<< CRYST >> max(subStep) ',maxval(crystallite_subStep) @@ -1052,9 +972,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) flush(6) endif endif +#endif ! --- integrate --- requires fully defined state array (basic + dependent state) - if (any(crystallite_todo)) call integrateState() where(.not. crystallite_converged .and. crystallite_subStep > subStepMinCryst) & ! do not try non-converged & fully cutbacked any further crystallite_todo = .true. @@ -1067,13 +987,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! --+>> CHECK FOR NON-CONVERGED CRYSTALLITES <<+-- elementLooping5: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do c = 1,myNcomponents + do c = 1,homogenization_Ngrains(mesh_element(3,e)) if (.not. crystallite_converged(c,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway) +#ifdef DEBUG if(iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> no convergence: respond fully elastic at el (elFE) ip ipc ', & e,'(',mesh_element(1,e),')',i,c +#endif invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,c,i,e)) Fe_guess = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,c,i,e), invFp), & math_inv33(crystallite_partionedFi0(1:3,1:3,c,i,e))) @@ -1081,6 +1002,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_P(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,c,i,e), invFp), & math_mul33x33(Tstar,transpose(invFp))) endif +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. c == debug_g) & .or. .not. iand(debug_level(debug_crystallite),debug_levelSelective) /= 0_pInt)) then @@ -1097,6 +1019,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) transpose(crystallite_Li(1:3,1:3,c,i,e)) flush(6) endif +#endif enddo enddo enddo elementLooping5 @@ -1106,11 +1029,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco) computeJacobian: if(updateJaco) then !$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,& - !$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,myNcomponents,error) + !$OMP rhs_3333,lhs_3333,temp_99,temp_33,temp_3333,error) elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) - myNcomponents = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - do c = 1_pInt,myNcomponents + do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) call constitutive_SandItsTangents(temp_33,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 @@ -1215,7 +1137,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo elementLooping6 !$OMP END PARALLEL DO endif computeJacobian -!why not OMP? end subroutine crystallite_stressAndItsTangent @@ -1223,20 +1144,20 @@ end subroutine crystallite_stressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 4th order explicit Runge Kutta method !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateRK4() +subroutine integrateStateRK4() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use FEsolving, only: & FEsolving_execElem, & FEsolving_execIP @@ -1404,7 +1325,7 @@ subroutine crystallite_integrateStateRK4() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1435,7 +1356,7 @@ subroutine crystallite_integrateStateRK4() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e,timeStepFraction(n)) ! fraction of original times step + crystallite_todo(g,i,e) = integrateStress(g,i,e,timeStepFraction(n)) ! fraction of original times step !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1506,27 +1427,27 @@ subroutine crystallite_integrateStateRK4() endif endif -end subroutine crystallite_integrateStateRK4 +end subroutine integrateStateRK4 !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 5th order Runge-Kutta Cash-Karp method with !> adaptive step size (use 5th order solution to advance = "local extrapolation") !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateRKCK45() +subroutine integrateStateRKCK45() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use numerics, only: & rTol_crystalliteState use FEsolving, only: & @@ -1602,8 +1523,10 @@ subroutine crystallite_integrateStateRKCK45() singleRun ! flag indicating computation for single (g,i,e) triple eIter = FEsolving_execElem(1:2) +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,1x,i1)') '<< CRYST >> Runge--Kutta step',1 +#endif ! --- LOOP ITERATOR FOR ELEMENT, GRAIN, IP --- do e = eIter(1),eIter(2) @@ -1723,7 +1646,7 @@ subroutine crystallite_integrateStateRKCK45() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1754,7 +1677,7 @@ subroutine crystallite_integrateStateRKCK45() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e,C(stage)) ! fraction of original time step + crystallite_todo(g,i,e) = integrateStress(g,i,e,C(stage)) ! fraction of original time step !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1901,9 +1824,6 @@ subroutine crystallite_integrateStateRKCK45() relSourceStateResiduum(s,mySource,g,i,e) = & sourceStateResiduum(s,mySource,g,i,e) / sourceState(p)%p(mySource)%state(s,cc) enddo - !$OMP FLUSH(relPlasticStateResiduum) - !$OMP FLUSH(relSourceStateResiduum) -! @Martin: do we need flushing? why..? crystallite_todo(g,i,e) = all(abs(relPlasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & rTol_crystalliteState .or. & abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & @@ -1943,7 +1863,7 @@ subroutine crystallite_integrateStateRKCK45() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1974,7 +1894,7 @@ subroutine crystallite_integrateStateRKCK45() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) + crystallite_todo(g,i,e) = integrateStress(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1998,32 +1918,33 @@ subroutine crystallite_integrateStateRKCK45() ! --- nonlocal convergence check --- - +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' ! if not requesting Integration of just a single IP +#endif if ((.not. singleRun) .and. any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged -end subroutine crystallite_integrateStateRKCK45 +end subroutine integrateStateRKCK45 !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateAdaptiveEuler() +subroutine integrateStateAdaptiveEuler() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use numerics, only: & rTol_crystalliteState use FEsolving, only: & @@ -2169,7 +2090,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -2201,7 +2122,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) + crystallite_todo(g,i,e) = integrateStress(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -2274,8 +2195,6 @@ subroutine crystallite_integrateStateAdaptiveEuler() + 0.5_pReal * sourceState(p)%p(mySource)%dotState(:,c) & * crystallite_subdt(g,i,e) ! contribution to absolute residuum in state enddo - !$OMP FLUSH(plasticStateResiduum) - !$OMP FLUSH(sourceStateResiduum) ! --- relative residui --- forall (s = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%dotState(s,c)) > 0.0_pReal) & @@ -2287,11 +2206,8 @@ subroutine crystallite_integrateStateAdaptiveEuler() relSourceStateResiduum(s,mySource,g,i,e) = & sourceStateResiduum(s,mySource,g,i,e) / sourceState(p)%p(mySource)%dotState(s,c) enddo - !$OMP FLUSH(relPlasticStateResiduum) - !$OMP FLUSH(relSourceStateResiduum) #ifdef DEBUG - if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)& .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then @@ -2319,13 +2235,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() abs(sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e)) < & sourceState(p)%p(mySource)%aTolState(1:mySizeSourceDotState)) enddo - if (converged) then - crystallite_converged(g,i,e) = .true. ! ... converged per definitionem - if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then - !$OMP CRITICAL (distributionState) - !$OMP END CRITICAL (distributionState) - endif - endif + if (converged) crystallite_converged(g,i,e) = .true. ! ... converged per definitionem endif enddo; enddo; enddo !$OMP ENDDO @@ -2333,33 +2243,34 @@ subroutine crystallite_integrateStateAdaptiveEuler() ! --- NONLOCAL CONVERGENCE CHECK --- - +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' +#endif if ((.not. singleRun) .and. any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged -end subroutine crystallite_integrateStateAdaptiveEuler +end subroutine integrateStateAdaptiveEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, and state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateEuler() +subroutine integrateStateEuler() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level, & debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use numerics, only: & numerics_timeSyncing use FEsolving, only: & @@ -2487,7 +2398,7 @@ eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local... .and. .not. numerics_timeSyncing) then @@ -2521,7 +2432,7 @@ eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) + crystallite_todo(g,i,e) = integrateStress(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e) & ! if broken non-local... .and. .not. numerics_timeSyncing) then @@ -2553,27 +2464,27 @@ eIter = FEsolving_execElem(1:2) crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged endif -end subroutine crystallite_integrateStateEuler +end subroutine integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine crystallite_integrateStateFPI() +subroutine integrateStateFPI() use, intrinsic :: & IEEE_arithmetic - use debug, only: & #ifdef DEBUG + use debug, only: & debug_e, & debug_i, & debug_g, & -#endif debug_level,& debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use numerics, only: & nState, & rTol_crystalliteState @@ -2637,8 +2548,10 @@ subroutine crystallite_integrateStateFPI() singleRun = (eIter(1) == eIter(2) .and. iIter(1,eIter(1)) == iIter(2,eIter(2))) +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo at start of state integration' +#endif !-------------------------------------------------------------------------------------------------- ! initialize dotState @@ -2693,8 +2606,10 @@ subroutine crystallite_integrateStateFPI() NaN = NaN .or. any(IEEE_is_NaN(sourceState(p)%p(mySource)%dotState(:,c))) enddo if (NaN) then ! NaN occured in any dotState +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,*) '<< CRYST >> dotstate ',plasticState(p)%dotState(:,c) +#endif if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... !$OMP CRITICAL (checkTodo) crystallite_todo = crystallite_todo .and. crystallite_localPlasticity ! ...all non-locals done (and broken) @@ -2708,9 +2623,10 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO ! --- UPDATE STATE --- +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after preguess of state' - +#endif !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains @@ -2766,14 +2682,15 @@ subroutine crystallite_integrateStateFPI() ! --- STRESS INTEGRATION --- +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo before stress integration' - +#endif !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - crystallite_todo(g,i,e) = crystallite_integrateStress(g,i,e) + crystallite_todo(g,i,e) = integrateStress(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e) .and. .not. crystallite_localPlasticity(g,i,e)) then ! broken non-local... !$OMP CRITICAL (checkTodo) @@ -2784,13 +2701,10 @@ subroutine crystallite_integrateStateFPI() enddo; enddo; enddo !$OMP ENDDO - !$OMP SINGLE - !$OMP CRITICAL (write2out) +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_todo(:,:,:)),' grains todo after stress integration' - !$OMP END CRITICAL (write2out) - !$OMP END SINGLE - +#endif ! --- DOT STATE --- @@ -2964,7 +2878,7 @@ subroutine crystallite_integrateStateFPI() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. crystallite_converged(g,i,e)) then ! converged and still alive... - crystallite_todo(g,i,e) = crystallite_stateJump(g,i,e) + crystallite_todo(g,i,e) = stateJump(g,i,e) !$OMP FLUSH(crystallite_todo) if (.not. crystallite_todo(g,i,e)) then ! if state jump fails, then convergence is broken crystallite_converged(g,i,e) = .false. @@ -2979,10 +2893,11 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO !$OMP END PARALLEL +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & write(6,'(a,i8,a,i2)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), & ' grains converged after state integration #', NiterationState - +#endif ! --- NON-LOCAL CONVERGENCE CHECK --- @@ -2991,12 +2906,15 @@ subroutine crystallite_integrateStateFPI() crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged endif +#ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) then write(6,'(a,i8,a)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), & ' grains converged after non-local check' write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_todo(:,:,:)), & ' grains todo after state integration #', NiterationState endif +#endif + ! --- CHECK IF DONE WITH INTEGRATION --- doneWithIntegration = .true. @@ -3010,14 +2928,14 @@ subroutine crystallite_integrateStateFPI() enddo elemLoop enddo crystalliteLooping -end subroutine crystallite_integrateStateFPI +end subroutine integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief calculates a jump in the state according to the current state and the current stress !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state !-------------------------------------------------------------------------------------------------- -logical function crystallite_stateJump(ipc,ip,el) +logical function stateJump(ipc,ip,el) use, intrinsic :: & IEEE_arithmetic use prec, only: & @@ -3067,7 +2985,7 @@ logical function crystallite_stateJump(ipc,ip,el) mySizePlasticDeltaState = plasticState(p)%sizeDeltaState if( any(IEEE_is_NaN(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c)))) then ! NaN occured in deltaState - crystallite_stateJump = .false. + stateJump = .false. return endif @@ -3081,7 +2999,7 @@ logical function crystallite_stateJump(ipc,ip,el) myOffsetSourceDeltaState = sourceState(p)%p(mySource)%offsetDeltaState mySizeSourceDeltaState = sourceState(p)%p(mySource)%sizeDeltaState if (any(IEEE_is_NaN(sourceState(p)%p(mySource)%deltaState(1:mySizeSourceDeltaState,c)))) then ! NaN occured in deltaState - crystallite_stateJump = .false. + stateJump = .false. return endif sourceState(p)%p(mySource)%state(myOffsetSourceDeltaState + 1_pInt : & @@ -3104,9 +3022,9 @@ logical function crystallite_stateJump(ipc,ip,el) endif #endif - crystallite_stateJump = .true. + stateJump = .true. -end function crystallite_stateJump +end function stateJump !-------------------------------------------------------------------------------------------------- @@ -3140,7 +3058,7 @@ end function crystallite_push33ToRef !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction !-------------------------------------------------------------------------------------------------- -logical function crystallite_integrateStress(& +logical function integrateStress(& ipc,& ! grain number ip,& ! integration point number el,& ! element number @@ -3157,16 +3075,16 @@ logical function crystallite_integrateStress(& iJacoLpresiduum, & subStepSizeLp, & subStepSizeLi - use debug, only: debug_level, & #ifdef DEBUG + use debug, only: debug_level, & debug_e, & debug_i, & debug_g, & -#endif debug_crystallite, & debug_levelBasic, & debug_levelExtensive, & debug_levelSelective +#endif use constitutive, only: constitutive_LpAndItsTangents, & constitutive_LiAndItsTangents, & @@ -3258,7 +3176,7 @@ logical function crystallite_integrateStress(& dgesv !* be pessimistic - crystallite_integrateStress = .false. + integrateStress = .false. #ifdef DEBUG if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & @@ -3597,7 +3515,7 @@ logical function crystallite_integrateStress(& !* set return flag to true - crystallite_integrateStress = .true. + integrateStress = .true. #ifdef DEBUG if (iand(debug_level(debug_crystallite),debug_levelExtensive) /= 0_pInt & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & @@ -3612,7 +3530,7 @@ logical function crystallite_integrateStress(& endif #endif -end function crystallite_integrateStress +end function integrateStress !-------------------------------------------------------------------------------------------------- @@ -3728,9 +3646,7 @@ function crystallite_postResults(ipc, ip, el) math_det33, & math_I3, & inDeg, & - math_Mandel6to33, & - math_qMul, & - math_qConj + math_Mandel6to33 use mesh, only: & mesh_element, & mesh_ipVolume, & @@ -3761,10 +3677,6 @@ function crystallite_postResults(ipc, ip, el) 1+plasticState(material_phase(ipc,ip,el))%sizePostResults + & sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: & crystallite_postResults - real(pReal), dimension(3,3) :: & - Ee - real(pReal), dimension(4) :: & - rotation real(pReal) :: & detF integer(pInt) :: & @@ -3808,18 +3720,6 @@ function crystallite_postResults(ipc, ip, el) crystallite_postResults(c+1:c+mySize) = & math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree - case (grainrotationx_ID) - mySize = 1_pInt - rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates - crystallite_postResults(c+1) = inDeg * rotation(1) * rotation(4) ! angle in degree - case (grainrotationy_ID) - mySize = 1_pInt - rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates - crystallite_postResults(c+1) = inDeg * rotation(2) * rotation(4) ! angle in degree - case (grainrotationz_ID) - mySize = 1_pInt - rotation = math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates - crystallite_postResults(c+1) = inDeg * rotation(3) * rotation(4) ! angle in degree ! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33 ! thus row index i is slow, while column index j is fast. reminder: "row is slow" @@ -3828,20 +3728,10 @@ function crystallite_postResults(ipc, ip, el) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)),[mySize]) - case (e_ID) - mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( & - transpose(crystallite_partionedF(1:3,1:3,ipc,ip,el)), & - crystallite_partionedF(1:3,1:3,ipc,ip,el)) - math_I3),[mySize]) case (fe_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & reshape(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)),[mySize]) - case (ee_ID) - Ee = 0.5_pReal *(math_mul33x33(transpose(crystallite_Fe(1:3,1:3,ipc,ip,el)), & - crystallite_Fe(1:3,1:3,ipc,ip,el)) - math_I3) - mySize = 9_pInt - crystallite_postResults(c+1:c+mySize) = reshape(Ee,[mySize]) case (fp_ID) mySize = 9_pInt crystallite_postResults(c+1:c+mySize) = & @@ -3887,8 +3777,8 @@ function crystallite_postResults(ipc, ip, el) c = c + 1_pInt if (size(crystallite_postResults)-c > 0_pInt) & crystallite_postResults(c+1:size(crystallite_postResults)) = & - constitutive_postResults(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe, & - ipc, ip, el) + constitutive_postResults(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fi(1:3,1:3,ipc,ip,el), & + crystallite_Fe, ipc, ip, el) end function crystallite_postResults diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 77d301400..82a97dc53 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -448,8 +448,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) subStepSizeHomog, & stepIncreaseHomog, & nMPstate - use math, only: & - math_transpose33 use FEsolving, only: & FEsolving_execElem, & FEsolving_execIP, & @@ -496,6 +494,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) crystallite_converged, & crystallite_stressAndItsTangent, & crystallite_orientations +#ifdef DEBUG use debug, only: & debug_level, & debug_homogenization, & @@ -504,6 +503,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) debug_levelSelective, & debug_e, & debug_i +#endif implicit none real(pReal), intent(in) :: dt !< time increment @@ -517,18 +517,16 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) mySource, & myNgrains -!-------------------------------------------------------------------------------------------------- -! initialize to starting condition +#ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then - !$OMP CRITICAL (write2out) write(6,'(/a,i5,1x,i2)') '<< HOMOG >> Material Point start at el ip ', debug_e, debug_i write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F0', & - math_transpose33(materialpoint_F0(1:3,1:3,debug_i,debug_e)) + transpose(materialpoint_F0(1:3,1:3,debug_i,debug_e)) write(6,'(a,/,3(12x,3(f14.9,1x)/))') '<< HOMOG >> F', & - math_transpose33(materialpoint_F(1:3,1:3,debug_i,debug_e)) - !$OMP END CRITICAL (write2out) + transpose(materialpoint_F(1:3,1:3,debug_i,debug_e)) endif +#endif !-------------------------------------------------------------------------------------------------- ! initialize restoration points of ... @@ -608,10 +606,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) !--------------------------------------------------------------------------------------------------- ! calculate new subStep and new subFrac materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e) - !$OMP FLUSH(materialpoint_subFrac) materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), & stepIncreaseHomog*materialpoint_subStep(i,e)) ! introduce flexibility for step increase/acceleration - !$OMP FLUSH(materialpoint_subStep) steppingNeeded: if (materialpoint_subStep(i,e) > subStepMinHomog) then @@ -671,7 +667,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) hydrogenfluxState(mappingHomogenization(2,i,e))%subState0(:,mappingHomogenization(1,i,e)) = & hydrogenfluxState(mappingHomogenization(2,i,e))%State( :,mappingHomogenization(1,i,e))! ...internal hydrogen transport state materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad - !$OMP FLUSH(materialpoint_subF0) endif steppingNeeded else converged @@ -689,7 +684,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) !$OMP END CRITICAL (setTerminallyIll) else ! cutback makes sense materialpoint_subStep(i,e) = subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback - !$OMP FLUSH(materialpoint_subStep) #ifdef DEBUG if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt & @@ -752,8 +746,9 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) if (materialpoint_subStep(i,e) > subStepMinHomog) then materialpoint_requested(i,e) = .true. - materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) + & - materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) - materialpoint_F0(1:3,1:3,i,e)) + materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) & + + materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) & + - materialpoint_F0(1:3,1:3,i,e)) materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt materialpoint_doneAndHappy(1:2,i,e) = [.false.,.true.] endif @@ -810,13 +805,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) materialpoint_doneAndHappy(1:2,i,e) = homogenization_updateState(i,e) materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy endif - !$OMP FLUSH(materialpoint_converged) - if (materialpoint_converged(i,e)) then - if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0_pInt) then - !$OMP CRITICAL (distributionMPState) - !$OMP END CRITICAL (distributionMPState) - endif - endif endif enddo IpLooping3 enddo elementLooping3 @@ -838,9 +826,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) enddo elementLooping4 !$OMP END PARALLEL DO else - !$OMP CRITICAL (write2out) write(6,'(/,a,/)') '<< HOMOG >> Material Point terminally ill' - !$OMP END CRITICAL (write2out) endif end subroutine materialpoint_stressAndItsTangent diff --git a/src/homogenization_RGC.f90 b/src/homogenization_RGC.f90 index 92ea5301d..1d7bc6f86 100644 --- a/src/homogenization_RGC.f90 +++ b/src/homogenization_RGC.f90 @@ -942,10 +942,10 @@ pure function homogenization_RGC_postResults(ip,el,avgP,avgF) do o = 1_pInt,homogenization_Noutput(mesh_element(3,el)) select case(homogenization_RGC_outputID(o,homID)) case (avgdefgrad_ID) - homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(avgF,[9]) + homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgF),[9]) c = c + 9_pInt case (avgfirstpiola_ID) - homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(avgP,[9]) + homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgP),[9]) c = c + 9_pInt case (ipcoords_ID) homogenization_RGC_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates diff --git a/src/homogenization_isostrain.f90 b/src/homogenization_isostrain.f90 index 8ee0df73d..24aedf75f 100644 --- a/src/homogenization_isostrain.f90 +++ b/src/homogenization_isostrain.f90 @@ -300,10 +300,10 @@ pure function homogenization_isostrain_postResults(ip,el,avgP,avgF) homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal) c = c + 1_pInt case (avgdefgrad_ID) - homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(avgF,[9]) + homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgF),[9]) c = c + 9_pInt case (avgfirstpiola_ID) - homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(avgP,[9]) + homogenization_isostrain_postResults(c+1_pInt:c+9_pInt) = reshape(transpose(avgP),[9]) c = c + 9_pInt case (ipcoords_ID) homogenization_isostrain_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 59a106435..f7c723521 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -57,17 +57,23 @@ module plastic_phenopowerlaw Nslip, & !< active number of slip systems per family Ntwin !< active number of twin systems per family real(pReal), dimension(:), allocatable :: & - tau0_slip, & !< initial critical shear stress for slip - tau0_twin, & !< initial critical shear stress for twin - tausat_slip, & !< maximum critical shear stress for slip + xi_slip_0, & !< initial critical shear stress for slip + xi_twin_0, & !< initial critical shear stress for twin + xi_slip_sat, & !< maximum critical shear stress for slip nonSchmidCoeff, & - H_int !< per family hardening activity (optional) + H_int, & !< per family hardening activity (optional) !ToDo: Better name! + gamma_twin_char !< characteristic shear for twins real(pReal), dimension(:,:), allocatable :: & interaction_SlipSlip, & !< slip resistance from slip activity interaction_SlipTwin, & !< slip resistance from twin activity interaction_TwinSlip, & !< twin resistance from slip activity interaction_TwinTwin !< twin resistance from twin activity - + real(pReal), dimension(:,:,:), allocatable :: & + Schmid_slip, & + Schmid_twin + real(pReal), dimension(:,:,:,:), allocatable :: & + nonSchmid_pos, & + nonSchmid_neg integer(kind(undefined_ID)), dimension(:), allocatable :: & outputID !< ID of each post result output end type @@ -76,10 +82,10 @@ module plastic_phenopowerlaw type, private :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & - s_slip, & - s_twin, & - accshear_slip, & - accshear_twin, & + xi_slip, & + xi_twin, & + gamma_slip, & + gamma_twin, & whole real(pReal), pointer, dimension(:) :: & sumGamma, & @@ -116,8 +122,6 @@ subroutine plastic_phenopowerlaw_init debug_constitutive,& debug_levelBasic use math, only: & - math_Mandel3333to66, & - math_Voigt66to3333, & math_expand use IO, only: & IO_warning, & @@ -127,7 +131,7 @@ subroutine plastic_phenopowerlaw_init phase_plasticity, & phase_plasticityInstance, & phase_Noutput, & - PLASTICITY_PHENOPOWERLAW_label, & + PLASTICITY_PHENOPOWERLAW_LABEL, & PLASTICITY_PHENOPOWERLAW_ID, & material_phase, & plasticState @@ -154,9 +158,13 @@ subroutine plastic_phenopowerlaw_init real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] - type(tParameters) :: prm + type(tParameters) :: & + prm + type(tPhenopowerlawState) :: & + stt, & + dot - integer(kind(undefined_ID)) :: & + integer(kind(undefined_ID)) :: & outputID !< ID of each post result output character(len=512) :: & @@ -182,57 +190,104 @@ subroutine plastic_phenopowerlaw_init do p = 1_pInt, size(phase_plasticityInstance) if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle instance = phase_plasticityInstance(p) - associate(prm => param(instance)) + associate(prm => param(instance),stt => state(instance),dot => dotState(instance)) + extmsg = '' - prm%Nslip = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray) - if (size(prm%Nslip) > count(lattice_NslipSystem(:,p) > 0_pInt)) call IO_error(150_pInt,ext_msg='Nslip') - if (any(lattice_NslipSystem(1:size(prm%Nslip),p)-prm%Nslip < 0_pInt)) call IO_error(150_pInt,ext_msg='Nslip') + prm%Nslip = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray) prm%totalNslip = sum(prm%Nslip) + if (size(prm%Nslip) > count(lattice_NslipSystem(:,p) > 0_pInt)) & + call IO_error(150_pInt,ext_msg='Nslip') + if (any(lattice_NslipSystem(1:size(prm%Nslip),p)-prm%Nslip < 0_pInt)) & + call IO_error(150_pInt,ext_msg='Nslip') - if (prm%totalNslip > 0_pInt) then - prm%tau0_slip = config_phase(p)%getFloats('tau0_slip') - prm%tausat_slip = config_phase(p)%getFloats('tausat_slip') - prm%interaction_SlipSlip = spread(config_phase(p)%getFloats('interaction_slipslip'),2,1) - prm%H_int = config_phase(p)%getFloats('h_int',& - defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) + slipActive: if (prm%totalNslip > 0_pInt) then + ! reading in slip related parameters + prm%xi_slip_0 = config_phase(p)%getFloats('tau0_slip', requiredShape=shape(prm%Nslip)) + prm%xi_slip_sat = config_phase(p)%getFloats('tausat_slip', requiredShape=shape(prm%Nslip)) + prm%interaction_SlipSlip = spread(config_phase(p)%getFloats('interaction_slipslip', & + requiredShape=shape(prm%Nslip)),2,1) + prm%H_int = config_phase(p)%getFloats('h_int', requiredShape=shape(prm%Nslip), & + defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',& - defaultVal = emptyRealArray ) + defaultVal = emptyRealArray ) prm%gdot0_slip = config_phase(p)%getFloat('gdot0_slip') prm%n_slip = config_phase(p)%getFloat('n_slip') prm%a_slip = config_phase(p)%getFloat('a_slip') prm%h0_SlipSlip = config_phase(p)%getFloat('h0_slipslip') - endif - prm%Ntwin = config_phase(p)%getInts('ntwin', defaultVal=emptyIntArray) - if (size(prm%Ntwin) > count(lattice_NtwinSystem(:,p) > 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') - if (any(lattice_NtwinSystem(1:size(prm%Ntwin),p)-prm%Ntwin < 0_pInt)) call IO_error(150_pInt,ext_msg='Ntwin') + ! sanity checks for slip related parameters + if (any(prm%xi_slip_0 < 0.0_pReal .and. prm%Nslip > 0_pInt)) & + extmsg = trim(extmsg)//"xi_slip_0 " + if (any(prm%xi_slip_sat < prm%xi_slip_0 .and. prm%Nslip > 0_pInt)) & + extmsg = trim(extmsg)//"xi_slip_sat " + + if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//"gdot0_slip " + if (dEq0(prm%a_slip)) extmsg = trim(extmsg)//"a_slip " ! ToDo: negative values ok? + if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//"n_slip " ! ToDo: negative values ok? + + ! expand slip related parameters from system => family + prm%xi_slip_0 = math_expand(prm%xi_slip_0,prm%Nslip) + prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip) + prm%H_int = math_expand(prm%H_int,prm%Nslip) + else slipActive + allocate(prm%xi_slip_0(0)) + endif slipActive + + prm%Ntwin = config_phase(p)%getInts('ntwin', defaultVal=emptyIntArray) prm%totalNtwin = sum(prm%Ntwin) + if (size(prm%Ntwin) > count(lattice_NtwinSystem(:,p) > 0_pInt)) & + call IO_error(150_pInt,ext_msg='Ntwin') + if (any(lattice_NtwinSystem(1:size(prm%Ntwin),p)-prm%Ntwin < 0_pInt)) & + call IO_error(150_pInt,ext_msg='Ntwin') - if (prm%totalNtwin > 0_pInt) then - prm%tau0_twin = config_phase(p)%getFloats('tau0_twin') - prm%interaction_TwinTwin = spread(config_phase(p)%getFloats('interaction_twintwin'),2,1) + twinActive: if (prm%totalNtwin > 0_pInt) then + ! reading in twin related parameters + prm%xi_twin_0 = config_phase(p)%getFloats('tau0_twin',requiredShape=shape(prm%Ntwin)) + prm%interaction_TwinTwin = spread(config_phase(p)%getFloats('interaction_twintwin', & + requiredShape=shape(prm%Ntwin)),2,1) prm%gdot0_twin = config_phase(p)%getFloat('gdot0_twin') prm%n_twin = config_phase(p)%getFloat('n_twin') prm%spr = config_phase(p)%getFloat('s_pr') - prm%twinB = config_phase(p)%getFloat('twin_b') - prm%twinC = config_phase(p)%getFloat('twin_c') - prm%twinD = config_phase(p)%getFloat('twin_d') - prm%twinE = config_phase(p)%getFloat('twin_e') prm%h0_TwinTwin = config_phase(p)%getFloat('h0_twintwin') - endif - if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then + ! sanity checks for twin related parameters + if (any(prm%xi_twin_0 < 0.0_pReal .and. prm%Ntwin > 0_pInt)) & + extmsg = trim(extmsg)//"xi_twin_0 " + if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//"gdot0_twin " + if (dEq0(prm%n_twin)) extmsg = trim(extmsg)//"n_twin " ! ToDo: negative values ok? + + ! expand slip related parameters from system => family + prm%xi_twin_0 = math_expand(prm%xi_twin_0,prm%Ntwin) + else twinActive + allocate(prm%xi_twin_0(0)) + endif twinActive + + slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then prm%interaction_SlipTwin = spread(config_phase(p)%getFloats('interaction_sliptwin'),2,1) prm%interaction_TwinSlip = spread(config_phase(p)%getFloats('interaction_twinslip'),2,1) prm%h0_TwinSlip = config_phase(p)%getFloat('h0_twinslip') - endif + else slipAndTwinActive + prm%h0_TwinSlip = 0.0_pReal + endif slipAndTwinActive + ! optional parameters that should be defined + prm%twinB = config_phase(p)%getFloat('twin_b',defaultVal=1.0_pReal) + prm%twinC = config_phase(p)%getFloat('twin_c',defaultVal=0.0_pReal) + prm%twinD = config_phase(p)%getFloat('twin_d',defaultVal=0.0_pReal) + prm%twinE = config_phase(p)%getFloat('twin_e',defaultVal=0.0_pReal) prm%aTolResistance = config_phase(p)%getFloat('atol_resistance',defaultVal=1.0_pReal) - prm%aTolShear = config_phase(p)%getFloat('atol_shear',defaultVal=1.0e-6_pReal) - prm%aTolTwinfrac = config_phase(p)%getFloat('atol_twinfrac',defaultVal=1.0e-6_pReal) + prm%aTolShear = config_phase(p)%getFloat('atol_shear', defaultVal=1.0e-6_pReal) + prm%aTolTwinfrac = config_phase(p)%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) + + if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//"aTolresistance " + if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//"aTolShear " + if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//"atoltwinfrac " + + if (extmsg /= '') call IO_error(211_pInt,ip=instance,& + ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') outputs = config_phase(p)%getStrings('(output)',defaultVal=emptyStringArray) allocate(prm%outputID(0)) @@ -240,36 +295,36 @@ subroutine plastic_phenopowerlaw_init outputID = undefined_ID select case(outputs(i)) case ('resistance_slip') - outputID = resistance_slip_ID + outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0_pInt) outputSize = prm%totalNslip case ('accumulatedshear_slip') - outputID = accumulatedshear_slip_ID + outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0_pInt) outputSize = prm%totalNslip case ('shearrate_slip') - outputID = shearrate_slip_ID + outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0_pInt) outputSize = prm%totalNslip case ('resolvedstress_slip') - outputID = resolvedstress_slip_ID + outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0_pInt) outputSize = prm%totalNslip case ('resistance_twin') - outputID = resistance_twin_ID + outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) outputSize = prm%totalNtwin case ('accumulatedshear_twin') - outputID = accumulatedshear_twin_ID + outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) outputSize = prm%totalNtwin case ('shearrate_twin') - outputID = shearrate_twin_ID + outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) outputSize = prm%totalNtwin case ('resolvedstress_twin') - outputID = resolvedstress_twin_ID + outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) outputSize = prm%totalNtwin - case ('totalvolfrac_twin') - outputID = totalvolfrac_twin_ID - outputSize = 1_pInt case ('totalshear') - outputID = totalshear_ID + outputID = merge(totalshear_ID,undefined_ID,prm%totalNslip>0_pInt) + outputSize = 1_pInt + case ('totalvolfrac_twin') + outputID = merge(totalvolfrac_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) outputSize = 1_pInt end select @@ -281,56 +336,19 @@ subroutine plastic_phenopowerlaw_init end do - extmsg = '' - if (sum(prm%Nslip) > 0_pInt) then - if (size(prm%tau0_slip) /= size(prm%Nslip)) call IO_error(211_pInt,ip=instance, & - ext_msg='shape(tau0_slip) ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (size(prm%tausat_slip) /= size(prm%Nslip)) call IO_error(211_pInt,ip=instance, & - ext_msg='shape(tausat_slip) ('//PLASTICITY_PHENOPOWERLAW_label//')') - if (size(prm%H_int) /= size(prm%Nslip)) call IO_error(211_pInt,ip=instance, & - ext_msg='shape(H_int) ('//PLASTICITY_PHENOPOWERLAW_label//')') - - if (any(prm%tau0_slip < 0.0_pReal .and. prm%Nslip > 0_pInt)) & - extmsg = trim(extmsg)//"tau0_slip " - if (any(prm%tausat_slip < prm%tau0_slip .and. prm%Nslip > 0_pInt)) & - extmsg = trim(extmsg)//"tausat_slip " - - if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//" gdot0_slip " - if (dEq0(prm%a_slip)) extmsg = trim(extmsg)//" a_slip " ! ToDo: negative values ok? - if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//" n_slip " ! ToDo: negative values ok? - endif - - if (sum(prm%Ntwin) > 0_pInt) then - if (size(prm%tau0_twin) /= size(prm%ntwin)) call IO_error(211_pInt,ip=instance,& - ext_msg='shape(tau0_twin) ('//PLASTICITY_PHENOPOWERLAW_label//')') - - if (any(prm%tau0_twin < 0.0_pReal .and. prm%Ntwin > 0_pInt)) & - extmsg = trim(extmsg)//"tau0_twin " - - if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//"gdot0_twin " - if (dEq0(prm%n_twin)) extmsg = trim(extmsg)//"n_twin " ! ToDo: negative values ok? - endif - - if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//"aTolresistance " - if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//"aTolShear " - if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//"atoltwinfrac " - - if (extmsg /= '') call IO_error(211_pInt,ip=instance,& - ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') - !-------------------------------------------------------------------------------------------------- ! allocate state arrays - NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase - sizeState = size(['tau_slip ','accshear_slip']) * prm%TotalNslip & - + size(['tau_twin ','accshear_twin']) * prm%TotalNtwin & + NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase + sizeState = size(['tau_slip ','gamma_slip']) * prm%TotalNslip & + + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin & + size(['sum(gamma)', 'sum(f) ']) sizeDotState = sizeState plasticState(p)%sizeState = sizeState plasticState(p)%sizeDotState = sizeDotState plasticState(p)%sizePostResults = sum(plastic_phenopowerlaw_sizePostResult(:,instance)) - plasticState(p)%nSlip = sum(prm%Nslip) - plasticState(p)%nTwin = sum(prm%Ntwin) + plasticState(p)%nSlip = prm%totalNslip + plasticState(p)%nTwin = prm%totalNtwin allocate(plasticState(p)%aTolState ( sizeState), source=0.0_pReal) allocate(plasticState(p)%state0 ( sizeState,NipcMyPhase), source=0.0_pReal) allocate(plasticState(p)%partionedState0 ( sizeState,NipcMyPhase), source=0.0_pReal) @@ -351,12 +369,24 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! calculate hardening matrices - allocate(temp1(sum(prm%Nslip),sum(prm%Nslip)),source =0.0_pReal) - allocate(temp2(sum(prm%Nslip),sum(prm%Ntwin)),source =0.0_pReal) + allocate(temp1(prm%totalNslip,prm%totalNslip),source = 0.0_pReal) + allocate(temp2(prm%totalNslip,prm%totalNtwin),source = 0.0_pReal) + allocate(prm%Schmid_slip(3,3,prm%totalNslip),source = 0.0_pReal) + allocate(prm%nonSchmid_pos(3,3,size(prm%nonSchmidCoeff),prm%totalNslip),source = 0.0_pReal) + allocate(prm%nonSchmid_neg(3,3,size(prm%nonSchmidCoeff),prm%totalNslip),source = 0.0_pReal) + i = 0_pInt mySlipFamilies: do f = 1_pInt,size(prm%Nslip,1) ! >>> interaction slip -- X index_myFamily = sum(prm%Nslip(1:f-1_pInt)) mySlipSystems: do j = 1_pInt,prm%Nslip(f) + i = i + 1_pInt + prm%Schmid_slip(1:3,1:3,i) = lattice_Sslip(1:3,1:3,1,sum(lattice_Nslipsystem(1:f-1,p))+j,p) + do k = 1,size(prm%nonSchmidCoeff) + prm%nonSchmid_pos(1:3,1:3,k,i) = lattice_Sslip(1:3,1:3,2*k, index_myFamily+j,p) & + * prm%nonSchmidCoeff(k) + prm%nonSchmid_neg(1:3,1:3,k,i) = lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+j,p) & + * prm%nonSchmidCoeff(k) + enddo otherSlipFamilies: do o = 1_pInt,size(prm%Nslip,1) index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) otherSlipSystems: do k = 1_pInt,prm%Nslip(o) @@ -380,13 +410,19 @@ subroutine plastic_phenopowerlaw_init enddo mySlipFamilies prm%interaction_SlipSlip = temp1; deallocate(temp1) prm%interaction_SlipTwin = temp2; deallocate(temp2) - - allocate(temp1(sum(prm%Ntwin),sum(prm%Nslip)),source =0.0_pReal) - allocate(temp2(sum(prm%Ntwin),sum(prm%Ntwin)),source =0.0_pReal) + + allocate(temp1(prm%totalNtwin,prm%totalNslip),source = 0.0_pReal) + allocate(temp2(prm%totalNtwin,prm%totalNtwin),source = 0.0_pReal) + allocate(prm%Schmid_twin(3,3,prm%totalNtwin),source = 0.0_pReal) + allocate(prm%gamma_twin_char(prm%totalNtwin),source = 0.0_pReal) + i = 0_pInt myTwinFamilies: do f = 1_pInt,size(prm%Ntwin,1) ! >>> interaction twin -- X index_myFamily = sum(prm%Ntwin(1:f-1_pInt)) myTwinSystems: do j = 1_pInt,prm%Ntwin(f) + i = i + 1_pInt + prm%Schmid_twin(1:3,1:3,i) = lattice_Stwin(1:3,1:3,sum(lattice_NTwinsystem(1:f-1,p))+j,p) + prm%gamma_twin_char(i) = lattice_shearTwin(sum(lattice_Ntwinsystem(1:f-1,p))+j,p) slipFamilies: do o = 1_pInt,size(prm%Nslip,1) index_otherFamily = sum(prm%Nslip(1:o-1_pInt)) slipSystems: do k = 1_pInt,prm%Nslip(o) @@ -414,50 +450,49 @@ subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1_pInt - endIndex = plasticState(p)%nSlip - state (instance)%s_slip=>plasticState(p)%state (startIndex:endIndex,:) - dotState(instance)%s_slip=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%state0(startIndex:endIndex,:) = & - spread(math_expand(prm%tau0_slip, prm%Nslip), 2, NipcMyPhase) + endIndex = prm%totalNslip + stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) + stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) + dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance startIndex = endIndex + 1_pInt - endIndex = endIndex + plasticState(p)%nTwin - state (instance)%s_twin=>plasticState(p)%state (startIndex:endIndex,:) - dotState(instance)%s_twin=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%state0(startIndex:endIndex,:) = & - spread(math_expand(prm%tau0_twin, prm%Ntwin), 2, NipcMyPhase) + endIndex = endIndex + prm%totalNtwin + stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) + stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) + dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance startIndex = endIndex + 1_pInt endIndex = endIndex + 1_pInt - state (instance)%sumGamma=>plasticState(p)%state (startIndex,:) - dotState(instance)%sumGamma=>plasticState(p)%dotState(startIndex,:) + stt%sumGamma => plasticState(p)%state (startIndex,:) + dot%sumGamma => plasticState(p)%dotState(startIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear startIndex = endIndex + 1_pInt endIndex = endIndex + 1_pInt - state (instance)%sumF=>plasticState(p)%state (startIndex,:) - dotState(instance)%sumF=>plasticState(p)%dotState(startIndex,:) + stt%sumF=>plasticState(p)%state (startIndex,:) + dot%sumF=>plasticState(p)%dotState(startIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTwinFrac startIndex = endIndex + 1_pInt - endIndex = endIndex + plasticState(p)%nSlip - state (instance)%accshear_slip=>plasticState(p)%state (startIndex:endIndex,:) - dotState(instance)%accshear_slip=>plasticState(p)%dotState(startIndex:endIndex,:) + endIndex = endIndex + prm%totalNslip + stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) + dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear ! global alias - plasticState(p)%slipRate =>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%accumulatedSlip =>plasticState(p)%state(startIndex:endIndex,:) + plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) startIndex = endIndex + 1_pInt - endIndex = endIndex + plasticState(p)%nTwin - state (instance)%accshear_twin=>plasticState(p)%state (startIndex:endIndex,:) - dotState(instance)%accshear_twin=>plasticState(p)%dotState(startIndex:endIndex,:) + endIndex = endIndex + prm%totalNtwin + stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) + dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - - dotState(instance)%whole =>plasticState(p)%dotState - + + plasticState(p)%state0 = plasticState(p)%state + dot%whole => plasticState(p)%dotState + end associate enddo @@ -467,194 +502,90 @@ end subroutine plastic_phenopowerlaw_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) - use prec, only: & - dNeq0 - use math, only: & - math_Plain3333to99, & - math_Mandel6to33 - use lattice, only: & - lattice_Sslip, & - lattice_Sslip_v, & - lattice_Stwin, & - lattice_Stwin_v, & - lattice_NslipSystem, & - lattice_NtwinSystem - use material, only: & - phasememberAt, & - material_phase, & - phase_plasticityInstance +subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) implicit none real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient - real(pReal), dimension(9,9), intent(out) :: & - dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + instance, & + of integer(pInt) :: & - index_myFamily, & - f,i,j,k,l,m,n, & - of, & - ph - real(pReal) :: & - tau_slip_pos,tau_slip_neg, & - gdot_slip_pos,gdot_slip_neg, & + i,k,l,m,n + real(pReal), dimension(param(instance)%totalNslip) :: & dgdot_dtauslip_pos,dgdot_dtauslip_neg, & - gdot_twin,dgdot_dtautwin,tau_twin - real(pReal), dimension(3,3,3,3) :: & - dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor - real(pReal), dimension(3,3,2) :: & - nonSchmid_tensor + gdot_slip_pos,gdot_slip_neg + real(pReal), dimension(param(instance)%totalNtwin) :: & + gdot_twin,dgdot_dtautwin + type(tParameters) :: prm type(tPhenopowerlawState) :: stt - of = phasememberAt(ipc,ip,el) - ph = material_phase(ipc,ip,el) - - associate(prm => param(phase_plasticityInstance(ph)), stt => state(phase_plasticityInstance(ph))) + associate(prm => param(instance), stt => state(instance)) Lp = 0.0_pReal - dLp_dTstar3333 = 0.0_pReal - dLp_dTstar99 = 0.0_pReal + dLp_dMp = 0.0_pReal -!-------------------------------------------------------------------------------------------------- -! Slip part - j = 0_pInt - slipFamilies: do f = 1_pInt,size(prm%Nslip,1) - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems: do i = 1_pInt,prm%Nslip(f) - j = j+1_pInt - - ! Calculation of Lp - tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) - tau_slip_neg = tau_slip_pos - nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) - nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) - do k = 1,size(prm%nonSchmidCoeff) - tau_slip_pos = tau_slip_pos + prm%nonSchmidCoeff(k)* & - dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) - tau_slip_neg = tau_slip_neg + prm%nonSchmidCoeff(k)* & - dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) - nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + prm%nonSchmidCoeff(k)*& - lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph) - nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + prm%nonSchmidCoeff(k)*& - lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph) - enddo - gdot_slip_pos = 0.5_pReal*prm%gdot0_slip* & - ((abs(tau_slip_pos)/(stt%s_slip(j,of)))**prm%n_slip)*sign(1.0_pReal,tau_slip_pos) - - gdot_slip_neg = 0.5_pReal*prm%gdot0_slip* & - ((abs(tau_slip_neg)/(stt%s_slip(j,of)))**prm%n_slip)*sign(1.0_pReal,tau_slip_neg) - - Lp = Lp + (1.0_pReal-stt%sumF(of))*& - (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) - - ! Calculation of the tangent of Lp - if (dNeq0(tau_slip_pos)) then - dgdot_dtauslip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & - dgdot_dtauslip_pos*lattice_Sslip(k,l,1,index_myFamily+i,ph)* & - nonSchmid_tensor(m,n,1) - endif - - if (dNeq0(tau_slip_neg)) then - dgdot_dtauslip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & - dgdot_dtauslip_neg*lattice_Sslip(k,l,1,index_myFamily+i,ph)* & - nonSchmid_tensor(m,n,2) - endif - enddo slipSystems - enddo slipFamilies - -!-------------------------------------------------------------------------------------------------- -! Twinning part - j = 0_pInt - twinFamilies: do f = 1_pInt,size(prm%Ntwin,1) - index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems: do i = 1_pInt,prm%Ntwin(f) - j = j+1_pInt - - ! Calculation of Lp - tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) - gdot_twin = (1.0_pReal-stt%sumF(of))*prm%gdot0_twin*& - (abs(tau_twin)/stt%s_twin(j,of))**& - prm%n_twin*max(0.0_pReal,sign(1.0_pReal,tau_twin)) - Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph) - - ! Calculation of the tangent of Lp - if (dNeq0(gdot_twin)) then - dgdot_dtautwin = gdot_twin*prm%n_twin/tau_twin - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & - dgdot_dtautwin*lattice_Stwin(k,l,index_myFamily+i,ph)* & - lattice_Stwin(m,n,index_myFamily+i,ph) - endif - enddo twinSystems - enddo twinFamilies - - dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) + call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + slipSystems: do i = 1_pInt, prm%totalNslip + Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) & + *(prm%Schmid_slip(m,n,i) + sum(prm%nonSchmid_pos(m,n,:,i))) + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) & + *(prm%Schmid_slip(m,n,i) + sum(prm%nonSchmid_neg(m,n,:,i))) + enddo slipSystems + call kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtautwin) + twinSystems: do i = 1_pInt, prm%totalNtwin + Lp = Lp + gdot_twin(i)*prm%Schmid_twin(1:3,1:3,i) + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) + enddo twinSystems end associate + + end subroutine plastic_phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) - use lattice, only: & - lattice_Sslip_v, & - lattice_Stwin_v, & - lattice_NslipSystem, & - lattice_NtwinSystem, & - lattice_shearTwin - use material, only: & - material_phase, & - phasememberAt, & - phase_plasticityInstance +subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) implicit none - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element !< microstructure state + instance, & + of integer(pInt) :: & - ph, & - f,i,j,k, & - index_myFamily, & - of + i,k real(pReal) :: & c_SlipSlip,c_TwinSlip,c_TwinTwin, & - ssat_offset, & - tau_slip_pos,tau_slip_neg,tau_twin + xi_slip_sat_offset - real(pReal), dimension(param(phase_plasticityInstance(material_phase(ipc,ip,el)))%totalNslip) :: & - gdot_slip,left_SlipSlip,right_SlipSlip - real(pReal), dimension(param(phase_plasticityInstance(material_phase(ipc,ip,el)))%totalNtwin) :: & - gdot_twin + real(pReal), dimension(param(instance)%totalNslip) :: & + left_SlipSlip,right_SlipSlip, & + gdot_slip_pos,gdot_slip_neg type(tParameters) :: prm - type(tPhenopowerlawState) :: dst,stt + type(tPhenopowerlawState) :: dot,stt - of = phasememberAt(ipc,ip,el) - ph = material_phase(ipc,ip,el) - associate( prm => param(phase_plasticityInstance(ph)), & - stt => state(phase_plasticityInstance(ph)), & - dst => dotState(phase_plasticityInstance(ph))) + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - dst%whole(:,of) = 0.0_pReal + dot%whole(:,of) = 0.0_pReal !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices @@ -663,216 +594,245 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) c_TwinTwin = prm%h0_TwinTwin * stt%sumF(of)**prm%twinD !-------------------------------------------------------------------------------------------------- -! calculate left and right vectors and calculate dot gammas - ssat_offset = prm%spr*sqrt(stt%sumF(of)) - j = 0_pInt - slipFamilies1: do f =1_pInt,size(prm%Nslip,1) - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems1: do i = 1_pInt,prm%Nslip(f) - j = j+1_pInt - left_SlipSlip(j) = 1.0_pReal + prm%H_int(f) ! modified no system-dependent left part - right_SlipSlip(j) = abs(1.0_pReal-stt%s_slip(j,of) / (prm%tausat_slip(f)+ssat_offset)) **prm%a_slip & - * sign(1.0_pReal,1.0_pReal-stt%s_slip(j,of) / (prm%tausat_slip(f)+ssat_offset)) +! calculate left and right vectors + left_SlipSlip = 1.0_pReal + prm%H_int + xi_slip_sat_offset = prm%spr*sqrt(stt%sumF(of)) + right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & + * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) !-------------------------------------------------------------------------------------------------- -! Calculation of dot gamma - tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) - tau_slip_neg = tau_slip_pos - nonSchmidSystems: do k = 1,size(prm%nonSchmidCoeff) - tau_slip_pos = tau_slip_pos + prm%nonSchmidCoeff(k)* & - dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k, index_myFamily+i,ph)) - tau_slip_neg = tau_slip_neg +prm%nonSchmidCoeff(k)* & - dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) - enddo nonSchmidSystems - gdot_slip(j) = prm%gdot0_slip*0.5_pReal* & - ( (abs(tau_slip_pos)/(stt%s_slip(j,of)))**prm%n_slip*sign(1.0_pReal,tau_slip_pos) & - +(abs(tau_slip_neg)/(stt%s_slip(j,of)))**prm%n_slip*sign(1.0_pReal,tau_slip_neg)) - enddo slipSystems1 - enddo slipFamilies1 - - j = 0_pInt - twinFamilies1: do f = 1_pInt,size(prm%Ntwin,1) - index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems1: do i = 1_pInt,prm%Ntwin(f) - j = j+1_pInt +! shear rates + call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg) + dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) + dot%sumGamma(of) = sum(dot%gamma_slip(:,of)) + call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of)) + if (stt%sumF(of) < 0.98_pReal) dot%sumF(of) = sum(dot%gamma_twin(:,of)/prm%gamma_twin_char) !-------------------------------------------------------------------------------------------------- -! Calculation of dot vol frac - tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) - gdot_twin(j) = (1.0_pReal-stt%sumF(of))*& ! 1-F - prm%gdot0_twin*& - (abs(tau_twin)/stt%s_twin(j,of))**& - prm%n_twin*max(0.0_pReal,sign(1.0_pReal,tau_twin)) - enddo twinSystems1 - enddo twinFamilies1 +! hardening + hardeningSlip: do i = 1_pInt, prm%totalNslip + dot%xi_slip(i,of) = & + c_SlipSlip * left_SlipSlip(i) & + * dot_product(prm%interaction_SlipSlip(i,:),right_SlipSlip*dot%gamma_slip(:,of)) & + + & + dot_product(prm%interaction_SlipTwin(i,:),dot%gamma_twin(:,of)) + enddo hardeningSlip -!-------------------------------------------------------------------------------------------------- -! calculate the overall hardening based on above - do j = 1_pInt,prm%totalNslip - dst%s_slip(j,of) = c_SlipSlip * left_SlipSlip(j) * & ! evolution of slip resistance j - dot_product(prm%interaction_SlipSlip(j,1:prm%totalNslip),right_SlipSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor - dot_product(prm%interaction_SlipTwin(j,1:prm%totalNtwin),gdot_twin) ! dot gamma_twin modulated by right-side twin factor - enddo - dst%sumGamma(of) = dst%sumGamma(of) + sum(abs(gdot_slip)) - dst%accshear_slip(1:prm%totalNslip,of) = abs(gdot_slip) + hardeningTwin: do i = 1_pInt, prm%totalNtwin + dot%xi_twin(i,of) = & + c_TwinSlip & + * dot_product(prm%interaction_TwinSlip(i,:),dot%gamma_slip(:,of)) & + + & + c_TwinTwin & + * dot_product(prm%interaction_TwinTwin(i,:),dot%gamma_twin(:,of)) + enddo hardeningTwin - j = 0_pInt - twinFamilies2: do f = 1_pInt,size(prm%Ntwin,1) - index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems2: do i = 1_pInt,prm%Ntwin(f) - j = j+1_pInt - dst%s_twin(j,of) = & ! evolution of twin resistance j - c_TwinSlip * dot_product(prm%interaction_TwinSlip(j,1:prm%totalNslip),abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor - c_TwinTwin * dot_product(prm%interaction_TwinTwin(j,1:prm%totalNtwin),gdot_twin) ! dot gamma_twin modulated by right-side twin factor - if (stt%sumF(of) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0 - dst%sumF(of) = dst%sumF(of) + gdot_twin(j)/lattice_shearTwin(index_myFamily+i,ph) - dst%accshear_twin(j,of) = abs(gdot_twin(j)) - enddo twinSystems2 - enddo twinFamilies2 end associate + end subroutine plastic_phenopowerlaw_dotState +!-------------------------------------------------------------------------------------------------- +!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress +!> @details: Shear rates are calculated only optionally. NOTE: Agains the common convention, the +!> result (i.e. intent(out)) variables are the last to have the optional arguments at the end +!-------------------------------------------------------------------------------------------------- +subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, & + dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) + use prec, only: & + dNeq0 + use math, only: & + math_mul33xx33 + + implicit none + type(tParameters), intent(in) :: & + prm + type(tPhenopowerlawState), intent(in) :: & + stt + integer(pInt), intent(in) :: & + of + real(pReal), dimension(prm%totalNslip), intent(out) :: & + gdot_slip_pos, & + gdot_slip_neg + real(pReal), dimension(prm%totalNslip), optional, intent(out) :: & + dgdot_dtau_slip_pos, & + dgdot_dtau_slip_neg + real(pReal), dimension(3,3), intent(in) :: & + Mp + + real(pReal), dimension(prm%totalNslip) :: & + tau_slip_pos, & + tau_slip_neg + + integer(pInt) :: i, j + + do i = 1_pInt, prm%totalNslip + tau_slip_pos(i) = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i)) + tau_slip_neg(i) = tau_slip_pos(i) + do j = 1,size(prm%nonSchmidCoeff) + tau_slip_pos(i) = tau_slip_pos(i) + math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j,i)) + tau_slip_neg(i) = tau_slip_neg(i) + math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j,i)) + enddo + enddo + + gdot_slip_pos = 0.5_pReal*prm%gdot0_slip & + * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) + gdot_slip_neg = 0.5_pReal*prm%gdot0_slip & + * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) + + if (present(dgdot_dtau_slip_pos)) then + where(dNeq0(tau_slip_pos)) + dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_slip/tau_slip_pos + else where + dgdot_dtau_slip_pos = 0.0_pReal + end where + endif + if (present(dgdot_dtau_slip_neg)) then + where(dNeq0(tau_slip_neg)) + dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_slip/tau_slip_neg + else where + dgdot_dtau_slip_neg = 0.0_pReal + end where + endif + +end subroutine kinetics_slip + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress +!> @details: Shear rates are calculated only optionally. NOTE: Agains the common convention, the +!> result (i.e. intent(out)) variables are the last to have the optional arguments at the end +!-------------------------------------------------------------------------------------------------- +subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin) + use prec, only: & + dNeq0 + use math, only: & + math_mul33xx33 + + implicit none + type(tParameters), intent(in) :: & + prm + type(tPhenopowerlawState), intent(in) :: & + stt + integer(pInt), intent(in) :: & + of + real(pReal), dimension(3,3), intent(in) :: & + Mp + real(pReal), dimension(prm%totalNtwin), intent(out) :: & + gdot_twin + real(pReal), dimension(prm%totalNtwin), optional, intent(out) :: & + dgdot_dtau_twin + + real(pReal), dimension(prm%totalNtwin) :: & + tau_twin + integer(pInt) :: i + + do i = 1_pInt, prm%totalNtwin + tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) + enddo + gdot_twin = merge((1.0_pReal-stt%sumF(of))*prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin, & + 0.0_pReal, tau_twin>0.0_pReal) + + if (present(dgdot_dtau_twin)) then + where(dNeq0(tau_twin)) + dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin + else where + dgdot_dtau_twin = 0.0_pReal + end where + endif + +end subroutine kinetics_twin + + !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- -function plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) +function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) use material, only: & material_phase, & plasticState, & phasememberAt, & phase_plasticityInstance - use lattice, only: & - lattice_Sslip_v, & - lattice_Stwin_v, & - lattice_NslipSystem, & - lattice_NtwinSystem, & - lattice_NnonSchmid + use math, only: & + math_mul33xx33, & + math_Mandel6to33 implicit none - real(pReal), dimension(6), intent(in) :: & - Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress integer(pInt), intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element !< microstructure state + instance, & + of - real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults) :: & - plastic_phenopowerlaw_postResults + real(pReal), dimension(sum(plastic_phenopowerlaw_sizePostResult(:,instance))) :: & + postResults integer(pInt) :: & - ph, of, & - o,f,i,c,j,k, & - index_myFamily + o,c,i,j real(pReal) :: & - tau_slip_pos,tau_slip_neg,tau + tau_slip_pos, tau_slip_neg + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_slip_pos,gdot_slip_neg type(tParameters) :: prm - type(tPhenopowerlawState) :: stt, dst + type(tPhenopowerlawState) :: stt - of = phasememberAt(ipc,ip,el) - ph = material_phase(ipc,ip,el) + associate( prm => param(instance), stt => state(instance)) - associate( prm => param(phase_plasticityInstance(ph)), & - stt => state(phase_plasticityInstance(ph)), & - dst => dotState(phase_plasticityInstance(ph))) - - plastic_phenopowerlaw_postResults = 0.0_pReal + postResults = 0.0_pReal c = 0_pInt outputsLoop: do o = 1_pInt,size(prm%outputID) select case(prm%outputID(o)) + case (resistance_slip_ID) - plastic_phenopowerlaw_postResults(c+1_pInt:c+prm%totalNslip) = stt%s_slip(1:prm%totalNslip,of) + postResults(c+1_pInt:c+prm%totalNslip) = stt%xi_slip(1:prm%totalNslip,of) c = c + prm%totalNslip - case (accumulatedshear_slip_ID) - plastic_phenopowerlaw_postResults(c+1_pInt:c+prm%totalNslip) = stt%accshear_slip(1:prm%totalNslip,of) + postResults(c+1_pInt:c+prm%totalNslip) = stt%gamma_slip(1:prm%totalNslip,of) c = c + prm%totalNslip - case (shearrate_slip_ID) - j = 0_pInt - slipFamilies1: do f = 1_pInt,size(prm%Nslip,1) - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems1: do i = 1_pInt,prm%Nslip(f) - j = j + 1_pInt - tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) - tau_slip_neg = tau_slip_pos - do k = 1,lattice_NnonSchmid(ph) - tau_slip_pos = tau_slip_pos +prm%nonSchmidCoeff(k)* & - dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) - tau_slip_neg = tau_slip_neg +prm%nonSchmidCoeff(k)* & - dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) - enddo - plastic_phenopowerlaw_postResults(c+j) = prm%gdot0_slip*0.5_pReal* & - ((abs(tau_slip_pos)/stt%s_slip(j,of))**prm%n_slip & - *sign(1.0_pReal,tau_slip_pos) & - +(abs(tau_slip_neg)/(stt%s_slip(j,of)))**prm%n_slip & - *sign(1.0_pReal,tau_slip_neg)) - enddo slipSystems1 - enddo slipFamilies1 + call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg) + postResults(c+1_pInt:c+prm%totalNslip) = gdot_slip_pos+gdot_slip_neg c = c + prm%totalNslip - case (resolvedstress_slip_ID) - j = 0_pInt - slipFamilies2: do f = 1_pInt,size(prm%Nslip,1) - index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family - slipSystems2: do i = 1_pInt,prm%Nslip(f) - j = j + 1_pInt - plastic_phenopowerlaw_postResults(c+j) = & - dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) - enddo slipSystems2 - enddo slipFamilies2 + do i = 1_pInt, prm%totalNslip + tau_slip_pos = math_mul33xx33(Mp,prm%Schmid_slip(1:3,1:3,i)) + tau_slip_neg = tau_slip_pos + !do j = 1,size(prm%nonSchmidCoeff) + ! tau_slip_pos = tau_slip_pos + math_mul33xx33(S,prm%nonSchmid_pos(1:3,1:3,j,i)) + ! tau_slip_neg = tau_slip_neg + math_mul33xx33(S,prm%nonSchmid_neg(1:3,1:3,j,i)) + !enddo + postResults(c+i) = 0.5_pReal*(tau_slip_pos+tau_slip_neg) + enddo c = c + prm%totalNslip - case (totalshear_ID) - plastic_phenopowerlaw_postResults(c+1_pInt) = stt%sumGamma(of) - c = c + 1_pInt - case (resistance_twin_ID) - plastic_phenopowerlaw_postResults(c+1_pInt:c+prm%totalNtwin) = & - stt%s_twin(1:prm%totalNtwin,of) + postResults(c+1_pInt:c+prm%totalNtwin) = stt%xi_twin(1:prm%totalNtwin,of) c = c + prm%totalNtwin - case (accumulatedshear_twin_ID) - plastic_phenopowerlaw_postResults(c+1_pInt:c+prm%totalNtwin) = & - stt%accshear_twin(1:prm%totalNtwin,of) + postResults(c+1_pInt:c+prm%totalNtwin) = stt%gamma_twin(1:prm%totalNtwin,of) c = c + prm%totalNtwin - case (shearrate_twin_ID) - j = 0_pInt - twinFamilies1: do f = 1_pInt,size(prm%Ntwin,1) - index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems1: do i = 1_pInt,prm%Ntwin(f) - j = j + 1_pInt - tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) - plastic_phenopowerlaw_postResults(c+j) = (1.0_pReal-stt%sumF(of))*& ! 1-F - prm%gdot0_twin*& - (abs(tau)/stt%s_twin(j,of))**& - prm%n_twin*max(0.0_pReal,sign(1.0_pReal,tau)) - enddo twinSystems1 - enddo twinFamilies1 + call kinetics_twin(prm,stt,of,Mp,postResults(c+1_pInt:c+prm%totalNtwin)) c = c + prm%totalNtwin - case (resolvedstress_twin_ID) - j = 0_pInt - twinFamilies2: do f = 1_pInt,size(prm%Ntwin,1) - index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family - twinSystems2: do i = 1_pInt,prm%Ntwin(f) - j = j + 1_pInt - plastic_phenopowerlaw_postResults(c+j) = & - dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) - enddo twinSystems2 - enddo twinFamilies2 + do i = 1_pInt, prm%totalNtwin + postResults(c+i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) + enddo c = c + prm%totalNtwin + case (totalshear_ID) + postResults(c+1_pInt) = stt%sumGamma(of) + c = c + 1_pInt case (totalvolfrac_twin_ID) - plastic_phenopowerlaw_postResults(c+1_pInt) = stt%sumF(of) + postResults(c+1_pInt) = stt%sumF(of) c = c + 1_pInt end select enddo outputsLoop end associate + end function plastic_phenopowerlaw_postResults end module plastic_phenopowerlaw