From 79663a7f7624aa7fffcc6e8346a592e06ae1785a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Feb 2012 16:00:00 +0000 Subject: [PATCH] polishing: adding pInt, removing unused use-statements etc --- code/constitutive.f90 | 85 +++--- code/constitutive_dislotwin.f90 | 398 +++++++++++++------------ code/constitutive_j2.f90 | 77 +++-- code/constitutive_phenopowerlaw.f90 | 114 +++---- code/constitutive_titanmod.f90 | 443 +++++++++++++--------------- 5 files changed, 544 insertions(+), 573 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 2d6135b46..1abb04377 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -68,8 +68,9 @@ CONTAINS !* Module initialization * !************************************** subroutine constitutive_init() -use prec, only: pReal,pInt -use debug, only: debug_verbosity, debug_selectiveDebugger, debug_e, debug_i, debug_g +use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) +use prec, only: pInt +use debug, only: debug_verbosity use numerics, only: numerics_integrator use IO, only: IO_error, IO_open_file, IO_open_jobFile_stat, IO_write_jobFile use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips @@ -82,7 +83,7 @@ use constitutive_nonlocal implicit none -integer(pInt), parameter :: fileunit = 200 +integer(pInt), parameter :: fileunit = 200_pInt integer(pInt) g, & ! grain number i, & ! integration point number e, & ! element number @@ -114,7 +115,7 @@ close(fileunit) ! --- WRITE DESCRIPTION FILE FOR CONSTITUTIVE PHASE OUTPUT --- call IO_write_jobFile(fileunit,'outputConstitutive') -do p = 1,material_Nphase +do p = 1_pInt,material_Nphase i = phase_constitutionInstance(p) ! which instance of a constitution is present phase knownConstitution = .true. ! assume valid select case(phase_constitution(p)) ! split per constitiution @@ -141,7 +142,7 @@ do p = 1,material_Nphase write(fileunit,*) if (knownConstitution) then write(fileunit,'(a)') '(constitution)'//char(9)//trim(phase_constitution(p)) - do e = 1,phase_Noutput(p) + do e = 1_pInt,phase_Noutput(p) write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) enddo endif @@ -166,22 +167,22 @@ allocate(constitutive_aTolState(gMax,iMax,eMax)) allocate(constitutive_sizeDotState(gMax,iMax,eMax)) ; constitutive_sizeDotState = 0_pInt allocate(constitutive_sizeState(gMax,iMax,eMax)) ; constitutive_sizeState = 0_pInt allocate(constitutive_sizePostResults(gMax,iMax,eMax)); constitutive_sizePostResults = 0_pInt -if (any(numerics_integrator == 1)) then +if (any(numerics_integrator == 1_pInt)) then allocate(constitutive_previousDotState(gMax,iMax,eMax)) allocate(constitutive_previousDotState2(gMax,iMax,eMax)) endif -if (any(numerics_integrator == 4)) then +if (any(numerics_integrator == 4_pInt)) then allocate(constitutive_RK4dotState(gMax,iMax,eMax)) endif -if (any(numerics_integrator == 5)) then +if (any(numerics_integrator == 5_pInt)) then allocate(constitutive_RKCK45dotState(6,gMax,iMax,eMax)) endif !$OMP PARALLEL DO PRIVATE(myNgrains,myInstance) - do e = 1,mesh_NcpElems ! loop over elements + do e = 1_pInt,mesh_NcpElems ! loop over elements myNgrains = homogenization_Ngrains(mesh_element(3,e)) - do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs - do g = 1,myNgrains ! loop over grains + do i = 1_pInt,FE_Nips(mesh_element(2,e)) ! loop over IPs + do g = 1_pInt,myNgrains ! loop over grains myInstance = phase_constitutionInstance(material_phase(g,i,e)) select case(phase_constitution(material_phase(g,i,e))) @@ -194,15 +195,15 @@ endif allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then + if (any(numerics_integrator == 1_pInt)) then allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 4)) then + if (any(numerics_integrator == 4_pInt)) then allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 5)) then - do s = 1,6 + if (any(numerics_integrator == 5_pInt)) then + do s = 1_pInt,6_pInt allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) enddo endif @@ -221,15 +222,15 @@ endif allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then + if (any(numerics_integrator == 1_pInt)) then allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 4)) then + if (any(numerics_integrator == 4_pInt)) then allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 5)) then - do s = 1,6 + if (any(numerics_integrator == 5_pInt)) then + do s = 1_pInt,6_pInt allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) enddo endif @@ -248,15 +249,15 @@ endif allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then + if (any(numerics_integrator == 1_pInt)) then allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 4)) then + if (any(numerics_integrator == 4_pInt)) then allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 5)) then - do s = 1,6 + if (any(numerics_integrator == 5_pInt)) then + do s = 1_pInt,6_pInt allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) enddo endif @@ -275,15 +276,15 @@ endif allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then + if (any(numerics_integrator == 1_pInt)) then allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 4)) then + if (any(numerics_integrator == 4_pInt)) then allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 5)) then - do s = 1,6 + if (any(numerics_integrator == 5_pInt)) then + do s = 1_pInt,6_pInt allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) enddo endif @@ -302,15 +303,15 @@ endif allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) - if (any(numerics_integrator == 1)) then + if (any(numerics_integrator == 1_pInt)) then allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 4)) then + if (any(numerics_integrator == 4_pInt)) then allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) endif - if (any(numerics_integrator == 5)) then - do s = 1,6 + if (any(numerics_integrator == 5_pInt)) then + do s = 1_pInt,6_pInt allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) enddo endif @@ -340,7 +341,7 @@ constitutive_maxSizePostResults = maxval(constitutive_sizePostResults) write(6,*) '<<<+- constitutive init -+>>>' write(6,*) '$Id$' #include "compilation_info.f90" - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then write(6,'(a32,1x,7(i8,1x))') 'constitutive_state0: ', shape(constitutive_state0) write(6,'(a32,1x,7(i8,1x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0) write(6,'(a32,1x,7(i8,1x))') 'constitutive_subState0: ', shape(constitutive_subState0) @@ -457,11 +458,7 @@ endfunction subroutine constitutive_microstructure(Temperature, Fe, Fp, ipc, ip, el) use prec, only: pReal,pInt use material, only: phase_constitution, & - material_phase, & - homogenization_maxNgrains -use mesh, only: mesh_NcpElems, & - mesh_maxNips, & - mesh_maxNipNeighbors + material_phase use constitutive_j2, only: constitutive_j2_label, & constitutive_j2_microstructure use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, & @@ -581,8 +578,7 @@ use debug, only: debug_cumDotStateCalls, & debug_cumDotStateTicks, & debug_verbosity use mesh, only: mesh_NcpElems, & - mesh_maxNips, & - mesh_maxNipNeighbors + mesh_maxNips use material, only: phase_constitution, & material_phase, & homogenization_maxNgrains @@ -620,7 +616,7 @@ integer(pLongInt) tick, tock, & tickrate, & maxticks -if (debug_verbosity > 0) then +if (debug_verbosity > 0_pInt) then call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) endif @@ -644,7 +640,7 @@ select case (phase_constitution(material_phase(ipc,ip,el))) end select -if (debug_verbosity > 6) then +if (debug_verbosity > 6_pInt) then call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) !$OMP CRITICAL (debugTimingDotState) debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt @@ -699,7 +695,7 @@ integer(pLongInt) tick, tock, & maxticks -if (debug_verbosity > 0) then +if (debug_verbosity > 0_pInt) then call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) endif @@ -722,7 +718,7 @@ select case (phase_constitution(material_phase(ipc,ip,el))) end select -if (debug_verbosity > 6) then +if (debug_verbosity > 6_pInt) then call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) !$OMP CRITICAL (debugTimingDotTemperature) debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt @@ -748,8 +744,7 @@ function constitutive_postResults(Tstar_v, Fe, Temperature, dt, ipc, ip, el) !********************************************************************* use prec, only: pReal,pInt use mesh, only: mesh_NcpElems, & - mesh_maxNips, & - mesh_maxNipNeighbors + mesh_maxNips use material, only: phase_constitution, & material_phase, & homogenization_maxNgrains diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index a0d459950..b011f0855 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -133,7 +133,7 @@ subroutine constitutive_dislotwin_init(file) !************************************** !* Module initialization * !************************************** -use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) +use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pInt,pReal use math, only: math_Mandel3333to66,math_Voigt66to3333,math_mul3x3 use mesh, only: mesh_maxNips, mesh_NcpElems @@ -144,9 +144,10 @@ use lattice !* Input variables integer(pInt), intent(in) :: file !* Local variables -integer(pInt), parameter :: maxNchunks = 21 +integer(pInt), parameter :: maxNchunks = 21_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions -integer(pInt) section,maxNinstance,f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t1,t2,ns,nt,mySize,myStructure,maxTotalNslip,maxTotalNtwin +integer(pInt) :: section, maxNinstance,mySize,myStructure,maxTotalNslip,maxTotalNtwin,& + f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t1,t2,ns,nt character(len=64) tag character(len=1024) line @@ -157,130 +158,130 @@ character(len=1024) line #include "compilation_info.f90" !$OMP END CRITICAL (write2out) -maxNinstance = count(phase_constitution == constitutive_dislotwin_label) -if (maxNinstance == 0) return +maxNinstance = int(count(phase_constitution == constitutive_dislotwin_label),pInt) +if (maxNinstance == 0_pInt) return !* Space allocation for global variables allocate(constitutive_dislotwin_sizeDotState(maxNinstance)) + constitutive_dislotwin_sizeDotState = 0_pInt allocate(constitutive_dislotwin_sizeState(maxNinstance)) + constitutive_dislotwin_sizeState = 0_pInt allocate(constitutive_dislotwin_sizePostResults(maxNinstance)) + constitutive_dislotwin_sizePostResults = 0_pInt allocate(constitutive_dislotwin_sizePostResult(maxval(phase_Noutput),maxNinstance)) + constitutive_dislotwin_sizePostResult = 0_pInt allocate(constitutive_dislotwin_output(maxval(phase_Noutput),maxNinstance)) + constitutive_dislotwin_output = '' allocate(constitutive_dislotwin_Noutput(maxNinstance)) -constitutive_dislotwin_sizeDotState = 0_pInt -constitutive_dislotwin_sizeState = 0_pInt -constitutive_dislotwin_sizePostResults = 0_pInt -constitutive_dislotwin_sizePostResult = 0_pInt -constitutive_dislotwin_output = '' -constitutive_dislotwin_Noutput = 0_pInt + constitutive_dislotwin_Noutput = 0_pInt allocate(constitutive_dislotwin_structureName(maxNinstance)) + constitutive_dislotwin_structureName = '' allocate(constitutive_dislotwin_structure(maxNinstance)) + constitutive_dislotwin_structure = 0_pInt allocate(constitutive_dislotwin_Nslip(lattice_maxNslipFamily,maxNinstance)) + constitutive_dislotwin_Nslip = 0_pInt allocate(constitutive_dislotwin_Ntwin(lattice_maxNtwinFamily,maxNinstance)) + constitutive_dislotwin_Ntwin = 0_pInt allocate(constitutive_dislotwin_slipFamily(lattice_maxNslip,maxNinstance)) + constitutive_dislotwin_slipFamily = 0_pInt allocate(constitutive_dislotwin_twinFamily(lattice_maxNtwin,maxNinstance)) + constitutive_dislotwin_twinFamily = 0_pInt allocate(constitutive_dislotwin_slipSystemLattice(lattice_maxNslip,maxNinstance)) + constitutive_dislotwin_slipSystemLattice = 0_pInt allocate(constitutive_dislotwin_twinSystemLattice(lattice_maxNtwin,maxNinstance)) + constitutive_dislotwin_twinSystemLattice = 0_pInt allocate(constitutive_dislotwin_totalNslip(maxNinstance)) + constitutive_dislotwin_totalNslip = 0_pInt allocate(constitutive_dislotwin_totalNtwin(maxNinstance)) -constitutive_dislotwin_structureName = '' -constitutive_dislotwin_structure = 0_pInt -constitutive_dislotwin_Nslip = 0_pInt -constitutive_dislotwin_Ntwin = 0_pInt -constitutive_dislotwin_slipFamily = 0_pInt -constitutive_dislotwin_twinFamily = 0_pInt -constitutive_dislotwin_slipSystemLattice = 0_pInt -constitutive_dislotwin_twinSystemLattice = 0_pInt -constitutive_dislotwin_totalNslip = 0_pInt -constitutive_dislotwin_totalNtwin = 0_pInt + constitutive_dislotwin_totalNtwin = 0_pInt allocate(constitutive_dislotwin_CoverA(maxNinstance)) + constitutive_dislotwin_CoverA = 0.0_pReal allocate(constitutive_dislotwin_C11(maxNinstance)) + constitutive_dislotwin_C11 = 0.0_pReal allocate(constitutive_dislotwin_C12(maxNinstance)) + constitutive_dislotwin_C12 = 0.0_pReal allocate(constitutive_dislotwin_C13(maxNinstance)) + constitutive_dislotwin_C13 = 0.0_pReal allocate(constitutive_dislotwin_C33(maxNinstance)) + constitutive_dislotwin_C33 = 0.0_pReal allocate(constitutive_dislotwin_C44(maxNinstance)) + constitutive_dislotwin_C44 = 0.0_pReal allocate(constitutive_dislotwin_Gmod(maxNinstance)) + constitutive_dislotwin_Gmod = 0.0_pReal allocate(constitutive_dislotwin_CAtomicVolume(maxNinstance)) + constitutive_dislotwin_CAtomicVolume = 0.0_pReal allocate(constitutive_dislotwin_D0(maxNinstance)) + constitutive_dislotwin_D0 = 0.0_pReal allocate(constitutive_dislotwin_Qsd(maxNinstance)) + constitutive_dislotwin_Qsd = 0.0_pReal allocate(constitutive_dislotwin_GrainSize(maxNinstance)) + constitutive_dislotwin_GrainSize = 0.0_pReal allocate(constitutive_dislotwin_p(maxNinstance)) + constitutive_dislotwin_p = 0.0_pReal allocate(constitutive_dislotwin_q(maxNinstance)) + constitutive_dislotwin_q = 0.0_pReal allocate(constitutive_dislotwin_MaxTwinFraction(maxNinstance)) + constitutive_dislotwin_MaxTwinFraction = 0.0_pReal allocate(constitutive_dislotwin_r(maxNinstance)) + constitutive_dislotwin_r = 0.0_pReal allocate(constitutive_dislotwin_CEdgeDipMinDistance(maxNinstance)) + constitutive_dislotwin_CEdgeDipMinDistance = 0.0_pReal allocate(constitutive_dislotwin_Cmfptwin(maxNinstance)) + constitutive_dislotwin_Cmfptwin = 0.0_pReal allocate(constitutive_dislotwin_Cthresholdtwin(maxNinstance)) + constitutive_dislotwin_Cthresholdtwin = 0.0_pReal allocate(constitutive_dislotwin_SolidSolutionStrength(maxNinstance)) + constitutive_dislotwin_SolidSolutionStrength = 0.0_pReal allocate(constitutive_dislotwin_L0(maxNinstance)) + constitutive_dislotwin_L0 = 0.0_pReal allocate(constitutive_dislotwin_aTolRho(maxNinstance)) + constitutive_dislotwin_aTolRho = 0.0_pReal allocate(constitutive_dislotwin_Cslip_66(6,6,maxNinstance)) + constitutive_dislotwin_Cslip_66 = 0.0_pReal allocate(constitutive_dislotwin_Cslip_3333(3,3,3,3,maxNinstance)) + constitutive_dislotwin_Cslip_3333 = 0.0_pReal allocate(constitutive_dislotwin_sbResistance(maxNinstance)) + constitutive_dislotwin_sbResistance = 0.0_pReal allocate(constitutive_dislotwin_sbVelocity(maxNinstance)) + constitutive_dislotwin_sbVelocity = 0.0_pReal allocate(constitutive_dislotwin_SFE_0K(maxNinstance)) + constitutive_dislotwin_SFE_0K = 0.0_pReal allocate(constitutive_dislotwin_dSFE_dT(maxNinstance)) -constitutive_dislotwin_CoverA = 0.0_pReal -constitutive_dislotwin_C11 = 0.0_pReal -constitutive_dislotwin_C12 = 0.0_pReal -constitutive_dislotwin_C13 = 0.0_pReal -constitutive_dislotwin_C33 = 0.0_pReal -constitutive_dislotwin_C44 = 0.0_pReal -constitutive_dislotwin_Gmod = 0.0_pReal -constitutive_dislotwin_CAtomicVolume = 0.0_pReal -constitutive_dislotwin_D0 = 0.0_pReal -constitutive_dislotwin_Qsd = 0.0_pReal -constitutive_dislotwin_GrainSize = 0.0_pReal -constitutive_dislotwin_p = 0.0_pReal -constitutive_dislotwin_q = 0.0_pReal -constitutive_dislotwin_MaxTwinFraction = 0.0_pReal -constitutive_dislotwin_r = 0.0_pReal -constitutive_dislotwin_CEdgeDipMinDistance = 0.0_pReal -constitutive_dislotwin_Cmfptwin = 0.0_pReal -constitutive_dislotwin_Cthresholdtwin = 0.0_pReal -constitutive_dislotwin_SolidSolutionStrength= 0.0_pReal -constitutive_dislotwin_L0 = 0.0_pReal -constitutive_dislotwin_aTolRho = 0.0_pReal -constitutive_dislotwin_Cslip_66 = 0.0_pReal -constitutive_dislotwin_Cslip_3333 = 0.0_pReal -constitutive_dislotwin_sbResistance = 0.0_pReal -constitutive_dislotwin_sbVelocity = 0.0_pReal -constitutive_dislotwin_SFE_0K = 0.0_pReal -constitutive_dislotwin_dSFE_dT = 0.0_pReal + constitutive_dislotwin_dSFE_dT = 0.0_pReal allocate(constitutive_dislotwin_rhoEdge0(lattice_maxNslipFamily,maxNinstance)) + constitutive_dislotwin_rhoEdge0 = 0.0_pReal allocate(constitutive_dislotwin_rhoEdgeDip0(lattice_maxNslipFamily,maxNinstance)) + constitutive_dislotwin_rhoEdgeDip0 = 0.0_pReal allocate(constitutive_dislotwin_burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_dislotwin_burgersPerSlipFamily = 0.0_pReal allocate(constitutive_dislotwin_burgersPerTwinFamily(lattice_maxNtwinFamily,maxNinstance)) + constitutive_dislotwin_burgersPerTwinFamily = 0.0_pReal allocate(constitutive_dislotwin_QedgePerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_dislotwin_QedgePerSlipFamily = 0.0_pReal allocate(constitutive_dislotwin_v0PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_dislotwin_v0PerSlipFamily = 0.0_pReal allocate(constitutive_dislotwin_Ndot0PerTwinFamily(lattice_maxNtwinFamily,maxNinstance)) + constitutive_dislotwin_Ndot0PerTwinFamily = 0.0_pReal allocate(constitutive_dislotwin_twinsizePerTwinFamily(lattice_maxNtwinFamily,maxNinstance)) + constitutive_dislotwin_twinsizePerTwinFamily = 0.0_pReal allocate(constitutive_dislotwin_CLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance)) -constitutive_dislotwin_rhoEdge0 = 0.0_pReal -constitutive_dislotwin_rhoEdgeDip0 = 0.0_pReal -constitutive_dislotwin_burgersPerSlipFamily = 0.0_pReal -constitutive_dislotwin_burgersPerTwinFamily = 0.0_pReal -constitutive_dislotwin_QedgePerSlipFamily = 0.0_pReal -constitutive_dislotwin_v0PerSlipFamily = 0.0_pReal -constitutive_dislotwin_Ndot0PerTwinFamily = 0.0_pReal -constitutive_dislotwin_twinsizePerTwinFamily = 0.0_pReal -constitutive_dislotwin_CLambdaSlipPerSlipFamily = 0.0_pReal + constitutive_dislotwin_CLambdaSlipPerSlipFamily = 0.0_pReal allocate(constitutive_dislotwin_interactionSlipSlip(lattice_maxNinteraction,maxNinstance)) + constitutive_dislotwin_interactionSlipSlip = 0.0_pReal allocate(constitutive_dislotwin_interactionSlipTwin(lattice_maxNinteraction,maxNinstance)) + constitutive_dislotwin_interactionSlipTwin = 0.0_pReal allocate(constitutive_dislotwin_interactionTwinSlip(lattice_maxNinteraction,maxNinstance)) + constitutive_dislotwin_interactionTwinSlip = 0.0_pReal allocate(constitutive_dislotwin_interactionTwinTwin(lattice_maxNinteraction,maxNinstance)) -constitutive_dislotwin_interactionSlipSlip = 0.0_pReal -constitutive_dislotwin_interactionSlipTwin = 0.0_pReal -constitutive_dislotwin_interactionTwinSlip = 0.0_pReal -constitutive_dislotwin_interactionTwinTwin = 0.0_pReal + constitutive_dislotwin_interactionTwinTwin = 0.0_pReal allocate(constitutive_dislotwin_sbSv(6,6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) -constitutive_dislotwin_sbSv = 0.0_pReal + constitutive_dislotwin_sbSv = 0.0_pReal !* Readout data from material.config file rewind(file) line = '' -section = 0 +section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to read(file,'(a1024)',END=100) line @@ -297,7 +298,7 @@ do ! read thru sections of if (section > 0_pInt .and. phase_constitution(section) == constitutive_dislotwin_label) then ! one of my sections i = phase_constitutionInstance(section) ! which instance of my constitution is present phase positions = IO_stringPos(line,maxNchunks) - tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key select case(tag) case ('constitution') cycle @@ -414,7 +415,7 @@ enddo if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt,e=i) if (sum(constitutive_dislotwin_Nslip(:,i)) <= 0_pInt) call IO_error(241_pInt,e=i,ext_msg='nslip') if (sum(constitutive_dislotwin_Ntwin(:,i)) < 0_pInt) call IO_error(241_pInt,e=i,ext_msg='ntwin') - do f = 1,lattice_maxNslipFamily + do f = 1_pInt,lattice_maxNslipFamily if (constitutive_dislotwin_Nslip(f,i) > 0_pInt) then if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='rhoEdge0') if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='rhoEdgeDip0') @@ -422,7 +423,7 @@ enddo if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='v0') endif enddo - do f = 1,lattice_maxNtwinFamily + do f = 1_pInt,lattice_maxNtwinFamily if (constitutive_dislotwin_Ntwin(f,i) > 0_pInt) then if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='twinburgers') if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(241_pInt,e=i,ext_msg='ndot0') @@ -450,35 +451,35 @@ maxTotalNslip = maxval(constitutive_dislotwin_totalNslip) maxTotalNtwin = maxval(constitutive_dislotwin_totalNtwin) allocate(constitutive_dislotwin_burgersPerSlipSystem(maxTotalNslip, maxNinstance)) + constitutive_dislotwin_burgersPerSlipSystem = 0.0_pReal allocate(constitutive_dislotwin_burgersPerTwinSystem(maxTotalNtwin, maxNinstance)) + constitutive_dislotwin_burgersPerTwinSystem= 0.0_pReal allocate(constitutive_dislotwin_QedgePerSlipSystem(maxTotalNslip, maxNinstance)) + constitutive_dislotwin_QedgePerSlipSystem = 0.0_pReal allocate(constitutive_dislotwin_v0PerSlipSystem(maxTotalNslip, maxNinstance)) + constitutive_dislotwin_v0PerSlipSystem = 0.0_pReal allocate(constitutive_dislotwin_Ndot0PerTwinSystem(maxTotalNtwin, maxNinstance)) + constitutive_dislotwin_Ndot0PerTwinSystem = 0.0_pReal allocate(constitutive_dislotwin_twinsizePerTwinSystem(maxTotalNtwin, maxNinstance)) + constitutive_dislotwin_twinsizePerTwinSystem = 0.0_pReal allocate(constitutive_dislotwin_CLambdaSlipPerSlipSystem(maxTotalNslip, maxNinstance)) -constitutive_dislotwin_burgersPerSlipSystem = 0.0_pReal -constitutive_dislotwin_burgersPerTwinSystem = 0.0_pReal -constitutive_dislotwin_QedgePerSlipSystem = 0.0_pReal -constitutive_dislotwin_v0PerSlipSystem = 0.0_pReal -constitutive_dislotwin_Ndot0PerTwinSystem = 0.0_pReal -constitutive_dislotwin_twinsizePerTwinSystem = 0.0_pReal -constitutive_dislotwin_CLambdaSlipPerSlipSystem = 0.0_pReal + constitutive_dislotwin_CLambdaSlipPerSlipSystem = 0.0_pReal allocate(constitutive_dislotwin_interactionMatrixSlipSlip(maxTotalNslip,maxTotalNslip,maxNinstance)) + constitutive_dislotwin_interactionMatrixSlipSlip = 0.0_pReal allocate(constitutive_dislotwin_interactionMatrixSlipTwin(maxTotalNslip,maxTotalNtwin,maxNinstance)) + constitutive_dislotwin_interactionMatrixSlipTwin = 0.0_pReal allocate(constitutive_dislotwin_interactionMatrixTwinSlip(maxTotalNtwin,maxTotalNslip,maxNinstance)) + constitutive_dislotwin_interactionMatrixTwinSlip = 0.0_pReal allocate(constitutive_dislotwin_interactionMatrixTwinTwin(maxTotalNtwin,maxTotalNtwin,maxNinstance)) + constitutive_dislotwin_interactionMatrixTwinTwin = 0.0_pReal allocate(constitutive_dislotwin_forestProjectionEdge(maxTotalNslip,maxTotalNslip,maxNinstance)) -constitutive_dislotwin_interactionMatrixSlipSlip = 0.0_pReal -constitutive_dislotwin_interactionMatrixSlipTwin = 0.0_pReal -constitutive_dislotwin_interactionMatrixTwinSlip = 0.0_pReal -constitutive_dislotwin_interactionMatrixTwinTwin = 0.0_pReal -constitutive_dislotwin_forestProjectionEdge = 0.0_pReal + constitutive_dislotwin_forestProjectionEdge = 0.0_pReal allocate(constitutive_dislotwin_Ctwin_66(6,6,maxTotalNtwin,maxNinstance)) + constitutive_dislotwin_Ctwin_66 = 0.0_pReal allocate(constitutive_dislotwin_Ctwin_3333(3,3,3,3,maxTotalNtwin,maxNinstance)) -constitutive_dislotwin_Ctwin_66 = 0.0_pReal -constitutive_dislotwin_Ctwin_3333 = 0.0_pReal + constitutive_dislotwin_Ctwin_3333 = 0.0_pReal do i = 1_pInt,maxNinstance myStructure = constitutive_dislotwin_structure(i) @@ -504,11 +505,11 @@ do i = 1_pInt,maxNinstance !* Determine size of state array ns = constitutive_dislotwin_totalNslip(i) nt = constitutive_dislotwin_totalNtwin(i) - constitutive_dislotwin_sizeDotState(i) = & - size(constitutive_dislotwin_listBasicSlipStates)*ns+size(constitutive_dislotwin_listBasicTwinStates)*nt - constitutive_dislotwin_sizeState(i) = & - constitutive_dislotwin_sizeDotState(i)+ & - size(constitutive_dislotwin_listDependentSlipStates)*ns+size(constitutive_dislotwin_listDependentTwinStates)*nt + constitutive_dislotwin_sizeDotState(i) =int(size(constitutive_dislotwin_listBasicSlipStates),pInt)*ns& + +int(size(constitutive_dislotwin_listBasicTwinStates),pInt)*nt + constitutive_dislotwin_sizeState(i) = constitutive_dislotwin_sizeDotState(i)& + + int(size(constitutive_dislotwin_listDependentSlipStates),pInt)*ns& + + int(size(constitutive_dislotwin_listDependentTwinStates),pInt)*nt !* Determine size of postResults array do o = 1_pInt,constitutive_dislotwin_Noutput(i) @@ -665,7 +666,7 @@ function constitutive_dislotwin_stateInit(myInstance) !********************************************************************* use prec, only: pReal,pInt use math, only: pi -use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily +use lattice, only: lattice_maxNslipFamily implicit none !* Input-Output variables @@ -701,7 +702,7 @@ constitutive_dislotwin_stateInit(ns+1:2_pInt*ns) = rhoEdgeDip0 forall (s = 1_pInt:ns) & invLambdaSlip0(s) = sqrt(dot_product((rhoEdge0+rhoEdgeDip0),constitutive_dislotwin_forestProjectionEdge(1:ns,s,myInstance)))/ & constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,myInstance) -constitutive_dislotwin_stateInit(2_pInt*ns+nt+1:3_pInt*ns+nt) = invLambdaSlip0 +constitutive_dislotwin_stateInit(2_pInt*ns+nt+1_pInt:3_pInt*ns+nt) = invLambdaSlip0 forall (s = 1_pInt:ns) & MeanFreePathSlip0(s) = & @@ -715,14 +716,14 @@ sqrt(dot_product((rhoEdge0+rhoEdgeDip0),constitutive_dislotwin_interactionMatrix constitutive_dislotwin_stateInit(5_pInt*ns+3_pInt*nt+1:6_pInt*ns+3_pInt*nt) = tauSlipThreshold0 !* Initialize dependent twin microstructural variables -forall (t = 1:nt) & +forall (t = 1_pInt:nt) & MeanFreePathTwin0(t) = constitutive_dislotwin_GrainSize(myInstance) -constitutive_dislotwin_stateInit(5*ns+2*nt+1:5*ns+3*nt) = MeanFreePathTwin0 +constitutive_dislotwin_stateInit(5_pInt*ns+2_pInt*nt+1_pInt:5_pInt*ns+3_pInt*nt) = MeanFreePathTwin0 -forall (t = 1:nt) & +forall (t = 1_pInt:nt) & TwinVolume0(t) = & (pi/6.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)*MeanFreePathTwin0(t)**(2.0_pReal) -constitutive_dislotwin_stateInit(6*ns+4*nt+1:6*ns+5*nt) = TwinVolume0 +constitutive_dislotwin_stateInit(6_pInt*ns+4_pInt*nt+1_pInt:6_pInt*ns+5_pInt*nt) = TwinVolume0 !write(6,*) '#STATEINIT#' !write(6,*) @@ -782,13 +783,13 @@ ns = constitutive_dislotwin_totalNslip(myInstance) nt = constitutive_dislotwin_totalNtwin(myInstance) !* Total twin volume fraction -sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 +sumf = sum(state(g,ip,el)%p((2_pInt*ns+1_pInt):(2_pInt*ns+nt))) ! safe for nt == 0 !* Homogenized elasticity matrix constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*constitutive_dislotwin_Cslip_66(:,:,myInstance) -do i=1,nt +do i=1_pInt,nt constitutive_dislotwin_homogenizedC = & - constitutive_dislotwin_homogenizedC + state(g,ip,el)%p(2*ns+i)*constitutive_dislotwin_Ctwin_66(:,:,i,myInstance) + constitutive_dislotwin_homogenizedC + state(g,ip,el)%p(2_pInt*ns+i)*constitutive_dislotwin_Ctwin_66(:,:,i,myInstance) enddo return @@ -808,7 +809,6 @@ use prec, only: pReal,pInt,p_vec use math, only: pi use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance -use lattice, only: lattice_interactionSlipTwin,lattice_interactionTwinTwin !use debug, only: debugger implicit none @@ -846,69 +846,69 @@ sfe = constitutive_dislotwin_SFE_0K(myInstance) + & constitutive_dislotwin_dSFE_dT(myInstance) * Temperature !* rescaled twin volume fraction for topology -forall (t = 1:nt) & +forall (t = 1_pInt:nt) & fOverStacksize(t) = & - state(g,ip,el)%p(2*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance) + state(g,ip,el)%p(2_pInt*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance) !* 1/mean free distance between 2 forest dislocations seen by a moving dislocation -forall (s = 1:ns) & - state(g,ip,el)%p(2*ns+nt+s) = & - sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),& +forall (s = 1_pInt:ns) & + state(g,ip,el)%p(2_pInt*ns+nt+s) = & + sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1_pInt:2_pInt*ns)),& constitutive_dislotwin_forestProjectionEdge(1:ns,s,myInstance)))/ & constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,myInstance) !* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation !$OMP CRITICAL (evilmatmul) -state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = 0.0_pReal +state(g,ip,el)%p((3_pInt*ns+nt+1_pInt):(4_pInt*ns+nt)) = 0.0_pReal if (nt > 0_pInt) & - state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = & + state(g,ip,el)%p((3_pInt*ns+nt+1):(4_pInt*ns+nt)) = & matmul(constitutive_dislotwin_interactionMatrixSlipTwin(1:ns,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf) !$OMP END CRITICAL (evilmatmul) !* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin !$OMP CRITICAL (evilmatmul) if (nt > 0_pInt) & - state(g,ip,el)%p((4*ns+nt+1):(4*ns+2*nt)) = & + state(g,ip,el)%p((4_pInt*ns+nt+1_pInt):(4_pInt*ns+2_pInt*nt)) = & matmul(constitutive_dislotwin_interactionMatrixTwinTwin(1:nt,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf) !$OMP END CRITICAL (evilmatmul) !* mean free path between 2 obstacles seen by a moving dislocation -do s = 1,ns +do s = 1_pInt,ns if (nt > 0_pInt) then - state(g,ip,el)%p(4*ns+2*nt+s) = & + state(g,ip,el)%p(4_pInt*ns+2_pInt*nt+s) = & constitutive_dislotwin_GrainSize(myInstance)/(1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*& - (state(g,ip,el)%p(2*ns+nt+s)+state(g,ip,el)%p(3*ns+nt+s))) + (state(g,ip,el)%p(2_pInt*ns+nt+s)+state(g,ip,el)%p(3_pInt*ns+nt+s))) else - state(g,ip,el)%p(4*ns+s) = & + state(g,ip,el)%p(4_pInt*ns+s) = & constitutive_dislotwin_GrainSize(myInstance)/& - (1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*(state(g,ip,el)%p(2*ns+s))) + (1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*(state(g,ip,el)%p(2_pInt*ns+s))) endif enddo !* mean free path between 2 obstacles seen by a growing twin -forall (t = 1:nt) & - state(g,ip,el)%p(5*ns+2*nt+t) = & +forall (t = 1_pInt:nt) & + state(g,ip,el)%p(5_pInt*ns+2_pInt*nt+t) = & (constitutive_dislotwin_Cmfptwin(myInstance)*constitutive_dislotwin_GrainSize(myInstance))/& - (1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*state(g,ip,el)%p(4*ns+nt+t)) + (1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*state(g,ip,el)%p(4_pInt*ns+nt+t)) !* threshold stress for dislocation motion -forall (s = 1:ns) & - state(g,ip,el)%p(5*ns+3*nt+s) = constitutive_dislotwin_SolidSolutionStrength(myInstance)+ & +forall (s = 1_pInt:ns) & + state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+s) = constitutive_dislotwin_SolidSolutionStrength(myInstance)+ & constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)*& - sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),& + sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1_pInt:2_pInt*ns)),& constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance))) !* threshold stress for growing twin -forall (t = 1:nt) & - state(g,ip,el)%p(6*ns+3*nt+t) = & +forall (t = 1_pInt:nt) & + state(g,ip,el)%p(6_pInt*ns+3_pInt*nt+t) = & constitutive_dislotwin_Cthresholdtwin(myInstance)*& (sfe/(3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,myInstance))+& 3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,myInstance)*constitutive_dislotwin_Gmod(myInstance)/& (constitutive_dislotwin_L0(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(t,myInstance))) !* final twin volume after growth -forall (t = 1:nt) & - state(g,ip,el)%p(6*ns+4*nt+t) = & +forall (t = 1_pInt:nt) & + state(g,ip,el)%p(6_pInt*ns+4_pInt*nt+t) = & (pi/6.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)*state(g,ip,el)%p(5*ns+2*nt+t)**(2.0_pReal) !if ((ip==1).and.(el==1)) then @@ -993,7 +993,7 @@ ns = constitutive_dislotwin_totalNslip(myInstance) nt = constitutive_dislotwin_totalNtwin(myInstance) !* Total twin volume fraction -sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 +sumf = sum(state(g,ip,el)%p((2_pInt*ns+1_pInt):(2_pInt*ns+nt))) ! safe for nt == 0 Lp = 0.0_pReal dLp_dTstar3333 = 0.0_pReal @@ -1003,9 +1003,9 @@ dLp_dTstar = 0.0_pReal gdot_slip = 0.0_pReal dgdot_dtauslip = 0.0_pReal j = 0_pInt -do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family +do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j+1_pInt !* Calculation of Lp @@ -1036,7 +1036,7 @@ do f = 1,lattice_maxNslipFamily ! loop over all Lp = Lp + (1.0_pReal - sumf)*gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,myStructure) !* Calculation of the tangent of Lp - forall (k=1:3,l=1:3,m=1:3,n=1:3) & + 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(j)*& lattice_Sslip(k,l,index_myFamily+i,myStructure)*& @@ -1049,7 +1049,7 @@ if(constitutive_dislotwin_sbVelocity(myInstance) /= 0.0_pReal) then gdot_sb = 0.0_pReal dgdot_dtausb = 0.0_pReal call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error) - do j = 1,6 + do j = 1_pInt,6_pInt sb_s = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_sComposition(1:3,j)) sb_m = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_mComposition(1:3,j)) sb_Smatrix = math_tensorproduct(sb_s,sb_m) @@ -1073,20 +1073,20 @@ if(constitutive_dislotwin_sbVelocity(myInstance) /= 0.0_pReal) then DotGamma0 = constitutive_dislotwin_sbVelocity(myInstance) !* Shear rates due to shearband - gdot_sb(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*& + gdot_sb(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_q(myInstance))*& sign(1.0_pReal,tau_sb(j)) !* Derivatives of shear rates dgdot_dtausb(j) = & ((abs(gdot_sb(j))*BoltzmannRatio*& constitutive_dislotwin_p(myInstance)*constitutive_dislotwin_q(myInstance))/constitutive_dislotwin_sbResistance(myInstance))*& - StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal) + StressRatio_pminus1*(1_pInt-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal) !* Plastic velocity gradient for shear banding Lp = Lp + gdot_sb(j)*sb_Smatrix !* Calculation of the tangent of Lp - forall (k=1:3,l=1:3,m=1:3,n=1:3) & + 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_dtausb(j)*& sb_Smatrix(k,l)*& @@ -1098,9 +1098,9 @@ end if gdot_twin = 0.0_pReal dgdot_dtautwin = 0.0_pReal j = 0_pInt -do f = 1,lattice_maxNtwinFamily ! loop over all slip families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family +do f = 1_pInt,lattice_maxNtwinFamily ! loop over all slip families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family j = j+1_pInt !* Calculation of Lp @@ -1122,7 +1122,7 @@ do f = 1,lattice_maxNtwinFamily ! loop over all Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,myStructure) !* Calculation of the tangent of Lp - forall (k=1:3,l=1:3,m=1:3,n=1:3) & + 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(j)*& lattice_Stwin(k,l,index_myFamily+i,myStructure)*& @@ -1164,11 +1164,9 @@ use prec, only: pReal,pInt,p_vec use math, only: pi use mesh, only: mesh_NcpElems, mesh_maxNips use material, only: homogenization_maxNgrains, material_phase, phase_constitutionInstance -use lattice, only: lattice_Sslip,lattice_Sslip_v, & - lattice_Stwin,lattice_Stwin_v, & +use lattice, only: lattice_Sslip_v, lattice_Stwin_v, & lattice_maxNslipFamily,lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem, & - lattice_shearTwin + lattice_NslipSystem,lattice_NtwinSystem implicit none !* Input-Output variables @@ -1196,24 +1194,26 @@ ns = constitutive_dislotwin_totalNslip(myInstance) nt = constitutive_dislotwin_totalNtwin(myInstance) !* Total twin volume fraction -sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 +sumf = sum(state(g,ip,el)%p((2_pInt*ns+1_pInt):(2_pInt*ns+nt))) ! safe for nt == 0 constitutive_dislotwin_dotState = 0.0_pReal !* Dislocation density evolution gdot_slip = 0.0_pReal j = 0_pInt -do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,MyStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family +do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,MyStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j+1_pInt !* Resolved shear stress on slip system tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) !* Stress ratios - StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance) - StressRatio_pminus1 = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) + StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& + constitutive_dislotwin_p(myInstance) + StressRatio_pminus1 = (abs(tau_slip(j))/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& + (constitutive_dislotwin_p(myInstance)-1.0_pReal) !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) !* Initial shear rates @@ -1222,7 +1222,7 @@ do f = 1,lattice_maxNslipFamily ! loop over all constitutive_dislotwin_v0PerSlipSystem(f,myInstance) !* Shear rates due to slip - gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*& + gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_q(myInstance))*& sign(1.0_pReal,tau_slip(j)) !* Multiplication @@ -1283,9 +1283,9 @@ enddo !* Twin volume fraction evolution j = 0_pInt -do f = 1,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,MyStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) twin system in family +do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,MyStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) twin system in family j = j+1_pInt !* Resolved shear stress on twin system @@ -1295,9 +1295,9 @@ do f = 1,lattice_maxNtwinFamily ! loop over all !* Shear rates and their derivatives due to twin if ( tau_twin(j) > 0.0_pReal ) then - constitutive_dislotwin_dotState(2*ns+j) = & + constitutive_dislotwin_dotState(2_pInt*ns+j) = & (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*& - state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) + state(g,ip,el)%p(6_pInt*ns+4_pInt*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) endif enddo @@ -1365,7 +1365,7 @@ use math, only: pi,math_Mandel6to33, math_spectralDecompositionSym33 use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin + lattice_NslipSystem,lattice_NtwinSystem implicit none !* Definition of variables @@ -1397,27 +1397,29 @@ constitutive_dislotwin_postResults = 0.0_pReal !* Spectral decomposition of stress call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error) -do o = 1,phase_Noutput(material_phase(g,ip,el)) +do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) select case(constitutive_dislotwin_output(o,myInstance)) case ('edge_density') - constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) + constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p(1:ns) c = c + ns case ('dipole_density') - constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p(ns+1:2*ns) + constitutive_dislotwin_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p(ns+1:2_pInt*ns) c = c + ns case ('shear_rate_slip') j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j + 1_pInt !* Resolved shear stress on slip system tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) !* Stress ratios - StressRatio_p = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance) - StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) + StressRatio_p = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& + constitutive_dislotwin_p(myInstance) + StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& + (constitutive_dislotwin_p(myInstance)-1.0_pReal) !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) !* Initial shear rates @@ -1427,29 +1429,33 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) !* Shear rates due to slip constitutive_dislotwin_postResults(c+j) = & - DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau) + DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& + constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau) enddo ; enddo c = c + ns case ('mfp_slip') - constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((4*ns+2*nt+1):(5*ns+2*nt)) + constitutive_dislotwin_postResults(c+1_pInt:c+ns) =& + state(g,ip,el)%p((4_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+2_pInt*nt)) c = c + ns case ('resolved_stress_slip') j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j + 1_pInt - constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) + constitutive_dislotwin_postResults(c+j) =& + dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) enddo; enddo c = c + ns case ('threshold_stress_slip') - constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((5*ns+3*nt+1):(6*ns+3*nt)) + constitutive_dislotwin_postResults(c+1_pInt:c+ns) = & + state(g,ip,el)%p((5_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+3_pInt*nt)) c = c + ns case ('edge_dipole_distance') j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j + 1_pInt constitutive_dislotwin_postResults(c+j) = & (3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/& @@ -1459,15 +1465,15 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) enddo; enddo c = c + ns case ('resolved_stress_shearband') - do j = 1,6 ! loop over all shearband families + do j = 1_pInt,6_pInt ! loop over all shearband families constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v, constitutive_dislotwin_sbSv(1:6,j,g,ip,el)) enddo - c = c + 6 + c = c + 6_pInt case ('schmid_factor_shearband') - constitutive_dislotwin_postResults(c+1:c+6) = constitutive_dislotwin_sbSv(1:6,j,g,ip,el) - c = c + 6 + constitutive_dislotwin_postResults(c+1_pInt:c+6_pInt) = constitutive_dislotwin_sbSv(1:6,j,g,ip,el) + c = c + 6_pInt case ('shear_rate_shearband') - do j = 1,6 ! loop over all shearband families + do j = 1_pInt,6_pInt ! loop over all shearband families !* Resolved shear stress on shearband system tau = dot_product(Tstar_v,constitutive_dislotwin_sbSv(1:6,j,g,ip,el)) !* Stress ratios @@ -1481,64 +1487,66 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) !* Shear rates due to slip constitutive_dislotwin_postResults(c+j) = & - DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau) + DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau) enddo - c = c + 6 + c = c + 6_pInt case ('twin_fraction') - constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((2*ns+1):(2*ns+nt)) + constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(g,ip,el)%p((2_pInt*ns+1_pInt):(2_pInt*ns+nt)) c = c + nt case ('shear_rate_twin') if (nt > 0_pInt) then j = 0_pInt - do f = 1,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family + do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) ! at which index starts my family do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) twin system in family j = j + 1_pInt !* Resolved shear stress on twin system tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) !* Stress ratios - StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau)**constitutive_dislotwin_r(myInstance) + StressRatio_r = (state(g,ip,el)%p(6_pInt*ns+3_pInt*nt+j)/tau)**constitutive_dislotwin_r(myInstance) !* Shear rates and their derivatives due to twin if ( tau > 0.0_pReal ) then constitutive_dislotwin_postResults(c+j) = & (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*& - state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) + state(g,ip,el)%p(6_pInt*ns+4_pInt*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) endif enddo ; enddo endif c = c + nt case ('mfp_twin') - constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((5*ns+2*nt+1):(5*ns+3*nt)) + constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(g,ip,el)%p((5_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+3_pInt*nt)) c = c + nt case ('resolved_stress_twin') if (nt > 0_pInt) then j = 0_pInt - do f = 1,lattice_maxNtwinFamily ! loop over all slip families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNtwinFamily ! loop over all slip families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family j = j + 1_pInt constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) enddo; enddo endif c = c + nt case ('threshold_stress_twin') - constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((6*ns+3*nt+1):(6*ns+4*nt)) + constitutive_dislotwin_postResults(c+1_pInt:c+nt) = state(g,ip,el)%p((6_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+4_pInt*nt)) c = c + nt case ('stress_exponent') j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j + 1_pInt !* Resolved shear stress on slip system tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) !* Stress ratios - StressRatio_p = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance) - StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) + StressRatio_p = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& + constitutive_dislotwin_p(myInstance) + StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**& + (constitutive_dislotwin_p(myInstance)-1.0_pReal) !* Boltzmann ratio BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) !* Initial shear rates @@ -1547,14 +1555,14 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) constitutive_dislotwin_v0PerSlipSystem(f,myInstance) !* Shear rates due to slip - gdot_slip = & - DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau) + gdot_slip = DotGamma0*exp(-BoltzmannRatio*(1_pInt-StressRatio_p)**& + constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau) !* Derivatives of shear rates dgdot_dtauslip = & ((abs(gdot_slip)*BoltzmannRatio*& constitutive_dislotwin_p(myInstance)*constitutive_dislotwin_q(myInstance))/state(g,ip,el)%p(5*ns+3*nt+j))*& - StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal) + StressRatio_pminus1*(1_pInt-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal) !* Stress exponent if (gdot_slip==0.0_pReal) then @@ -1565,12 +1573,12 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) enddo ; enddo c = c + ns case ('sb_eigenvalues') - forall (j = 1:3) & + forall (j = 1_pInt:3_pInt) & constitutive_dislotwin_postResults(c+j) = eigValues(j) - c = c + 3 + c = c + 3_pInt case ('sb_eigenvectors') - constitutive_dislotwin_postResults(c+1:c+9) = reshape(eigVectors,(/9/)) - c = c + 9 + constitutive_dislotwin_postResults(c+1_pInt:c+9_pInt) = reshape(eigVectors,(/9/)) + c = c + 9_pInt end select enddo diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index a6dd53d83..33dca8f5b 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -84,14 +84,15 @@ subroutine constitutive_j2_init(file) !************************************** !* Module initialization * !************************************** + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pInt, pReal use math, only: math_Mandel3333to66, math_Voigt66to3333 use IO use material use debug, only: debug_verbosity integer(pInt), intent(in) :: file - integer(pInt), parameter :: maxNchunks = 7 - integer(pInt), dimension(1+2*maxNchunks) :: positions + integer(pInt), parameter :: maxNchunks = 7_pInt + integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions integer(pInt) section, maxNinstance, i,j,k, mySize character(len=64) tag character(len=1024) line @@ -103,10 +104,10 @@ subroutine constitutive_j2_init(file) #include "compilation_info.f90" !$OMP END CRITICAL (write2out) - maxNinstance = count(phase_constitution == constitutive_j2_label) - if (maxNinstance == 0) return + maxNinstance = int(count(phase_constitution == constitutive_j2_label),pInt) + if (maxNinstance == 0_pInt) return - if (debug_verbosity > 0) then + if (debug_verbosity > 0_pInt) then !$OMP CRITICAL (write2out) write(6,'(a16,1x,i5)') '# instances:',maxNinstance write(6,*) @@ -133,7 +134,7 @@ subroutine constitutive_j2_init(file) rewind(file) line = '' - section = 0 + section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to read(file,'(a1024)',END=100) line @@ -156,34 +157,34 @@ subroutine constitutive_j2_init(file) cycle case ('(output)') constitutive_j2_Noutput(i) = constitutive_j2_Noutput(i) + 1_pInt - constitutive_j2_output(constitutive_j2_Noutput(i),i) = IO_lc(IO_stringValue(line,positions,2)) + constitutive_j2_output(constitutive_j2_Noutput(i),i) = IO_lc(IO_stringValue(line,positions,2_pInt)) case ('c11') - constitutive_j2_C11(i) = IO_floatValue(line,positions,2) + constitutive_j2_C11(i) = IO_floatValue(line,positions,2_pInt) case ('c12') - constitutive_j2_C12(i) = IO_floatValue(line,positions,2) + constitutive_j2_C12(i) = IO_floatValue(line,positions,2_pInt) case ('tau0') - constitutive_j2_tau0(i) = IO_floatValue(line,positions,2) + constitutive_j2_tau0(i) = IO_floatValue(line,positions,2_pInt) case ('gdot0') - constitutive_j2_gdot0(i) = IO_floatValue(line,positions,2) + constitutive_j2_gdot0(i) = IO_floatValue(line,positions,2_pInt) case ('n') - constitutive_j2_n(i) = IO_floatValue(line,positions,2) + constitutive_j2_n(i) = IO_floatValue(line,positions,2_pInt) case ('h0') - constitutive_j2_h0(i) = IO_floatValue(line,positions,2) + constitutive_j2_h0(i) = IO_floatValue(line,positions,2_pInt) case ('tausat') - constitutive_j2_tausat(i) = IO_floatValue(line,positions,2) + constitutive_j2_tausat(i) = IO_floatValue(line,positions,2_pInt) case ('a', 'w0') - constitutive_j2_a(i) = IO_floatValue(line,positions,2) + constitutive_j2_a(i) = IO_floatValue(line,positions,2_pInt) case ('taylorfactor') - constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2) + constitutive_j2_fTaylor(i) = IO_floatValue(line,positions,2_pInt) case ('atol_resistance') - constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2) + constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2_pInt) case default call IO_error(210_pInt,ext_msg=tag) end select endif enddo -100 do i = 1,maxNinstance ! sanity checks +100 do i = 1_pInt,maxNinstance ! sanity checks if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(211_pInt,ext_msg='tau0') if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='gdot0') if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='n') @@ -193,8 +194,8 @@ subroutine constitutive_j2_init(file) if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='aTol_resistance') enddo - do i = 1,maxNinstance - do j = 1,constitutive_j2_Noutput(i) + do i = 1_pInt,maxNinstance + do j = 1_pInt,constitutive_j2_Noutput(i) select case(constitutive_j2_output(j,i)) case('flowstress') mySize = 1_pInt @@ -211,13 +212,14 @@ subroutine constitutive_j2_init(file) endif enddo - constitutive_j2_sizeDotState(i) = 1 - constitutive_j2_sizeState(i) = 1 + constitutive_j2_sizeDotState(i) = 1_pInt + constitutive_j2_sizeState(i) = 1_pInt - forall(k=1:3) - forall(j=1:3) & - constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i) - constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i) + forall(k=1_pInt:3_pInt) + forall(j=1_pInt:3_pInt) + constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i) + end forall + constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i) constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i)) end forall constitutive_j2_Cslip_66(1:6,1:6,i) = & @@ -333,8 +335,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v, use math, only: math_mul6x6, & math_Mandel6to33, & math_Plain3333to99 - use lattice, only: lattice_Sslip, & - lattice_Sslip_v use mesh, only: mesh_NcpElems, & mesh_maxNips use material, only: homogenization_maxNgrains, & @@ -379,7 +379,7 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v, dLp_dTstar_99 = 0.0_pReal ! for Tstar==0 both Lp and dLp_dTstar are zero (if not n==1) - if (norm_Tstar_dev > 0) then + if (norm_Tstar_dev > 0_pInt) then ! Calculation of gamma_dot gamma_dot = constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & @@ -390,9 +390,9 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v, Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/constitutive_j2_fTaylor(matID) !* Calculation of the tangent of Lp - forall (k=1:3,l=1:3,m=1:3,n=1:3) & + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar_3333(k,l,m,n) = (constitutive_j2_n(matID)-1.0_pReal) * Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev - forall (k=1:3,l=1:3) & + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal dLp_dTstar_99 = math_Plain3333to99(gamma_dot / constitutive_j2_fTaylor(matID) * dLp_dTstar_3333 / norm_Tstar_dev) end if @@ -412,7 +412,6 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el) pInt, & p_vec use math, only: math_mul6x6 - use lattice, only: lattice_Sslip_v use mesh, only: mesh_NcpElems, & mesh_maxNips use material, only: homogenization_maxNgrains, & @@ -471,9 +470,8 @@ pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip, !*** variables and functions from other modules ***! use prec, only: pReal,pInt,p_vec - use lattice, only: lattice_Sslip_v use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance + use material, only: homogenization_maxNgrains implicit none !*** input variables ***! @@ -504,7 +502,6 @@ pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip pInt, & p_vec use math, only: math_mul6x6 - use lattice, only: lattice_Sslip_v use mesh, only: mesh_NcpElems, & mesh_maxNips use material, only: homogenization_maxNgrains, & @@ -550,17 +547,17 @@ pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip c = 0_pInt constitutive_j2_postResults = 0.0_pReal - do o = 1,phase_Noutput(material_phase(g,ip,el)) + do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) select case(constitutive_j2_output(o,matID)) case ('flowstress') - constitutive_j2_postResults(c+1) = state(g,ip,el)%p(1) - c = c + 1 + constitutive_j2_postResults(c+1_pInt) = state(g,ip,el)%p(1) + c = c + 1_pInt case ('strainrate') - constitutive_j2_postResults(c+1) = & + constitutive_j2_postResults(c+1_pInt) = & constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & / &!--------------------------------------------------- (constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) ** constitutive_j2_n(matID) - c = c + 1 + c = c + 1_pInt end select enddo diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index 8b8c7af86..697df4dc0 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -145,10 +145,11 @@ CONTAINS !**************************************** -subroutine constitutive_phenopowerlaw_init(file) +subroutine constitutive_phenopowerlaw_init(myFile) !************************************** !* Module initialization * !************************************** + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pInt, pReal use math, only: math_Mandel3333to66, math_Voigt66to3333 use IO @@ -163,7 +164,7 @@ subroutine constitutive_phenopowerlaw_init(file) lattice_interactionTwinSlip, & lattice_interactionTwinTwin - integer(pInt), intent(in) :: file + integer(pInt), intent(in) :: myFile integer(pInt), parameter :: maxNchunks = lattice_maxNinteraction + 1_pInt integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt) section, maxNinstance, i,j,k, f,o, & @@ -178,7 +179,7 @@ subroutine constitutive_phenopowerlaw_init(file) #include "compilation_info.f90" !$OMP END CRITICAL (write2out) - maxNinstance = count(phase_constitution == constitutive_phenopowerlaw_label) + maxNinstance = int(count(phase_constitution == constitutive_phenopowerlaw_label),pInt) if (maxNinstance == 0) return if (debug_verbosity > 0) then @@ -252,16 +253,16 @@ subroutine constitutive_phenopowerlaw_init(file) allocate(constitutive_phenopowerlaw_aTolResistance(maxNinstance)) constitutive_phenopowerlaw_aTolResistance = 0.0_pReal - rewind(file) + rewind(myFile) line = '' - section = 0 + section = 0_pInt do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line enddo do ! read thru sections of phase part - read(file,'(a1024)',END=100) line + read(myFile,'(a1024)',END=100) line if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'[',']') /= '') then ! next section @@ -542,14 +543,14 @@ function constitutive_phenopowerlaw_stateInit(myInstance) constitutive_phenopowerlaw_stateInit = 0.0_pReal - do i = 1,lattice_maxNslipFamily + do i = 1_pInt,lattice_maxNslipFamily constitutive_phenopowerlaw_stateInit(1+& sum(constitutive_phenopowerlaw_Nslip(1:i-1,myInstance)) : & sum(constitutive_phenopowerlaw_Nslip(1:i ,myInstance))) = & constitutive_phenopowerlaw_tau0_slip(i,myInstance) enddo - do i = 1,lattice_maxNtwinFamily + do i = 1_pInt,lattice_maxNtwinFamily constitutive_phenopowerlaw_stateInit(1+sum(constitutive_phenopowerlaw_Nslip(:,myInstance))+& sum(constitutive_phenopowerlaw_Ntwin(1:i-1,myInstance)) : & sum(constitutive_phenopowerlaw_Nslip(:,myInstance))+& @@ -677,17 +678,17 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp nSlip = constitutive_phenopowerlaw_totalNslip(matID) nTwin = constitutive_phenopowerlaw_totalNtwin(matID) - index_Gamma = nSlip + nTwin + 1 - index_F = nSlip + nTwin + 2 + index_Gamma = nSlip + nTwin + 1_pInt + index_F = nSlip + nTwin + 2_pInt Lp = 0.0_pReal dLp_dTstar3333 = 0.0_pReal dLp_dTstar = 0.0_pReal j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family - do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family j = j+1_pInt !* Calculation of Lp @@ -702,7 +703,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp if (gdot_slip(j) /= 0.0_pReal) then dgdot_dtauslip(j) = gdot_slip(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip(j) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & + 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(j)*lattice_Sslip(k,l,index_myFamily+i,structID)* & lattice_Sslip(m,n,index_myFamily+i,structID) @@ -711,9 +712,9 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp enddo j = 0_pInt - do f = 1,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family - do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family j = j+1_pInt !* Calculation of Lp @@ -729,7 +730,7 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp if (gdot_twin(j) /= 0.0_pReal) then dgdot_dtautwin(j) = gdot_twin(j)*constitutive_phenopowerlaw_n_twin(matID)/tau_twin(j) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & + 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(j)*lattice_Stwin(k,l,index_myFamily+i,structID)* & lattice_Stwin(m,n,index_myFamily+i,structID) @@ -755,7 +756,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el !* - constitutive_dotState : evolution of state variable * !********************************************************************* use prec, only: pReal,pInt,p_vec - use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & + use lattice, only: lattice_Sslip_v, lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance @@ -780,8 +781,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el nSlip = constitutive_phenopowerlaw_totalNslip(matID) nTwin = constitutive_phenopowerlaw_totalNtwin(matID) - index_Gamma = nSlip + nTwin + 1 - index_F = nSlip + nTwin + 2 + index_Gamma = nSlip + nTwin + 1_pInt + index_F = nSlip + nTwin + 2_pInt constitutive_phenopowerlaw_dotState = 0.0_pReal @@ -800,9 +801,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el ssat_offset = constitutive_phenopowerlaw_spr(matID)*sqrt(state(ipc,ip,el)%p(index_F)) j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family - do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family j = j+1_pInt h_slipslip(j) = c_slipslip*(1.0_pReal-state(ipc,ip,el)%p(j) / & ! system-dependent prefactor for slip--slip interaction (constitutive_phenopowerlaw_tausat_slip(f,matID)+ssat_offset))** & @@ -819,9 +820,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el enddo j = 0_pInt - do f = 1,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family - do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family j = j+1_pInt h_twinslip(j) = c_twinslip ! no system-dependent parts h_twintwin(j) = c_twintwin @@ -839,8 +840,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el !-- calculate the overall hardening based on above j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family j = j+1_pInt constitutive_phenopowerlaw_dotState(j) = & ! evolution of slip resistance j h_slipslip(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_slipslip(1:nSlip,j,matID),abs(gdot_slip)) + & ! dot gamma_slip @@ -851,9 +852,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el enddo j = 0_pInt - do f = 1,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family - do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family j = j+1_pInt constitutive_phenopowerlaw_dotState(j+nSlip) = & ! evolution of twin resistance j h_twinslip(j) * dot_product(constitutive_phenopowerlaw_hardeningMatrix_twinslip(1:nSlip,j,matID),abs(gdot_slip)) + & ! dot gamma_slip @@ -875,9 +876,8 @@ pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,stat !*** variables and functions from other modules ***! use prec, only: pReal,pInt,p_vec - use lattice, only: lattice_Sslip_v - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance + use mesh, only: mesh_NcpElems, mesh_maxNips + use material, only: homogenization_maxNgrains implicit none !*** input variables ***! @@ -932,23 +932,23 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat nSlip = constitutive_phenopowerlaw_totalNslip(matID) nTwin = constitutive_phenopowerlaw_totalNtwin(matID) - index_Gamma = nSlip + nTwin + 1 - index_F = nSlip + nTwin + 2 + index_Gamma = nSlip + nTwin + 1_pInt + index_F = nSlip + nTwin + 2_pInt constitutive_phenopowerlaw_postResults = 0.0_pReal c = 0_pInt - do o = 1,phase_Noutput(material_phase(ipc,ip,el)) + do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) select case(constitutive_phenopowerlaw_output(o,matID)) case ('resistance_slip') - constitutive_phenopowerlaw_postResults(c+1:c+nSlip) = state(ipc,ip,el)%p(1:nSlip) + constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(ipc,ip,el)%p(1:nSlip) c = c + nSlip case ('shearrate_slip') j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family - do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family j = j + 1_pInt tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID)) constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(matID)*& @@ -959,27 +959,27 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat case ('resolvedstress_slip') j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family - do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family j = j + 1_pInt constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID)) enddo; enddo c = c + nSlip case ('totalshear') - constitutive_phenopowerlaw_postResults(c+1) = state(ipc,ip,el)%p(index_Gamma) - c = c + 1 + constitutive_phenopowerlaw_postResults(c+1_pInt) = state(ipc,ip,el)%p(index_Gamma) + c = c + 1_pInt case ('resistance_twin') - constitutive_phenopowerlaw_postResults(c+1:c+nTwin) = state(ipc,ip,el)%p(1+nSlip:nTwin+nSlip) + constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = state(ipc,ip,el)%p(1_pInt+nSlip:nTwin+nSlip) c = c + nTwin case ('shearrate_twin') j = 0_pInt - do f = 1,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family - do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family j = j + 1_pInt tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID)) constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F @@ -991,17 +991,17 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat case ('resolvedstress_twin') j = 0_pInt - do f = 1,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family - do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family + do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family j = j + 1_pInt constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID)) enddo; enddo c = c + nTwin case ('totalvolfrac') - constitutive_phenopowerlaw_postResults(c+1) = state(ipc,ip,el)%p(index_F) - c = c + 1 + constitutive_phenopowerlaw_postResults(c+1_pInt) = state(ipc,ip,el)%p(index_F) + c = c + 1_pInt end select enddo diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index 2e9e49c46..a1990fff6 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -203,7 +203,7 @@ subroutine constitutive_titanmod_init(file) !************************************** !* Module initialization * !************************************** -use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) +use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use prec, only: pInt,pReal use math, only: math_Mandel3333to66,math_Voigt66to3333,math_mul3x3 use IO @@ -213,10 +213,11 @@ use lattice !* Input variables integer(pInt), intent(in) :: file !* Local variables -integer(pInt), parameter :: maxNchunks = 21 -integer(pInt), dimension(1+2*maxNchunks) :: positions -integer(pInt) section,maxNinstance,f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t,t1,t2,ns,nt,mySize,myStructure,maxTotalNslip, & +integer(pInt), parameter :: maxNchunks = 21_pInt +integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions +integer(pInt) section,f,i,j,k,l,m,n,o,p,q,r,s,s1,s2,t,t1,t2,ns,nt,mySize,myStructure,maxTotalNslip, & maxTotalNtwin +integer :: maxNinstance !no pInt character(len=64) tag character(len=1024) line @@ -229,158 +230,156 @@ maxNinstance = count(phase_constitution == constitutive_titanmod_label) if (maxNinstance == 0) return !* Space allocation for global variables -allocate(constitutive_titanmod_sizeDotState(maxNinstance)) +allocate(constitutive_titanmod_sizeDotState(maxNinstance)) + constitutive_titanmod_sizeDotState = 0_pInt allocate(constitutive_titanmod_sizeState(maxNinstance)) + constitutive_titanmod_sizeState = 0_pInt allocate(constitutive_titanmod_sizePostResults(maxNinstance)) + constitutive_titanmod_sizePostResults = 0_pInt allocate(constitutive_titanmod_sizePostResult(maxval(phase_Noutput),maxNinstance)) + constitutive_titanmod_sizePostResult = 0_pInt allocate(constitutive_titanmod_output(maxval(phase_Noutput),maxNinstance)) + constitutive_titanmod_output = '' allocate(constitutive_titanmod_Noutput(maxNinstance)) -constitutive_titanmod_sizeDotState = 0_pInt -constitutive_titanmod_sizeState = 0_pInt -constitutive_titanmod_sizePostResults = 0_pInt -constitutive_titanmod_sizePostResult = 0_pInt -constitutive_titanmod_output = '' -constitutive_titanmod_Noutput = 0_pInt + constitutive_titanmod_Noutput = 0_pInt allocate(constitutive_titanmod_structureName(maxNinstance)) + constitutive_titanmod_structureName = '' allocate(constitutive_titanmod_structure(maxNinstance)) + constitutive_titanmod_structure = 0_pInt allocate(constitutive_titanmod_Nslip(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_Nslip = 0_pInt allocate(constitutive_titanmod_Ntwin(lattice_maxNtwinFamily,maxNinstance)) + constitutive_titanmod_Ntwin = 0_pInt allocate(constitutive_titanmod_slipFamily(lattice_maxNslip,maxNinstance)) + constitutive_titanmod_slipFamily = 0_pInt allocate(constitutive_titanmod_twinFamily(lattice_maxNtwin,maxNinstance)) + constitutive_titanmod_twinFamily = 0_pInt allocate(constitutive_titanmod_slipSystemLattice(lattice_maxNslip,maxNinstance)) + constitutive_titanmod_slipSystemLattice = 0_pInt allocate(constitutive_titanmod_twinSystemLattice(lattice_maxNtwin,maxNinstance)) + constitutive_titanmod_twinSystemLattice = 0_pInt allocate(constitutive_titanmod_totalNslip(maxNinstance)) + constitutive_titanmod_totalNslip = 0_pInt allocate(constitutive_titanmod_totalNtwin(maxNinstance)) -constitutive_titanmod_structureName = '' -constitutive_titanmod_structure = 0_pInt -constitutive_titanmod_Nslip = 0_pInt -constitutive_titanmod_Ntwin = 0_pInt -constitutive_titanmod_slipFamily = 0_pInt -constitutive_titanmod_twinFamily = 0_pInt -constitutive_titanmod_slipSystemLattice = 0_pInt -constitutive_titanmod_twinSystemLattice = 0_pInt -constitutive_titanmod_totalNslip = 0_pInt -constitutive_titanmod_totalNtwin = 0_pInt + constitutive_titanmod_totalNtwin = 0_pInt allocate(constitutive_titanmod_CoverA(maxNinstance)) + constitutive_titanmod_CoverA = 0.0_pReal allocate(constitutive_titanmod_C11(maxNinstance)) + constitutive_titanmod_C11 = 0.0_pReal allocate(constitutive_titanmod_C12(maxNinstance)) + constitutive_titanmod_C12 = 0.0_pReal allocate(constitutive_titanmod_C13(maxNinstance)) + constitutive_titanmod_C13 = 0.0_pReal allocate(constitutive_titanmod_C33(maxNinstance)) + constitutive_titanmod_C33 = 0.0_pReal allocate(constitutive_titanmod_C44(maxNinstance)) + constitutive_titanmod_C44 = 0.0_pReal allocate(constitutive_titanmod_debyefrequency(maxNinstance)) + constitutive_titanmod_debyefrequency = 0.0_pReal allocate(constitutive_titanmod_kinkf0(maxNinstance)) + constitutive_titanmod_kinkf0 = 0.0_pReal allocate(constitutive_titanmod_Gmod(maxNinstance)) + constitutive_titanmod_Gmod = 0.0_pReal allocate(constitutive_titanmod_CAtomicVolume(maxNinstance)) + constitutive_titanmod_CAtomicVolume = 0.0_pReal allocate(constitutive_titanmod_dc(maxNinstance)) + constitutive_titanmod_dc = 0.0_pReal allocate(constitutive_titanmod_twinhpconstant(maxNinstance)) + constitutive_titanmod_twinhpconstant = 0.0_pReal allocate(constitutive_titanmod_GrainSize(maxNinstance)) + constitutive_titanmod_GrainSize = 0.0_pReal allocate(constitutive_titanmod_MaxTwinFraction(maxNinstance)) + constitutive_titanmod_MaxTwinFraction = 0.0_pReal allocate(constitutive_titanmod_r(maxNinstance)) + constitutive_titanmod_r = 0.0_pReal allocate(constitutive_titanmod_CEdgeDipMinDistance(maxNinstance)) + constitutive_titanmod_CEdgeDipMinDistance = 0.0_pReal allocate(constitutive_titanmod_Cmfptwin(maxNinstance)) + constitutive_titanmod_Cmfptwin = 0.0_pReal allocate(constitutive_titanmod_Cthresholdtwin(maxNinstance)) + constitutive_titanmod_Cthresholdtwin = 0.0_pReal allocate(constitutive_titanmod_aTolRho(maxNinstance)) + constitutive_titanmod_aTolRho = 0.0_pReal allocate(constitutive_titanmod_Cslip_66(6,6,maxNinstance)) + constitutive_titanmod_Cslip_66 = 0.0_pReal allocate(constitutive_titanmod_Cslip_3333(3,3,3,3,maxNinstance)) -constitutive_titanmod_CoverA = 0.0_pReal -constitutive_titanmod_C11 = 0.0_pReal -constitutive_titanmod_C12 = 0.0_pReal -constitutive_titanmod_C13 = 0.0_pReal -constitutive_titanmod_C33 = 0.0_pReal -constitutive_titanmod_C44 = 0.0_pReal -constitutive_titanmod_debyefrequency = 0.0_pReal -constitutive_titanmod_kinkf0 = 0.0_pReal -constitutive_titanmod_Gmod = 0.0_pReal -constitutive_titanmod_CAtomicVolume = 0.0_pReal -constitutive_titanmod_dc = 0.0_pReal -constitutive_titanmod_twinhpconstant = 0.0_pReal -constitutive_titanmod_GrainSize = 0.0_pReal -constitutive_titanmod_MaxTwinFraction = 0.0_pReal -constitutive_titanmod_r = 0.0_pReal -constitutive_titanmod_CEdgeDipMinDistance = 0.0_pReal -constitutive_titanmod_Cmfptwin = 0.0_pReal -constitutive_titanmod_Cthresholdtwin = 0.0_pReal -constitutive_titanmod_aTolRho = 0.0_pReal -constitutive_titanmod_Cslip_66 = 0.0_pReal -constitutive_titanmod_Cslip_3333 = 0.0_pReal + constitutive_titanmod_Cslip_3333 = 0.0_pReal allocate(constitutive_titanmod_rho_edge0(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_rho_edge0 = 0.0_pReal allocate(constitutive_titanmod_rho_screw0(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_rho_screw0 = 0.0_pReal allocate(constitutive_titanmod_shear_system0(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_shear_system0 = 0.0_pReal allocate(constitutive_titanmod_burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_burgersPerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_burgersPerTwinFamily(lattice_maxNtwinFamily,maxNinstance)) + constitutive_titanmod_burgersPerTwinFamily = 0.0_pReal allocate(constitutive_titanmod_f0_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_f0_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_tau0e_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_tau0e_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_tau0s_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_tau0s_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_capre_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_capre_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_caprs_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_caprs_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_pe_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_pe_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_ps_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_ps_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_qe_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_qe_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_qs_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_qs_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_v0e_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_v0e_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_v0s_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_v0s_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_kinkcriticallength_PerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_kinkcriticallength_PerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_twinsizePerTwinFamily(lattice_maxNtwinFamily,maxNinstance)) + constitutive_titanmod_twinsizePerTwinFamily = 0.0_pReal allocate(constitutive_titanmod_CeLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_CeLambdaSlipPerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_CsLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance)) + constitutive_titanmod_CsLambdaSlipPerSlipFamily = 0.0_pReal allocate(constitutive_titanmod_twinf0_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance)) + constitutive_titanmod_twinf0_PerTwinFamily = 0.0_pReal allocate(constitutive_titanmod_twinshearconstant_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance)) + constitutive_titanmod_twinshearconstant_PerTwinFamily = 0.0_pReal allocate(constitutive_titanmod_twintau0_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance)) + constitutive_titanmod_twintau0_PerTwinFamily = 0.0_pReal allocate(constitutive_titanmod_twinp_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance)) + constitutive_titanmod_twingamma0_PerTwinFamily = 0.0_pReal allocate(constitutive_titanmod_twinq_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance)) + constitutive_titanmod_twinLambdaSlipPerTwinFamily = 0.0_pReal allocate(constitutive_titanmod_twingamma0_PerTwinFamily(lattice_maxNTwinFamily,maxNinstance)) + constitutive_titanmod_twinp_PerTwinFamily = 0.0_pReal allocate(constitutive_titanmod_twinLambdaSlipPerTwinFamily(lattice_maxNTwinFamily,maxNinstance)) - -constitutive_titanmod_rho_edge0 = 0.0_pReal -constitutive_titanmod_rho_screw0 = 0.0_pReal -constitutive_titanmod_shear_system0 = 0.0_pReal -constitutive_titanmod_burgersPerSlipFamily = 0.0_pReal -constitutive_titanmod_burgersPerTwinFamily = 0.0_pReal -constitutive_titanmod_f0_PerSlipFamily = 0.0_pReal -constitutive_titanmod_tau0e_PerSlipFamily = 0.0_pReal -constitutive_titanmod_tau0s_PerSlipFamily = 0.0_pReal -constitutive_titanmod_capre_PerSlipFamily = 0.0_pReal -constitutive_titanmod_caprs_PerSlipFamily = 0.0_pReal -constitutive_titanmod_v0e_PerSlipFamily = 0.0_pReal -constitutive_titanmod_v0s_PerSlipFamily = 0.0_pReal -constitutive_titanmod_kinkcriticallength_PerSlipFamily = 0.0_pReal -constitutive_titanmod_twinsizePerTwinFamily = 0.0_pReal -constitutive_titanmod_CeLambdaSlipPerSlipFamily = 0.0_pReal -constitutive_titanmod_CsLambdaSlipPerSlipFamily = 0.0_pReal -constitutive_titanmod_pe_PerSlipFamily = 0.0_pReal -constitutive_titanmod_ps_PerSlipFamily = 0.0_pReal -constitutive_titanmod_qe_PerSlipFamily = 0.0_pReal -constitutive_titanmod_qs_PerSlipFamily = 0.0_pReal - -constitutive_titanmod_twinf0_PerTwinFamily = 0.0_pReal -constitutive_titanmod_twinshearconstant_PerTwinFamily = 0.0_pReal -constitutive_titanmod_twintau0_PerTwinFamily = 0.0_pReal -constitutive_titanmod_twingamma0_PerTwinFamily = 0.0_pReal -constitutive_titanmod_twinLambdaSlipPerTwinFamily = 0.0_pReal -constitutive_titanmod_twinp_PerTwinFamily = 0.0_pReal -constitutive_titanmod_twinq_PerTwinFamily = 0.0_pReal + constitutive_titanmod_twinq_PerTwinFamily = 0.0_pReal allocate(constitutive_titanmod_interactionSlipSlip(lattice_maxNinteraction,maxNinstance)) + constitutive_titanmod_interactionSlipSlip = 0.0_pReal allocate(constitutive_titanmod_interaction_ee(lattice_maxNinteraction,maxNinstance)) + constitutive_titanmod_interaction_ee = 0.0_pReal allocate(constitutive_titanmod_interaction_ss(lattice_maxNinteraction,maxNinstance)) + constitutive_titanmod_interaction_ss = 0.0_pReal allocate(constitutive_titanmod_interaction_es(lattice_maxNinteraction,maxNinstance)) + constitutive_titanmod_interaction_ss = 0.0_pReal allocate(constitutive_titanmod_interactionSlipTwin(lattice_maxNinteraction,maxNinstance)) + constitutive_titanmod_interactionSlipTwin = 0.0_pReal allocate(constitutive_titanmod_interactionTwinSlip(lattice_maxNinteraction,maxNinstance)) + constitutive_titanmod_interactionTwinSlip = 0.0_pReal allocate(constitutive_titanmod_interactionTwinTwin(lattice_maxNinteraction,maxNinstance)) -constitutive_titanmod_interactionSlipSlip = 0.0_pReal -constitutive_titanmod_interaction_ee = 0.0_pReal -constitutive_titanmod_interaction_ss = 0.0_pReal -constitutive_titanmod_interaction_ss = 0.0_pReal -constitutive_titanmod_interactionSlipTwin = 0.0_pReal -constitutive_titanmod_interactionTwinSlip = 0.0_pReal -constitutive_titanmod_interactionTwinTwin = 0.0_pReal + constitutive_titanmod_interactionTwinTwin = 0.0_pReal !* Readout data from material.config file rewind(file) line = '' -section = 0 +section = 0_pInt write(6,*) 'titanmod: Reading material parameters from material config file' @@ -447,7 +446,7 @@ enddo constitutive_titanmod_rho_edge0(j,i) = IO_floatValue(line,positions,1_pInt+j) write(6,*) tag,constitutive_titanmod_rho_edge0(1:4,i) case ('rho_screw0') - forall (j = 1:lattice_maxNslipFamily) & + forall (j = 1_pInt:lattice_maxNslipFamily) & constitutive_titanmod_rho_screw0(j,i) = IO_floatValue(line,positions,1_pInt+j) write(6,*) tag,constitutive_titanmod_rho_screw0(1:4,i) case ('slipburgers') @@ -606,7 +605,7 @@ write(6,*) 'Material Property reading done' if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt,e=i) if (sum(constitutive_titanmod_Nslip(:,i)) <= 0_pInt) call IO_error(231_pInt,e=i,ext_msg='nslip') if (sum(constitutive_titanmod_Ntwin(:,i)) < 0_pInt) call IO_error(231_pInt,e=i,ext_msg='ntwin') - do f = 1,lattice_maxNslipFamily + do f = 1_pInt,lattice_maxNslipFamily if (constitutive_titanmod_Nslip(f,i) > 0_pInt) then if (constitutive_titanmod_rho_edge0(f,i) < 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='rho_edge0') if (constitutive_titanmod_rho_screw0(f,i) < 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='rho_screw0') @@ -622,7 +621,7 @@ write(6,*) 'Material Property reading done' call IO_error(231_pInt,e=i,ext_msg='kinkCriticalLength') endif enddo - do f = 1,lattice_maxNtwinFamily + do f = 1_pInt,lattice_maxNtwinFamily if (constitutive_titanmod_Ntwin(f,i) > 0_pInt) then if (constitutive_titanmod_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinburgers') if (constitutive_titanmod_twinf0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(231_pInt,e=i,ext_msg='twinf0') @@ -765,7 +764,7 @@ do i = 1_pInt,maxNinstance write(6,*) 'Determined size of state and dot state' !* Determine size of postResults array - do o = 1,constitutive_titanmod_Noutput(i) + do o = 1_pInt,constitutive_titanmod_Noutput(i) select case(constitutive_titanmod_output(o,i)) case('rhoedge', & 'rhoscrew', & @@ -849,20 +848,20 @@ write(6,*) 'Determining elasticity matrix' do k=1_pInt,constitutive_titanmod_Ntwin(j,i) do l=1_pInt,3_pInt ; do m=1_pInt,3_pInt ; do n=1_pInt,3_pInt ; do o=1_pInt,3_pInt do p=1_pInt,3_pInt ; do q=1_pInt,3_pInt ; do r=1_pInt,3_pInt ; do s=1_pInt,3_pInt - constitutive_titanmod_Ctwin_3333(l,m,n,o,sum(constitutive_titanmod_Nslip(1:j-1,i))+k,i) = & - constitutive_titanmod_Ctwin_3333(l,m,n,o,sum(constitutive_titanmod_Nslip(1:j-1,i))+k,i) + & + constitutive_titanmod_Ctwin_3333(l,m,n,o,sum(constitutive_titanmod_Nslip(1:j-1_pInt,i))+k,i) = & + constitutive_titanmod_Ctwin_3333(l,m,n,o,sum(constitutive_titanmod_Nslip(1:j-1_pInt,i))+k,i) + & constitutive_titanmod_Cslip_3333(p,q,r,s,i)*& - lattice_Qtwin(l,p,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* & - lattice_Qtwin(m,q,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* & - lattice_Qtwin(n,r,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* & - lattice_Qtwin(o,s,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure) + lattice_Qtwin(l,p,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure)* & + lattice_Qtwin(m,q,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure)* & + lattice_Qtwin(n,r,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure)* & + lattice_Qtwin(o,s,sum(lattice_NslipSystem(1:j-1_pInt,myStructure))+k,myStructure) enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo constitutive_titanmod_Ctwin_66(:,:,k,i) = math_Mandel3333to66(constitutive_titanmod_Ctwin_3333(:,:,:,:,k,i)) enddo enddo !* Burgers vector, dislocation velocity prefactor for each slip system - do s = 1,constitutive_titanmod_totalNslip(i) + do s = 1_pInt,constitutive_titanmod_totalNslip(i) f = constitutive_titanmod_slipFamily(s,i) constitutive_titanmod_burgersPerSlipSystem(s,i) = constitutive_titanmod_burgersPerSlipFamily(f,i) constitutive_titanmod_f0_PerSlipSystem(s,i) = constitutive_titanmod_f0_PerSlipFamily(f,i) @@ -882,7 +881,7 @@ write(6,*) 'Determining elasticity matrix' enddo !* Burgers vector, nucleation rate prefactor and twin size for each twin system - do t = 1,constitutive_titanmod_totalNtwin(i) + do t = 1_pInt,constitutive_titanmod_totalNtwin(i) f = constitutive_titanmod_twinFamily(t,i) constitutive_titanmod_burgersPerTwinSystem(t,i) = constitutive_titanmod_burgersPerTwinFamily(f,i) constitutive_titanmod_twinsizePerTwinSystem(t,i) = constitutive_titanmod_twinsizePerTwinFamily(f,i) @@ -897,56 +896,44 @@ write(6,*) 'Determining elasticity matrix' enddo !* Construction of interaction matrices - do s1 = 1,constitutive_titanmod_totalNslip(i) - do s2 = 1,constitutive_titanmod_totalNslip(i) + do s1 = 1_pInt,constitutive_titanmod_totalNslip(i) + do s2 = 1_pInt,constitutive_titanmod_totalNslip(i) constitutive_titanmod_interactionMatrixSlipSlip(s1,s2,i) = & constitutive_titanmod_interactionSlipSlip(lattice_interactionSlipSlip(constitutive_titanmod_slipSystemLattice(s1,i), & constitutive_titanmod_slipSystemLattice(s2,i), & myStructure),i) - enddo; enddo - - do s1 = 1,constitutive_titanmod_totalNslip(i) - do s2 = 1,constitutive_titanmod_totalNslip(i) constitutive_titanmod_interactionMatrix_ee(s1,s2,i) = & constitutive_titanmod_interaction_ee(lattice_interactionSlipSlip(constitutive_titanmod_slipSystemLattice(s1,i), & constitutive_titanmod_slipSystemLattice(s2,i), & myStructure),i) - enddo; enddo - - do s1 = 1,constitutive_titanmod_totalNslip(i) - do s2 = 1,constitutive_titanmod_totalNslip(i) constitutive_titanmod_interactionMatrix_ss(s1,s2,i) = & constitutive_titanmod_interaction_ss(lattice_interactionSlipSlip(constitutive_titanmod_slipSystemLattice(s1,i), & constitutive_titanmod_slipSystemLattice(s2,i), & myStructure),i) - enddo; enddo - - do s1 = 1,constitutive_titanmod_totalNslip(i) - do s2 = 1,constitutive_titanmod_totalNslip(i) constitutive_titanmod_interactionMatrix_es(s1,s2,i) = & constitutive_titanmod_interaction_es(lattice_interactionSlipSlip(constitutive_titanmod_slipSystemLattice(s1,i), & constitutive_titanmod_slipSystemLattice(s2,i), & myStructure),i) enddo; enddo - do s1 = 1,constitutive_titanmod_totalNslip(i) - do t2 = 1,constitutive_titanmod_totalNtwin(i) + do s1 = 1_pInt,constitutive_titanmod_totalNslip(i) + do t2 = 1_pInt,constitutive_titanmod_totalNtwin(i) constitutive_titanmod_interactionMatrixSlipTwin(s1,t2,i) = & constitutive_titanmod_interactionSlipTwin(lattice_interactionSlipTwin(constitutive_titanmod_slipSystemLattice(s1,i), & constitutive_titanmod_twinSystemLattice(t2,i), & myStructure),i) enddo; enddo - do t1 = 1,constitutive_titanmod_totalNtwin(i) - do s2 = 1,constitutive_titanmod_totalNslip(i) + do t1 = 1_pInt,constitutive_titanmod_totalNtwin(i) + do s2 = 1_pInt,constitutive_titanmod_totalNslip(i) constitutive_titanmod_interactionMatrixTwinSlip(t1,s2,i) = & constitutive_titanmod_interactionTwinSlip(lattice_interactionTwinSlip(constitutive_titanmod_twinSystemLattice(t1,i), & constitutive_titanmod_slipSystemLattice(s2,i), & myStructure),i) enddo; enddo - do t1 = 1,constitutive_titanmod_totalNtwin(i) - do t2 = 1,constitutive_titanmod_totalNtwin(i) + do t1 = 1_pInt,constitutive_titanmod_totalNtwin(i) + do t2 = 1_pInt,constitutive_titanmod_totalNtwin(i) constitutive_titanmod_interactionMatrixTwinTwin(t1,t2,i) = & constitutive_titanmod_interactionTwinTwin(lattice_interactionTwinTwin(constitutive_titanmod_twinSystemLattice(t1,i), & constitutive_titanmod_twinSystemLattice(t2,i), & @@ -954,16 +941,12 @@ write(6,*) 'Determining elasticity matrix' enddo; enddo !* Calculation of forest projections for edge dislocations - do s1 = 1,constitutive_titanmod_totalNslip(i) - do s2 = 1,constitutive_titanmod_totalNslip(i) + do s1 = 1_pInt,constitutive_titanmod_totalNslip(i) + do s2 = 1_pInt,constitutive_titanmod_totalNslip(i) constitutive_titanmod_forestProjectionEdge(s1,s2,i) = & abs(math_mul3x3(lattice_sn(:,constitutive_titanmod_slipSystemLattice(s1,i),myStructure), & lattice_st(:,constitutive_titanmod_slipSystemLattice(s2,i),myStructure))) - enddo; enddo - !* Calculation of forest projections for screw dislocations - do s1 = 1,constitutive_titanmod_totalNslip(i) - do s2 = 1,constitutive_titanmod_totalNslip(i) constitutive_titanmod_forestProjectionScrew(s1,s2,i) = & abs(math_mul3x3(lattice_sn(:,constitutive_titanmod_slipSystemLattice(s1,i),myStructure), & lattice_sd(:,constitutive_titanmod_slipSystemLattice(s2,i),myStructure))) @@ -971,22 +954,17 @@ write(6,*) 'Determining elasticity matrix' !* Calculation of forest projections for edge dislocations in twin system - do t1 = 1,constitutive_titanmod_totalNtwin(i) - do t2 = 1,constitutive_titanmod_totalNtwin(i) + do t1 = 1_pInt,constitutive_titanmod_totalNtwin(i) + do t2 = 1_pInt,constitutive_titanmod_totalNtwin(i) constitutive_titanmod_TwinforestProjectionEdge(t1,t2,i) = & abs(math_mul3x3(lattice_tn(:,constitutive_titanmod_twinSystemLattice(t1,i),myStructure), & lattice_tt(:,constitutive_titanmod_twinSystemLattice(t2,i),myStructure))) - enddo; enddo - !* Calculation of forest projections for screw dislocations in twin system - do t1 = 1,constitutive_titanmod_totalNtwin(i) - do t2 = 1,constitutive_titanmod_totalNtwin(i) constitutive_titanmod_TwinforestProjectionScrew(t1,t2,i) = & abs(math_mul3x3(lattice_tn(:,constitutive_titanmod_twinSystemLattice(t1,i),myStructure), & lattice_td(:,constitutive_titanmod_twinSystemLattice(t2,i),myStructure))) enddo; enddo - enddo write(6,*) 'Init All done' return @@ -998,7 +976,6 @@ function constitutive_titanmod_stateInit(myInstance) !* initial microstructural state * !********************************************************************* use prec, only: pReal,pInt -use math, only: pi use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily implicit none @@ -1024,7 +1001,7 @@ constitutive_titanmod_stateInit = 0.0_pReal !* Initialize basic slip state variables ! For slip s1 = 0_pInt -do f = 1,lattice_maxNslipFamily +do f = 1_pInt,lattice_maxNslipFamily s0 = s1 + 1_pInt s1 = s0 + constitutive_titanmod_Nslip(f,myInstance) - 1_pInt do s = s0,s1 @@ -1037,7 +1014,7 @@ enddo !* Initialize basic slip state variables ! For twin ts1 = 0_pInt -do tf = 1,lattice_maxNtwinFamily +do tf = 1_pInt,lattice_maxNtwinFamily ts0 = ts1 + 1_pInt ts1 = ts0 + constitutive_titanmod_Ntwin(tf,myInstance) - 1_pInt do ts = ts0,ts1 @@ -1046,44 +1023,44 @@ do tf = 1,lattice_maxNtwinFamily enddo constitutive_titanmod_stateInit(1:ns) = rho_edge0 -constitutive_titanmod_stateInit(ns+1:2*ns) = rho_screw0 -constitutive_titanmod_stateInit(2*ns+1:3*ns) = shear_system0 -constitutive_titanmod_stateInit(3*ns+1:3*ns+nt) = twingamma_dot0 +constitutive_titanmod_stateInit(ns+1_pInt:2_pInt*ns) = rho_screw0 +constitutive_titanmod_stateInit(2_pInt*ns+1_pInt:3_pInt*ns) = shear_system0 +constitutive_titanmod_stateInit(3_pInt*ns+1_pInt:3_pInt*ns+nt) = twingamma_dot0 !* Initialize dependent slip microstructural variables -forall (s = 1:ns) & +forall (s = 1_pInt:ns) & segment_edge0(s) = constitutive_titanmod_CeLambdaSlipPerSlipSystem(s,myInstance)/ & sqrt(dot_product((rho_edge0),constitutive_titanmod_forestProjectionEdge(1:ns,s,myInstance))+ & dot_product((rho_screw0),constitutive_titanmod_forestProjectionScrew(1:ns,s,myInstance))) -constitutive_titanmod_stateInit(3*ns+nt+1:4*ns+nt) = segment_edge0 +constitutive_titanmod_stateInit(3_pInt*ns+nt+1_pInt:4_pInt*ns+nt) = segment_edge0 -forall (s = 1:ns) & +forall (s = 1_pInt:ns) & segment_screw0(s) = constitutive_titanmod_CsLambdaSlipPerSlipSystem(s,myInstance)/ & sqrt(dot_product((rho_edge0),constitutive_titanmod_forestProjectionEdge(1:ns,s,myInstance))+ & dot_product((rho_screw0),constitutive_titanmod_forestProjectionScrew(1:ns,s,myInstance))) -constitutive_titanmod_stateInit(4*ns+nt+1:5*ns+nt) = segment_screw0 +constitutive_titanmod_stateInit(4_pInt*ns+nt+1_pInt:5_pInt*ns+nt) = segment_screw0 -forall (s = 1:ns) & +forall (s = 1_pInt:ns) & resistance_edge0(s) = & constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerSlipSystem(s,myInstance)* & sqrt(dot_product((rho_edge0),constitutive_titanmod_interactionMatrix_ee(1:ns,s,myInstance))+dot_product((rho_screw0), & constitutive_titanmod_interactionMatrix_es(1:ns,s,myInstance))) -constitutive_titanmod_stateInit(5*ns+nt+1:6*ns+nt) = resistance_edge0 +constitutive_titanmod_stateInit(5_pInt*ns+nt+1_pInt:6_pInt*ns+nt) = resistance_edge0 -forall (s = 1:ns) & +forall (s = 1_pInt:ns) & resistance_screw0(s) = & constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerSlipSystem(s,myInstance)* & sqrt(dot_product((rho_edge0),constitutive_titanmod_interactionMatrix_es(1:ns,s,myInstance))+dot_product((rho_screw0), & constitutive_titanmod_interactionMatrix_ss(1:ns,s,myInstance))) -constitutive_titanmod_stateInit(6*ns+nt+1:7*ns+nt) = resistance_screw0 +constitutive_titanmod_stateInit(6_pInt*ns+nt+1_pInt:7_pInt*ns+nt) = resistance_screw0 -forall (t = 1:nt) & +forall (t = 1_pInt:nt) & resistance_twin0(t) = 0.0_pReal -constitutive_titanmod_stateInit(7*ns+nt+1:7*ns+2*nt)=resistance_twin0 +constitutive_titanmod_stateInit(7_pInt*ns+nt+1_pInt:7_pInt*ns+2_pInt*nt)=resistance_twin0 return end function @@ -1133,8 +1110,8 @@ ns = constitutive_titanmod_totalNslip(myInstance) nt = constitutive_titanmod_totalNtwin(myInstance) !* Total twin volume fraction -do i=1,nt -volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ & +do i=1_pInt,nt +volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance) enddo !sumf = sum(state(g,ip,el)%p((6*ns+7*nt+1):(6*ns+8*nt))) ! safe for nt == 0 @@ -1142,7 +1119,7 @@ sumf = sum(abs(volumefraction_pertwinsystem(1:nt))) ! safe for nt == 0 !* Homogenized elasticity matrix constitutive_titanmod_homogenizedC = (1.0_pReal-sumf)*constitutive_titanmod_Cslip_66(:,:,myInstance) -do i=1,nt +do i=1_pInt,nt constitutive_titanmod_homogenizedC = & ! constitutive_titanmod_homogenizedC + state(g,ip,el)%p(6*ns+7*nt+i)*constitutive_titanmod_Ctwin_66(:,:,i,myInstance) constitutive_titanmod_homogenizedC + volumefraction_pertwinsystem(i)*constitutive_titanmod_Ctwin_66(:,:,i,myInstance) @@ -1163,10 +1140,8 @@ subroutine constitutive_titanmod_microstructure(Temperature,state,g,ip,el) !* - el : current element * !********************************************************************* use prec, only: pReal,pInt,p_vec -use math, only: pi use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance -use lattice, only: lattice_interactionSlipTwin,lattice_interactionTwinTwin !use debug, only: debugger implicit none @@ -1205,8 +1180,8 @@ nt = constitutive_titanmod_totalNtwin(myInstance) !* State: 13*ns+2*nt+1 : 14*ns+2*nt StressRatio_screw_p !* Total twin volume fraction -do i=1,nt -volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ & +do i=1_pInt,nt +volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance) enddo @@ -1224,45 +1199,45 @@ sfe = 0.0002_pReal*Temperature-0.0396_pReal ! state(g,ip,el)%p(2*ns+t)/constitutive_titanmod_twinsizePerTwinSystem(t,myInstance) ! average segment length for edge dislocations in matrix -forall (s = 1:ns) & - state(g,ip,el)%p(3*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSystem(s,myInstance)/ & +forall (s = 1_pInt:ns) & + state(g,ip,el)%p(3_pInt*ns+nt+s) = constitutive_titanmod_CeLambdaSlipPerSlipSystem(s,myInstance)/ & sqrt(dot_product(state(g,ip,el)%p(1:ns), & constitutive_titanmod_forestProjectionEdge(1:ns,s,myInstance))+ & - dot_product(state(g,ip,el)%p(ns+1:2*ns), & + dot_product(state(g,ip,el)%p(ns+1_pInt:2_pInt*ns), & constitutive_titanmod_forestProjectionScrew(1:ns,s,myInstance))) ! average segment length for screw dislocations in matrix -forall (s = 1:ns) & - state(g,ip,el)%p(4*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSystem(s,myInstance)/ & +forall (s = 1_pInt:ns) & + state(g,ip,el)%p(4_pInt*ns+nt+s) = constitutive_titanmod_CsLambdaSlipPerSlipSystem(s,myInstance)/ & sqrt(dot_product(state(g,ip,el)%p(1:ns), & constitutive_titanmod_forestProjectionEdge(1:ns,s,myInstance))+ & - dot_product(state(g,ip,el)%p(ns+1:2*ns), & + dot_product(state(g,ip,el)%p(ns+1_pInt:2_pInt*ns), & constitutive_titanmod_forestProjectionScrew(1:ns,s,myInstance))) !* threshold stress or slip resistance for edge dislocation motion -forall (s = 1:ns) & - state(g,ip,el)%p(5*ns+nt+s) = & +forall (s = 1_pInt:ns) & + state(g,ip,el)%p(5_pInt*ns+nt+s) = & constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerSlipSystem(s,myInstance)*& sqrt(dot_product((state(g,ip,el)%p(1:ns)),& constitutive_titanmod_interactionMatrix_ee(1:ns,s,myInstance))+ & - dot_product((state(g,ip,el)%p(ns+1:2*ns)),& + dot_product((state(g,ip,el)%p(ns+1_pInt:2_pInt*ns)),& constitutive_titanmod_interactionMatrix_es(1:ns,s,myInstance))) !* threshold stress or slip resistance for screw dislocation motion -forall (s = 1:ns) & - state(g,ip,el)%p(6*ns+nt+s) = & +forall (s = 1_pInt:ns) & + state(g,ip,el)%p(6_pInt*ns+nt+s) = & constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerSlipSystem(s,myInstance)*& sqrt(dot_product((state(g,ip,el)%p(1:ns)),& constitutive_titanmod_interactionMatrix_es(1:ns,s,myInstance))+ & - dot_product((state(g,ip,el)%p(ns+1:2*ns)),& + dot_product((state(g,ip,el)%p(ns+1_pInt:2_pInt*ns)),& constitutive_titanmod_interactionMatrix_ss(1:ns,s,myInstance))) !* threshold stress or slip resistance for dislocation motion in twin -forall (t = 1:nt) & - state(g,ip,el)%p(7*ns+nt+t) = & +forall (t = 1_pInt:nt) & + state(g,ip,el)%p(7_pInt*ns+nt+t) = & constitutive_titanmod_Gmod(myInstance)*constitutive_titanmod_burgersPerTwinSystem(t,myInstance)*& - (dot_product((abs(state(g,ip,el)%p(2*ns+1:2*ns+nt))),& + (dot_product((abs(state(g,ip,el)%p(2_pInt*ns+1_pInt:2_pInt*ns+nt))),& constitutive_titanmod_interactionMatrixTwinTwin(1:nt,t,myInstance))) @@ -1288,8 +1263,8 @@ use prec, only: pReal,pInt,p_vec use math, only: math_Plain3333to99 use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance -use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin +use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & + lattice_NslipSystem,lattice_NtwinSystem, lattice_Stwin implicit none !* Input-Output variables @@ -1317,8 +1292,8 @@ myStructure = constitutive_titanmod_structure(myInstance) ns = constitutive_titanmod_totalNslip(myInstance) nt = constitutive_titanmod_totalNtwin(myInstance) -do i=1,nt -volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ & +do i=1_pInt,nt +volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance) enddo @@ -1336,9 +1311,9 @@ gdot_slip_edge = 0.0_pReal gdot_slip_screw = 0.0_pReal dgdot_dtauslip = 0.0_pReal j = 0_pInt -do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_titanmod_Nslip(f,myInstance) ! process each (active) slip system in family +do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_titanmod_Nslip(f,myInstance) ! process each (active) slip system in family j = j+1_pInt !* Calculation of Lp @@ -1347,9 +1322,9 @@ do f = 1,lattice_maxNslipFamily ! loop over all !************************************************* !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! if(myStructure>=3.and.j>3) then ! for all non-basal slip systems - if(myStructure==3) then ! only for prismatic and pyr systems in hex + if(myStructure==3_pInt) then ! only for prismatic and pyr systems in hex screwvelocity_prefactor=constitutive_titanmod_debyefrequency(myInstance)* & - state(g,ip,el)%p(4*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSystem(j,myInstance)/ & + state(g,ip,el)%p(4_pInt*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSystem(j,myInstance)/ & constitutive_titanmod_kinkcriticallength_PerSlipSystem(j,myInstance))**2 !* Stress ratio for screw ! No slip resistance for screw dislocations, only Peierls stress @@ -1472,7 +1447,7 @@ do f = 1,lattice_maxNslipFamily ! loop over all Lp = Lp + (1.0_pReal - sumf)*gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,myStructure) !* Calculation of the tangent of Lp - forall (k=1:3,l=1:3,m=1:3,n=1:3) & + 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(j)*& lattice_Sslip(k,l,index_myFamily+i,myStructure)*& @@ -1484,9 +1459,9 @@ enddo gdot_twin = 0.0_pReal dgdot_dtautwin = 0.0_pReal j = 0_pInt -do f = 1,lattice_maxNtwinFamily ! loop over all slip families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_titanmod_Ntwin(f,myInstance) ! process each (active) slip system in family +do f = 1_pInt,lattice_maxNtwinFamily ! loop over all slip families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_titanmod_Ntwin(f,myInstance) ! process each (active) slip system in family j = j+1_pInt !* Calculation of Lp @@ -1553,7 +1528,7 @@ do f = 1,lattice_maxNtwinFamily ! loop over all Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,myStructure) !* Calculation of the tangent of Lp - forall (k=1:3,l=1:3,m=1:3,n=1:3) & + 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(j)*& lattice_Stwin(k,l,index_myFamily+i,myStructure)*& @@ -1592,11 +1567,10 @@ function constitutive_titanmod_dotState(Tstar_v,Temperature,state,g,ip,el) !********************************************************************* use prec, only: pReal,pInt,p_vec -use math, only: pi use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance -use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin +use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily, & + lattice_NslipSystem,lattice_NtwinSystem, lattice_Stwin_v implicit none !* Input-Output variables @@ -1623,20 +1597,20 @@ MyStructure = constitutive_titanmod_structure(myInstance) ns = constitutive_titanmod_totalNslip(myInstance) nt = constitutive_titanmod_totalNtwin(myInstance) -do i=1,nt -volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ & +do i=1_pInt,nt +volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance) enddo -sumf = sum(abs(volumefraction_pertwinsystem(1:nt))) ! safe for nt == 0 +sumf = sum(abs(volumefraction_pertwinsystem(1_pInt:nt))) ! safe for nt == 0 constitutive_titanmod_dotState = 0.0_pReal j = 0_pInt - do f = 1,lattice_maxNslipFamily ! loop over all slip families - index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_titanmod_Nslip(f,myInstance) ! process each (active) slip system in family + do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_titanmod_Nslip(f,myInstance) ! process each (active) slip system in family j = j+1_pInt !* Multiplication of edge dislocations @@ -1667,9 +1641,9 @@ constitutive_titanmod_dotState = 0.0_pReal !* Twin fraction evolution j = 0_pInt -do f = 1,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,MyStructure)) ! at which index starts my family - do i = 1,constitutive_titanmod_Ntwin(f,myInstance) ! process each (active) twin system in family +do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,MyStructure)) ! at which index starts my family + do i = 1_pInt,constitutive_titanmod_Ntwin(f,myInstance) ! process each (active) twin system in family j = j+1_pInt !************************************************************************* @@ -1769,11 +1743,8 @@ pure function constitutive_titanmod_postResults(Tstar_v,Temperature,dt,state,g,i !* - el : current element * !********************************************************************* use prec, only: pReal,pInt,p_vec -use math, only: pi use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput -use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin implicit none !* Definition of variables @@ -1794,8 +1765,8 @@ myStructure = constitutive_titanmod_structure(myInstance) ns = constitutive_titanmod_totalNslip(myInstance) nt = constitutive_titanmod_totalNtwin(myInstance) -do i=1,nt -volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3*ns+i)/ & +do i=1_pInt,nt +volumefraction_pertwinsystem(i)=state(g,ip,el)%p(3_pInt*ns+i)/ & constitutive_titanmod_twinshearconstant_PerTwinSystem(i,myInstance) enddo @@ -1809,99 +1780,99 @@ sumf = sum(abs(volumefraction_pertwinsystem(1:nt))) ! safe for nt == 0 c = 0_pInt constitutive_titanmod_postResults = 0.0_pReal -do o = 1,phase_Noutput(material_phase(g,ip,el)) +do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) select case(constitutive_titanmod_output(o,myInstance)) case ('rhoedge') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p(1_pInt:ns) c = c + ns case ('rhoscrew') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p(ns+1:2*ns) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p(ns+1_pInt:2_pInt*ns) c = c + ns case ('segment_edge') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((3_pInt*ns+nt+1_pInt):(4_pInt*ns+nt)) c = c + ns case ('segment_screw') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((4*ns+nt+1):(5*ns+nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((4_pInt*ns+nt+1_pInt):(5_pInt*ns+nt)) c = c + ns case ('resistance_edge') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((5*ns+nt+1):(6*ns+nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((5_pInt*ns+nt+1_pInt):(6_pInt*ns+nt)) c = c + ns case ('resistance_screw') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((6*ns+nt+1):(7*ns+nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((6_pInt*ns+nt+1_pInt):(7_pInt*ns+nt)) c = c + ns case ('velocity_edge') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((7*ns+2*nt+1):(8*ns+2*nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((7*ns+2*nt+1):(8*ns+2*nt)) c = c + ns case ('velocity_screw') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((8*ns+2*nt+1):(9*ns+2*nt)) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = state(g,ip,el)%p((8*ns+2*nt+1):(9*ns+2*nt)) c = c + ns case ('tau_slip') - constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((9*ns+2*nt+1):(10*ns+2*nt))) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((9*ns+2*nt+1):(10*ns+2*nt))) c = c + ns case ('gdot_slip_edge') - constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((10*ns+2*nt+1):(11*ns+2*nt))) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((10*ns+2*nt+1):(11*ns+2*nt))) c = c + ns case ('gdot_slip_screw') - constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((11*ns+2*nt+1):(12*ns+2*nt))) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((11*ns+2*nt+1):(12*ns+2*nt))) c = c + ns case ('gdot_slip') - constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((10*ns+2*nt+1):(11*ns+2*nt))) + & + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((10*ns+2*nt+1):(11*ns+2*nt))) + & abs(state(g,ip,el)%p((11*ns+2*nt+1):(12*ns+2*nt))) c = c + ns case ('stressratio_edge_p') - constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((12*ns+2*nt+1):(13*ns+2*nt))) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((12*ns+2*nt+1):(13*ns+2*nt))) c = c + ns case ('stressratio_screw_p') - constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((13*ns+2*nt+1):(14*ns+2*nt))) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((13*ns+2*nt+1):(14*ns+2*nt))) c = c + ns case ('shear_system') - constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((2*ns+1):(3*ns))) + constitutive_titanmod_postResults(c+1_pInt:c+ns) = abs(state(g,ip,el)%p((2*ns+1):(3*ns))) c = c + ns case ('shear_basal') - constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+1):(2*ns+3)))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+1):(2*ns+3)))) + c = c + 1_pInt case ('shear_prism') - constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+4):(2*ns+6)))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+4):(2*ns+6)))) + c = c + 1_pInt case ('shear_pyra') - constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+7):(2*ns+12)))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+7):(2*ns+12)))) + c = c + 1_pInt case ('shear_pyrca') - constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+13):(2*ns+24)))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+13):(2*ns+24)))) + c = c + 1_pInt case ('rhoedge_basal') - constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((1):(3))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((1):(3))) + c = c + 1_pInt case ('rhoedge_prism') - constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((4):(6))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((4):(6))) + c = c + 1_pInt case ('rhoedge_pyra') - constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((7):(12))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((7):(12))) + c = c + 1_pInt case ('rhoedge_pyrca') - constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((13):(24))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((13):(24))) + c = c + 1_pInt case ('rhoscrew_basal') - constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((ns+1):(ns+3))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((ns+1):(ns+3))) + c = c + 1_pInt case ('rhoscrew_prism') - constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((ns+4):(ns+6))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((ns+4):(ns+6))) + c = c + 1_pInt case ('rhoscrew_pyra') - constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((ns+7):(ns+12))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((ns+7):(ns+12))) + c = c + 1_pInt case ('rhoscrew_pyrca') - constitutive_titanmod_postResults(c+1:c+1) = sum(state(g,ip,el)%p((ns+13):(ns+24))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(state(g,ip,el)%p((ns+13):(ns+24))) + c = c + 1_pInt case ('shear_total') - constitutive_titanmod_postResults(c+1:c+1) = sum(abs(state(g,ip,el)%p((2*ns+1):(3*ns)))) - c = c + 1 + constitutive_titanmod_postResults(c+1_pInt:c+1_pInt) = sum(abs(state(g,ip,el)%p((2*ns+1):(3*ns)))) + c = c + 1_pInt case ('twin_fraction') - constitutive_titanmod_postResults(c+1:c+nt) = abs(volumefraction_pertwinsystem(1:nt)) + constitutive_titanmod_postResults(c+1_pInt:c+nt) = abs(volumefraction_pertwinsystem(1:nt)) c = c + nt ! 'rhoedge', & ! 'rhoscrew', &