223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

again: adding pInt, removing redundant use statments, chang in dble to real(,pReal)

This commit is contained in:
Martin Diehl 2012-02-21 16:31:37 +00:00
parent 9dc730dea4
commit d8ffc29236
5 changed files with 300 additions and 317 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
code/CPFEM.f90
View File

@ -109,7 +109,8 @@ end subroutine
!*********************************************************
subroutine CPFEM_init()
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 IO, only: IO_read_jobBinaryFile
@ -128,9 +129,7 @@ subroutine CPFEM_init()
crystallite_dPdF0, &
crystallite_Tstar0_v
use homogenization, only: homogenization_sizeState, &
homogenization_state0, &
materialpoint_F, &
materialpoint_F0
homogenization_state0
implicit none
@ -231,12 +230,10 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
!*** variables and functions from other modules ***!
use prec, only: pReal, &
pInt
use numerics, only: relevantStrain, &
defgradTolerance, &
use numerics, only: defgradTolerance, &
iJacoStiffness
use debug, only: debug_e, &
debug_i, &
debug_g, &
debug_selectiveDebugger, &
debug_verbosity, &
debug_stressMaxLocation, &
@ -282,8 +279,8 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
microstructure_elemhomo, &
material_phase
use constitutive, only: constitutive_state0,constitutive_state
use crystallite, only: crystallite_F0, &
crystallite_partionedF, &
use crystallite, only: crystallite_partionedF,&
crystallite_F0, &
crystallite_Fp0, &
crystallite_Fp, &
crystallite_Lp0, &
@ -291,8 +288,7 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
crystallite_dPdF0, &
crystallite_dPdF, &
crystallite_Tstar0_v, &
crystallite_Tstar_v, &
crystallite_localConstitution
crystallite_Tstar_v
use homogenization, only: homogenization_sizeState, &
homogenization_state, &
homogenization_state0, &

164
code/crystallite.f90
View File

@ -103,14 +103,12 @@ CONTAINS
subroutine crystallite_init(Temperature)
!*** variables and functions from other modules ***!
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pInt, &
pReal
use debug, only: debug_info, &
debug_reset, &
debug_verbosity
use numerics, only: subStepSizeCryst, &
stepIncreaseCryst
use math, only: math_I3, &
math_EulerToR, &
math_inv33, &
@ -125,14 +123,11 @@ use mesh, only: mesh_element, &
mesh_maxNipNeighbors
use IO
use material
use lattice, only: lattice_symmetryType, &
lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem
use lattice, only: lattice_symmetryType
use constitutive, only: constitutive_microstructure
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, &
constitutive_phenopowerlaw_structure, &
constitutive_phenopowerlaw_Nslip
constitutive_phenopowerlaw_structure
use constitutive_titanmod, only: constitutive_titanmod_label, &
constitutive_titanmod_structure
use constitutive_dislotwin, only: constitutive_dislotwin_label, &
@ -141,8 +136,8 @@ use constitutive_nonlocal, only: constitutive_nonlocal_label, &
constitutive_nonlocal_structure
implicit none
integer(pInt), parameter :: file = 200, &
maxNchunks = 2
integer(pInt), parameter :: myFile = 200_pInt, &
maxNchunks = 2_pInt
!*** input variables ***!
real(pReal) Temperature
@ -234,18 +229,18 @@ allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), &
material_Ncrystallite)) ; crystallite_sizePostResult = 0_pInt
if (.not. IO_open_jobFile_stat(file,material_localFileExt)) then ! no local material configuration present...
call IO_open_file(file,material_configFile) ! ...open material.config file
if (.not. IO_open_jobFile_stat(myFile,material_localFileExt)) then ! no local material configuration present...
call IO_open_file(myFile,material_configFile) ! ...open material.config file
endif
line = ''
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= material_partCrystallite) ! wind forward to <crystallite>
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
@ -263,7 +258,7 @@ do ! read thru sections of
endif
enddo
100 close(file)
100 close(myFile)
do i = 1_pInt,material_Ncrystallite ! sanity checks
enddo
@ -299,18 +294,18 @@ enddo
! write description file for crystallite output
call IO_write_jobFile(file,'outputCrystallite')
call IO_write_jobFile(myFile,'outputCrystallite')
do p = 1_pInt,material_Ncrystallite
write(file,*)
write(file,'(a)') '['//trim(crystallite_name(p))//']'
write(file,*)
do e = 1,crystallite_Noutput(p)
write(file,'(a,i4)') trim(crystallite_output(e,p))//char(9),crystallite_sizePostResult(e,p)
write(myFile,*)
write(myFile,'(a)') '['//trim(crystallite_name(p))//']'
write(myFile,*)
do e = 1_pInt,crystallite_Noutput(p)
write(myFile,'(a,i4)') trim(crystallite_output(e,p))//char(9),crystallite_sizePostResult(e,p)
enddo
enddo
close(file)
close(myFile)
!$OMP PARALLEL PRIVATE(myNgrains,myPhase,myMat,myStructure)
@ -318,7 +313,7 @@ close(file)
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements
myNgrains = homogenization_Ngrains(mesh_element(3,e)) ! look up homogenization-->grainCount
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element
do g = 1,myNgrains
do g = 1_pInt,myNgrains
crystallite_Fp0(1:3,1:3,g,i,e) = math_EulerToR(material_EulerAngles(1:3,g,i,e)) ! plastic def gradient reflects init orientation
crystallite_F0(1:3,1:3,g,i,e) = math_I3
crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e))
@ -342,7 +337,7 @@ crystallite_requested = .true.
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,myNgrains
do g = 1_pInt,myNgrains
myPhase = material_phase(g,i,e)
myMat = phase_constitutionInstance(myPhase)
select case (phase_constitution(myPhase))
@ -374,7 +369,7 @@ crystallite_orientation0 = crystallite_orientation ! Store initial o
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,myNgrains
do g = 1_pInt,myNgrains
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), g, i, e) ! update dependent state variables to be consistent with basic states
enddo
@ -459,7 +454,6 @@ use numerics, only: subStepMinCryst, &
pert_Fg, &
pert_method, &
nCryst, &
iJacoStiffness, &
numerics_integrator, &
numerics_integrationMode
use debug, only: debug_verbosity, &
@ -557,7 +551,7 @@ crystallite_subStep = 0.0_pReal
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
do g = 1_pInt,myNgrains
if (crystallite_requested(g,i,e)) then ! initialize restoration point of ...
crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature
constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure
@ -710,7 +704,7 @@ enddo
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write (6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write (6,*)
write (6,'(a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e)) / 1e6
write (6,'(a,/,3(12x,3(f12.4,1x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal
write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Fp', math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write (6,'(a,/,3(12x,3(f14.9,1x)/))') '<< CRYST >> Lp', math_transpose33(crystallite_Lp(1:3,1:3,g,i,e))
write (6,*)
@ -922,10 +916,10 @@ use prec, only: pInt, &
pReal
use numerics, only: numerics_integrationMode
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_e, &
debug_i, &
debug_g, &
debug_selectiveDebugger, &
debug_StateLoopDistribution
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
@ -1029,7 +1023,7 @@ RK4dotTemperature = 0.0_pReal
! --- SECOND TO FOURTH RUNGE KUTTA STEP PLUS FINAL INTEGRATION ---
do n = 1,4
do n = 1_pInt,4_pInt
! --- state update ---
@ -1177,19 +1171,16 @@ subroutine crystallite_integrateStateRKCK45(gg,ii,ee)
use prec, only: pInt, &
pReal
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_e, &
debug_i, &
debug_g, &
debug_selectiveDebugger, &
debug_StateLoopDistribution
use numerics, only: rTol_crystalliteState, &
rTol_crystalliteTemperature, &
subStepSizeCryst, &
stepIncreaseCryst, &
numerics_integrationMode
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP, &
theInc
FEsolving_execIP
use mesh, only: mesh_element, &
mesh_NcpElems, &
mesh_maxNips
@ -1292,7 +1283,7 @@ else
eIter = FEsolving_execElem(1:2)
do e = eIter(1),eIter(2)
iIter(1:2,e) = FEsolving_execIP(1:2,e)
gIter(1:2,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/)
gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
enddo
singleRun = .false.
endif
@ -1345,7 +1336,7 @@ endif
! --- SECOND TO SIXTH RUNGE KUTTA STEP ---
do n = 1,5
do n = 1_pInt,5_pInt
! --- state update ---
@ -1558,7 +1549,7 @@ relTemperatureResiduum = 0.0_pReal
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then
mySizeDotState = constitutive_sizeDotState(g,i,e)
forall (s = 1:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) &
forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) &
relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s)
if (crystallite_Temperature(g,i,e) > 0) &
relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e)
@ -1571,7 +1562,7 @@ relTemperatureResiduum = 0.0_pReal
.and. abs(relTemperatureResiduum(g,i,e)) < rTol_crystalliteTemperature )
#ifndef _OPENMP
if (debug_verbosity > 5 &
if (debug_verbosity > 5_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,i8,1x,i3,1x,i3)') '<< CRYST >> updateState at el ip g ',e,i,g
write(6,*)
@ -1614,7 +1605,8 @@ relTemperatureResiduum = 0.0_pReal
crystallite_todo(g,i,e) = .false. ! ... integration done
if (debug_verbosity > 4) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(6,numerics_integrationMode) = debug_StateLoopDistribution(6,numerics_integrationMode) + 1
debug_StateLoopDistribution(6,numerics_integrationMode) =&
debug_StateLoopDistribution(6,numerics_integrationMode) + 1_pInt
!$OMP END CRITICAL (distributionState)
endif
else
@ -1666,8 +1658,6 @@ use debug, only: debug_verbosity, &
debug_StateLoopDistribution
use numerics, only: rTol_crystalliteState, &
rTol_crystalliteTemperature, &
subStepSizeCryst, &
stepIncreaseCryst, &
numerics_integrationMode
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
@ -1876,9 +1866,9 @@ relTemperatureResiduum = 0.0_pReal
! --- relative residui ---
forall (s = 1:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) &
forall (s = 1_pInt:mySizeDotState, abs(constitutive_state(g,i,e)%p(s)) > 0.0_pReal) &
relStateResiduum(s,g,i,e) = stateResiduum(s,g,i,e) / constitutive_state(g,i,e)%p(s)
if (crystallite_Temperature(g,i,e) > 0) &
if (crystallite_Temperature(g,i,e) > 0_pInt) &
relTemperatureResiduum(g,i,e) = temperatureResiduum(g,i,e) / crystallite_Temperature(g,i,e)
!$OMP FLUSH(relStateResiduum,relTemperatureResiduum)
@ -2136,10 +2126,6 @@ subroutine crystallite_integrateStateFPI(gg,ii,ee)
use prec, only: pInt, &
pReal
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_e, &
debug_i, &
debug_g, &
debug_StateLoopDistribution
use numerics, only: nState, &
numerics_integrationMode
@ -2148,9 +2134,7 @@ use FEsolving, only: FEsolving_execElem, &
use mesh, only: mesh_element, &
mesh_NcpElems
use material, only: homogenization_Ngrains
use constitutive, only: constitutive_sizeDotState, &
constitutive_state, &
constitutive_dotState, &
use constitutive, only: constitutive_dotState, &
constitutive_collectDotState, &
constitutive_dotTemperature, &
constitutive_microstructure, &
@ -2192,7 +2176,7 @@ else
eIter = FEsolving_execElem(1:2)
do e = eIter(1),eIter(2)
iIter(1:2,e) = FEsolving_execIP(1:2,e)
gIter(1:2,e) = (/1,homogenization_Ngrains(mesh_element(3,e))/)
gIter(1:2,e) = [1_pInt,homogenization_Ngrains(mesh_element(3,e))]
enddo
singleRun = .false.
endif
@ -2346,7 +2330,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
if (debug_verbosity > 4) then
!$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(NiterationState,numerics_integrationMode) = &
debug_StateLoopDistribution(NiterationState,numerics_integrationMode) + 1
debug_StateLoopDistribution(NiterationState,numerics_integrationMode) + 1_pInt
!$OMP END CRITICAL (distributionState)
endif
endif
@ -2411,11 +2395,8 @@ endsubroutine
subroutine crystallite_updateState(done, converged, g, i, e)
!*** variables and functions from other modules ***!
use prec, only: pReal, &
pInt, &
pLongInt
use numerics, only: rTol_crystalliteState, &
numerics_integrationMode
use prec, only: pInt
use numerics, only: rTol_crystalliteState
use constitutive, only: constitutive_dotState, &
constitutive_previousDotState, &
constitutive_sizeDotState, &
@ -2424,10 +2405,10 @@ use constitutive, only: constitutive_dotState, &
constitutive_aTolState, &
constitutive_microstructure
use debug, only: debug_verbosity, &
debug_g, &
debug_i, &
debug_selectiveDebugger, &
debug_e, &
debug_selectiveDebugger
debug_i, &
debug_g
!*** input variables ***!
integer(pInt), intent(in):: e, & ! element index
@ -2513,13 +2494,10 @@ endsubroutine
subroutine crystallite_updateTemperature(done, converged, g, i, e)
!*** variables and functions from other modules ***!
use prec, only: pReal, &
pInt, &
pLongInt
use prec, only: pInt
use numerics, only: rTol_crystalliteTemperature
use constitutive, only: constitutive_dotTemperature
use debug, only: debug_verbosity
use debug, only: debug_verbosity
!*** input variables ***!
integer(pInt), intent(in):: e, & ! element index
i, & ! integration point index
@ -2591,17 +2569,16 @@ use numerics, only: nStress, &
relevantStrain, &
numerics_integrationMode
use debug, only: debug_verbosity, &
debug_g, &
debug_i, &
debug_e, &
debug_selectiveDebugger, &
debug_e, &
debug_i, &
debug_g, &
debug_cumLpCalls, &
debug_cumLpTicks, &
debug_StressLoopDistribution, &
debug_LeapfrogBreakDistribution
use constitutive, only: constitutive_homogenizedC, &
constitutive_LpAndItsTangent, &
constitutive_state
use constitutive, only: constitutive_LpAndItsTangent, &
constitutive_homogenizedC
use math, only: math_mul33x33, &
math_mul33xx33, &
math_mul66x6, &
@ -2739,14 +2716,14 @@ steplength_max = 16.0_pReal
jacoCounter = 0_pInt
LpLoop: do
NiterationStress = NiterationStress + 1
NiterationStress = NiterationStress + 1_pInt
!* too many loops required ?
if (NiterationStress > nStress) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
if (debug_verbosity > 4_pInt) then
write(6,'(a,i8,1x,i2,1x,i3)') '<< CRYST >> integrateStress reached loop limit at el ip g ',e,i,g
write(6,*)
endif
@ -2764,7 +2741,7 @@ LpLoop: do
Tstar_v = 0.5_pReal * math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB) - math_I3))
p_hydro = sum(Tstar_v(1:3)) / 3.0_pReal
forall(n=1:3) Tstar_v(n) = Tstar_v(n) - p_hydro ! get deviatoric stress tensor
forall(n=1_pInt:3_pInt) Tstar_v(n) = Tstar_v(n) - p_hydro ! get deviatoric stress tensor
!* calculate plastic velocity gradient and its tangent according to constitutive law
@ -2875,7 +2852,7 @@ LpLoop: do
if (debug_verbosity > 4) then
!$OMP CRITICAL (distributionLeapfrogBreak)
debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) = &
debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) + 1
debug_LeapfrogBreakDistribution(NiterationStress,numerics_integrationMode) + 1_pInt
!$OMP END CRITICAL (distributionLeapfrogBreak)
endif
endif
@ -2887,7 +2864,7 @@ LpLoop: do
if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then
dT_dLp = 0.0_pReal
do h=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3
do h=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt; do m=1_pInt,3_pInt
dT_dLp(3*(h-1)+j,3*(k-1)+l) = dT_dLp(3*(h-1)+j,3*(k-1)+l) + C(h,j,l,m) * AB(k,m) + C(h,j,m,l) * BTA(m,k)
enddo; enddo; enddo; enddo; enddo
dT_dLp = -0.5_pReal * dt * dT_dLp
@ -2914,13 +2891,13 @@ LpLoop: do
return
endif
deltaLp = 0.0_pReal
do k=1,3; do l=1,3; do m=1,3; do n=1,3
deltaLp(k,l) = deltaLp(k,l) - inv_dR_dLp(3*(k-1)+l,3*(m-1)+n) * residuum(m,n)
do k=1_pInt,3_pInt; do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt
deltaLp(k,l) = deltaLp(k,l) - inv_dR_dLp(3_pInt*(k-1_pInt)+l,3_pInt*(m-1_pInt)+n) * residuum(m,n)
enddo; enddo; enddo; enddo
gradientR = 0.0_pReal
do k=1,3; do l=1,3; do m=1,3; do n=1,3
gradientR(k,l) = gradientR(k,l) + dR_dLp(3*(k-1)+l,3*(m-1)+n) * residuum(m,n)
do k=1_pInt,3_pInt; do l=1_pInt,3_pInt; do m=1_pInt,3_pInt; do n=1_pInt,3_pInt
gradientR(k,l) = gradientR(k,l) + dR_dLp(3*(k-1)+l,3_pInt*(m-1_pInt)+n) * residuum(m,n)
enddo; enddo; enddo; enddo
gradientR = gradientR / math_norm33(gradientR)
expectedImprovement = math_mul33xx33(deltaLp, gradientR)
@ -2942,8 +2919,8 @@ call math_invert33(invFp_new,Fp_new,det,error)
if (error) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',e,i,g, &
' ; iteration ', NiterationStress
write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',&
e,i,g, ' ; iteration ', NiterationStress
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new)
@ -2957,7 +2934,7 @@ Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resu
!* add volumetric component to 2nd Piola-Kirchhoff stress and calculate 1st Piola-Kirchhoff stress
forall (n=1:3) Tstar_v(n) = Tstar_v(n) + p_hydro
forall (n=1_pInt:3_pInt) Tstar_v(n) = Tstar_v(n) + p_hydro
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), math_transpose33(invFp_new)))
@ -2976,11 +2953,11 @@ crystallite_integrateStress = .true.
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) &
.and. numerics_integrationMode == 1_pInt) then
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1e6
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1.0e6_pReal
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Cauchy / MPa', &
math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose33(Fg_new)) / 1e6 / math_det33(Fg_new)
math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose33(Fg_new)) / 1.0e6_pReal / math_det33(Fg_new)
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fe Lp Fe^-1', &
math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order
math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
endif
#endif
@ -2988,7 +2965,7 @@ endif
if (debug_verbosity > 4) then
!$OMP CRITICAL (distributionStress)
debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) = &
debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1
debug_StressLoopDistribution(NiterationStress,numerics_integrationMode) + 1_pInt
!$OMP END CRITICAL (distributionStress)
endif
@ -3007,22 +2984,17 @@ use prec, only: pInt, &
use math, only: math_pDecomposition, &
math_RtoQuaternion, &
math_QuaternionDisorientation, &
inDeg, &
math_qConj
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
use IO, only: IO_warning
use material, only: material_phase, &
homogenization_Ngrains, &
phase_constitution, &
phase_localConstitution, &
phase_constitutionInstance
use mesh, only: mesh_element, &
mesh_ipNeighborhood, &
FE_NipNeighbors
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_e, debug_i, debug_g
use constitutive_nonlocal, only: constitutive_nonlocal_structure, &
constitutive_nonlocal_updateCompatibility
@ -3054,7 +3026,7 @@ logical error
!$OMP PARALLEL DO PRIVATE(error,U,R,orientation)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e))
do g = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
call math_pDecomposition(crystallite_Fe(1:3,1:3,g,i,e), U, R, error) ! polar decomposition of Fe
if (error) then
@ -3088,7 +3060,7 @@ logical error
! --- calculate disorientation between me and my neighbor ---
do n = 1,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors
do n = 1_pInt,FE_NipNeighbors(mesh_element(2,e)) ! loop through my neighbors
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
if ((neighboring_e > 0) .and. (neighboring_i > 0)) then ! if neighbor exists
@ -3181,7 +3153,7 @@ function crystallite_postResults(&
crystallite_postResults(c+1) = real(crystallite_sizePostResults(crystID),pReal) ! size of results from cryst
c = c + 1_pInt
do o = 1,crystallite_Noutput(crystID)
do o = 1_pInt,crystallite_Noutput(crystID)
mySize = 0_pInt
select case(crystallite_output(o,crystID))
case ('phase')

17
code/homogenization.f90
View File

@ -72,6 +72,7 @@ CONTAINS
!* Module initialization *
!**************************************
subroutine homogenization_init(Temperature)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: pReal,pInt
use math, only: math_I3
use debug, only: debug_verbosity
@ -288,10 +289,10 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_converged, &
crystallite_stressAndItsTangent, &
crystallite_orientations
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
use debug, only: debug_verbosity, &
debug_e, &
debug_i, &
debug_selectiveDebugger, &
debug_MaterialpointLoopDistribution, &
debug_MaterialpointStateLoopDistribution
use math, only: math_pDecomposition
@ -358,7 +359,7 @@ subroutine materialpoint_stressAndItsTangent(&
if ( materialpoint_converged(i,e) ) then
#ifndef _OPENMP
if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,1x,f10.8,1x,a,1x,f10.8,1x,a,/)') '<< HOMOG >> winding forward from', &
write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,/)') '<< HOMOG >> winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent'
endif
@ -409,7 +410,7 @@ subroutine materialpoint_stressAndItsTangent(&
#ifndef _OPENMP
if (debug_verbosity > 2 .and. ((e == debug_e .and. i == debug_i) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,1x,f10.8,/)') &
write(6,'(a,1x,f12.8,/)') &
'<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',&
materialpoint_subStep(i,e)
endif
@ -593,7 +594,7 @@ subroutine homogenization_partitionDeformation(&
el & ! element
)
use prec, only: pReal,pInt
use prec, only: pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_partionedF0,crystallite_partionedF
@ -634,7 +635,7 @@ function homogenization_updateState(&
ip, & ! integration point
el & ! element
)
use prec, only: pReal,pInt
use prec, only: pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_P,crystallite_dPdF,crystallite_partionedF,crystallite_partionedF0 ! modified <<<updated 31.07.2009>>>
@ -682,7 +683,7 @@ subroutine homogenization_averageStressAndItsTangent(&
ip, & ! integration point
el & ! element
)
use prec, only: pReal,pInt
use prec, only: pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_P,crystallite_dPdF
@ -724,7 +725,7 @@ subroutine homogenization_averageTemperature(&
ip, & ! integration point
el & ! element
)
use prec, only: pReal,pInt
use prec, only: pInt
use mesh, only: mesh_element
use material, only: homogenization_type, homogenization_maxNgrains
use crystallite, only: crystallite_Temperature

364
code/homogenization_RGC.f90
View File

@ -63,18 +63,19 @@ CONTAINS
!* Module initialization *
!**************************************
subroutine homogenization_RGC_init(&
file & ! file pointer to material configuration
myFile & ! file pointer to material configuration
)
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 debug, only: debug_verbosity
use math, only: math_Mandel3333to66, math_Voigt66to3333,math_I3,math_sampleRandomOri,math_EulerToR,inRad
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use IO
use material
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 4
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 4_pInt
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,e, output, mySize, myInstance
character(len=64) tag
character(len=1024) line
@ -86,7 +87,7 @@ subroutine homogenization_RGC_init(&
#include "compilation_info.f90"
!$OMP END CRITICAL (write2out)
maxNinstance = count(homogenization_type == homogenization_RGC_label)
maxNinstance = int(count(homogenization_type == homogenization_RGC_label),pInt)
if (maxNinstance == 0) return
allocate(homogenization_RGC_sizeState(maxNinstance)); homogenization_RGC_sizeState = 0_pInt
@ -100,20 +101,20 @@ subroutine homogenization_RGC_init(&
allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),maxNinstance))
homogenization_RGC_sizePostResult = 0_pInt
allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems))
forall (i = 1:mesh_maxNips,e = 1:mesh_NcpElems)
forall (i = 1_pInt:mesh_maxNips,e = 1_pInt:mesh_NcpElems)
homogenization_RGC_orientation(:,:,i,e) = math_I3
end forall
rewind(file)
rewind(myFile)
line = ''
section = 0
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to <homogenization>
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
@ -212,9 +213,9 @@ subroutine homogenization_RGC_init(&
homogenization_RGC_sizeState(i) &
= 3*(homogenization_RGC_Ngrains(1,i)-1_pInt)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) &
+ 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1_pInt)*homogenization_RGC_Ngrains(3,i) &
+ 3*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1_pInt) &
= 3_pInt*(homogenization_RGC_Ngrains(1,i)-1_pInt)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) &
+ 3_pInt*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1_pInt)*homogenization_RGC_Ngrains(3,i) &
+ 3_pInt*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1_pInt) &
+ 8_pInt ! (1) Average constitutive work, (2-4) Overall mismatch, (5) Average penalty energy,
! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component
enddo
@ -254,9 +255,9 @@ subroutine homogenization_RGC_partitionDeformation(&
)
use prec, only: pReal,pInt,p_vec
use debug, only: debug_verbosity
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance
use FEsolving, only: theInc,cycleCounter,theTime
use FEsolving, only: theInc,cycleCounter
implicit none
@ -291,7 +292,7 @@ subroutine homogenization_RGC_partitionDeformation(&
!* Compute the deformation gradient of individual grains due to relaxations
homID = homogenization_typeInstance(mesh_element(3,el))
F = 0.0_pReal
do iGrain = 1,homogenization_Ngrains(mesh_element(3,el))
do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el))
iGrain3 = homogenization_RGC_grain1to3(iGrain,homID)
do iFace = 1_pInt,nFace
intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain
@ -339,12 +340,11 @@ function homogenization_RGC_updateState(&
use prec, only: pReal,pInt,p_vec
use debug, only: debug_verbosity, debug_e, debug_i
use math, only: math_invert
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains,homogenization_typeInstance, &
homogenization_Ngrains
use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC, &
maxdRelax_RGC,viscPower_RGC,viscModus_RGC,refRelaxRate_RGC
use FEsolving, only: theInc,cycleCounter,theTime
implicit none
@ -379,8 +379,8 @@ function homogenization_RGC_updateState(&
homID = homogenization_typeInstance(mesh_element(3,el))
nGDim = homogenization_RGC_Ngrains(:,homID)
nGrain = homogenization_Ngrains(mesh_element(3,el))
nIntFaceTot = (nGDim(1)-1)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1)*nGDim(3) &
+ nGDim(1)*nGDim(2)*(nGDim(3)-1)
nIntFaceTot = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) &
+ nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt)
!* Allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
allocate(resid(3_pInt*nIntFaceTot)); resid = 0.0_pReal
@ -390,10 +390,10 @@ function homogenization_RGC_updateState(&
drelax = state%p(1:3_pInt*nIntFaceTot) - state0%p(1:3_pInt*nIntFaceTot)
!* Debugging the obtained state
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Obtained state: '
do i = 1,3*nIntFaceTot
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,2(e15.8,1x))')state%p(i)
enddo
write(6,*)' '
@ -407,16 +407,16 @@ function homogenization_RGC_updateState(&
call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el,homID)
!* Debugging the mismatch, stress and penalties of grains
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
do iGrain = 1,nGrain
do iGrain = 1_pInt,nGrain
write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e15.8))')'Mismatch magnitude of grain(',iGrain,') :',NN(1,iGrain),NN(2,iGrain),NN(3,iGrain)
write(6,*)' '
write(6,'(1x,a30,1x,i3)')'Stress and penalties of grain: ',iGrain
do i = 1,3
write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1,3), &
(R(i,j,iGrain), j = 1,3), &
(D(i,j,iGrain), j = 1,3)
do i = 1_pInt,3_pInt
write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1_pInt,3_pInt), &
(R(i,j,iGrain), j = 1_pInt,3_pInt), &
(D(i,j,iGrain), j = 1_pInt,3_pInt)
enddo
write(6,*)' '
enddo
@ -426,7 +426,7 @@ function homogenization_RGC_updateState(&
!* -------------------------------------------------------------------------------------------------------------
!*** Computing the residual stress from the balance of traction at all (interior) interfaces
do iNum = 1,nIntFaceTot
do iNum = 1_pInt,nIntFaceTot
faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index)
!* Identify the left/bottom/back grain (-|N)
@ -443,23 +443,23 @@ function homogenization_RGC_updateState(&
normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the interface normal
!* Compute the residual of traction at the interface (in local system, 4-dimensional index)
do i = 1,3
do i = 1_pInt,3_pInt
tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1_pInt)))/(refRelaxRate_RGC*dt))**viscPower_RGC, &
drelax(i+3*(iNum-1))) ! contribution from the relaxation viscosity
do j = 1,3
drelax(i+3*(iNum-1_pInt))) ! contribution from the relaxation viscosity
do j = 1_pInt,3_pInt
tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) &
+ (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j)
! contribution from material stress P, mismatch penalty R, and volume penalty D
! projected into the interface
resid(i+3*(iNum-1)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array
resid(i+3_pInt*(iNum-1_pInt)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array
enddo
enddo
!* Debugging the residual stress
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum
write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1,3)
write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1_pInt,3_pInt)
write(6,*)' '
!$OMP END CRITICAL (write2out)
endif
@ -468,13 +468,13 @@ function homogenization_RGC_updateState(&
!* -------------------------------------------------------------------------------------------------------------
!*** Convergence check for stress residual
stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress
stresLoc = maxloc(abs(P)) ! get the location of the maximum stress
residMax = maxval(abs(tract)) ! get the maximum of the residual
residLoc = maxloc(abs(tract)) ! get the position of the maximum residual
stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress
stresLoc = int(maxloc(abs(P)),pInt) ! get the location of the maximum stress
residMax = maxval(abs(tract)) ! get the maximum of the residual
residLoc = int(maxloc(abs(tract)),pInt) ! get the position of the maximum residual
!* Debugging the convergent criteria
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
if (debug_verbosity == 4_pInt .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a)')' '
write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el
@ -504,34 +504,34 @@ function homogenization_RGC_updateState(&
!* ... all energy densities computed by time-integration
constitutiveWork = state%p(3*nIntFaceTot+1)
penaltyEnergy = state%p(3*nIntFaceTot+5)
do iGrain = 1,homogenization_Ngrains(mesh_element(3,el)) ! time-integration loop for the calculating the work and energy
do i = 1,3
do j = 1,3
constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/dble(nGrain)
penaltyEnergy = penaltyEnergy + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/dble(nGrain)
do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) ! time-integration loop for the calculating the work and energy
do i = 1_pInt,3_pInt
do j = 1_pInt,3_pInt
constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
penaltyEnergy = penaltyEnergy + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
enddo
enddo
enddo
state%p(3*nIntFaceTot+1) = constitutiveWork ! the bulk mechanical/constitutive work
state%p(3*nIntFaceTot+2) = sum(NN(1,:))/dble(nGrain) ! the overall mismatch of all interface normal to e1-direction
state%p(3*nIntFaceTot+3) = sum(NN(2,:))/dble(nGrain) ! the overall mismatch of all interface normal to e2-direction
state%p(3*nIntFaceTot+4) = sum(NN(3,:))/dble(nGrain) ! the overall mismatch of all interface normal to e3-direction
state%p(3*nIntFaceTot+2) = sum(NN(1,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e1-direction
state%p(3*nIntFaceTot+3) = sum(NN(2,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e2-direction
state%p(3*nIntFaceTot+4) = sum(NN(3,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e3-direction
state%p(3*nIntFaceTot+5) = penaltyEnergy ! the overall penalty energy
state%p(3*nIntFaceTot+6) = volDiscrep ! the overall volume discrepancy
state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/dble(3*nIntFaceTot) ! the average rate of relaxation vectors
state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) ! the average rate of relaxation vectors
state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors
if (debug_verbosity == 4 .and. debug_e == el .and. debug_i == ip) then
if (debug_verbosity == 4_pInt .and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30,1x,e15.8)')'Constitutive work: ',constitutiveWork
write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/dble(nGrain), &
sum(NN(2,:))/dble(nGrain), &
sum(NN(3,:))/dble(nGrain)
write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/real(nGrain,pReal), &
sum(NN(2,:))/real(nGrain,pReal), &
sum(NN(3,:))/real(nGrain,pReal)
write(6,'(1x,a30,1x,e15.8)')'Penalty energy: ',penaltyEnergy
write(6,'(1x,a30,1x,e15.8)')'Volume discrepancy: ',volDiscrep
write(6,*)''
write(6,'(1x,a30,1x,e15.8)')'Maximum relaxation rate: ',maxval(abs(drelax))/dt
write(6,'(1x,a30,1x,e15.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/dble(3*nIntFaceTot)
write(6,'(1x,a30,1x,e15.8)')'Average relaxation rate: ',sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal)
write(6,*)''
call flush(6)
!$OMP END CRITICAL (write2out)
@ -575,20 +575,20 @@ function homogenization_RGC_updateState(&
!* ... of the constitutive stress tangent,
!* assembled from dPdF or material constitutive model "smatrix"
allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot)); smatrix = 0.0_pReal
do iNum = 1,nIntFaceTot
do iNum = 1_pInt,nIntFaceTot
faceID = homogenization_RGC_interface1to4(iNum,homID) ! assembling of local dPdF into global Jacobian matrix
!* Identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem
iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate into global grain ID
intFaceN = homogenization_RGC_getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system
intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system
normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the interface normal
do iFace = 1,nFace
do iFace = 1_pInt,nFace
intFaceN = homogenization_RGC_getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface
mornN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get normal of the interfaces
iMun = homogenization_RGC_interface4to1(intFaceN,homID) ! translate the interfaces ID into local 4-dimensional index
if (iMun .gt. 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3
do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l)
enddo;enddo;enddo;enddo
! projecting the material tangent dPdF into the interface
@ -598,16 +598,16 @@ function homogenization_RGC_updateState(&
!* Identify the right/up/front grain (+|P)
iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem
iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate sytem
iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate into global grain ID
intFaceP = homogenization_RGC_getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system
intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the connecting interface in local coordinate system
normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the interface normal
do iFace = 1,nFace
do iFace = 1_pInt,nFace
intFaceP = homogenization_RGC_getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface
mornP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get normal of the interfaces
iMun = homogenization_RGC_interface4to1(intFaceP,homID) ! translate the interfaces ID into local 4-dimensional index
if (iMun .gt. 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3
do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l)
enddo;enddo;enddo;enddo
endif
@ -615,11 +615,11 @@ function homogenization_RGC_updateState(&
enddo
!* Debugging the global Jacobian matrix of stress tangent
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of stress'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1,3*nIntFaceTot)
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -631,7 +631,7 @@ function homogenization_RGC_updateState(&
allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot)); pmatrix = 0.0_pReal
allocate(p_relax(3*nIntFaceTot)); p_relax = 0.0_pReal
allocate(p_resid(3*nIntFaceTot)); p_resid = 0.0_pReal
do ipert = 1,3*nIntFaceTot
do ipert = 1_pInt,3_pInt*nIntFaceTot
p_relax = relax
p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector
state%p(1:3*nIntFaceTot) = p_relax
@ -641,25 +641,25 @@ function homogenization_RGC_updateState(&
!* Computing the global stress residual array from the perturbed state
p_resid = 0.0_pReal
do iNum = 1,nIntFaceTot
do iNum = 1_pInt,nIntFaceTot
faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index)
!* Identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = homogenization_RGC_getInterface(2*faceID(1),iGr3N) ! identifying the interface ID of the grain
intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N) ! identifying the interface ID of the grain
normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the corresponding interface normal
!* Identify the right/up/front grain (+|P)
iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index)
iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = homogenization_RGC_getInterface(2*faceID(1)-1,iGr3P) ! identifying the interface ID of the grain
intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the interface ID of the grain
normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the corresponding normal
!* Compute the residual stress (contribution of mismatch and volume penalties) from perturbed state at all interfaces
do i = 1,3
do j = 1,3
do i = 1_pInt,3_pInt
do j = 1_pInt,3_pInt
p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) &
+ (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) &
+ (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) &
@ -674,8 +674,8 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1,3*nIntFaceTot)
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -684,18 +684,18 @@ function homogenization_RGC_updateState(&
!* ... of the numerical viscosity traction "rmatrix"
allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot)); rmatrix = 0.0_pReal
forall (i=1:3*nIntFaceTot) &
forall (i=1_pInt:3_pInt*nIntFaceTot) &
rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* &
(abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal)
! tangent due to numerical viscosity traction appears
! only in the main diagonal term
!* Debugging the global Jacobian matrix of numerical viscosity tangent
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1,3*nIntFaceTot)
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -708,8 +708,8 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix (total)'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1,3*nIntFaceTot)
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -727,7 +727,7 @@ function homogenization_RGC_updateState(&
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian inverse'
do i = 1,3*nIntFaceTot
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
@ -737,8 +737,8 @@ function homogenization_RGC_updateState(&
!* Calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
drelax = 0.0_pReal
do i = 1,3*nIntFaceTot
do j = 1,3*nIntFaceTot
do i = 1_pInt,3_pInt*nIntFaceTot
do j = 1_pInt,3_pInt*nIntFaceTot
drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable
enddo
enddo
@ -754,10 +754,10 @@ function homogenization_RGC_updateState(&
endif
!* Debugging the return state
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Returned state: '
do i = 1,3*nIntFaceTot
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,2(e15.8,1x))')state%p(i)
enddo
write(6,*)' '
@ -785,7 +785,7 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
use prec, only: pReal,pInt,p_vec
use debug, only: debug_verbosity
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains,homogenization_Ngrains,homogenization_typeInstance
use math, only: math_Plain3333to99
implicit none
@ -804,13 +804,13 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
Ngrains = homogenization_Ngrains(mesh_element(3,el))
!* Debugging the grain tangent
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
do iGrain = 1,Ngrains
do iGrain = 1_pInt,Ngrains
dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain))
write(6,'(1x,a30,1x,i3)')'Stress tangent of grain: ',iGrain
do i = 1,9
write(6,'(1x,(e15.8,1x))') (dPdF99(i,j), j = 1,9)
do i = 1_pInt,9_pInt
write(6,'(1x,(e15.8,1x))') (dPdF99(i,j), j = 1_pInt,9_pInt)
enddo
write(6,*)' '
enddo
@ -819,8 +819,8 @@ subroutine homogenization_RGC_averageStressAndItsTangent(&
endif
!* Computing the average first Piola-Kirchhoff stress P and the average tangent dPdF
avgP = sum(P,3)/dble(Ngrains)
dAvgPdAvgF = sum(dPdF,5)/dble(Ngrains)
avgP = sum(P,3)/real(Ngrains,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal)
endsubroutine
@ -834,7 +834,7 @@ function homogenization_RGC_averageTemperature(&
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains, homogenization_Ngrains
implicit none
@ -842,11 +842,11 @@ function homogenization_RGC_averageTemperature(&
real(pReal), dimension (homogenization_maxNgrains), intent(in) :: Temperature
integer(pInt), intent(in) :: ip,el
real(pReal) homogenization_RGC_averageTemperature
integer(pInt) homID, Ngrains
integer(pInt) :: Ngrains
!* Computing the average temperature
Ngrains = homogenization_Ngrains(mesh_element(3,el))
homogenization_RGC_averageTemperature = sum(Temperature(1:Ngrains))/dble(Ngrains)
homogenization_RGC_averageTemperature = sum(Temperature(1:Ngrains))/real(Ngrains,pReal)
endfunction
@ -861,7 +861,7 @@ pure function homogenization_RGC_postResults(&
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use material, only: homogenization_typeInstance,homogenization_Noutput,homogenization_Ngrains
use material, only: homogenization_typeInstance,homogenization_Noutput
implicit none
!* Definition of variables
@ -873,34 +873,34 @@ pure function homogenization_RGC_postResults(&
homogenization_RGC_postResults
homID = homogenization_typeInstance(mesh_element(3,el))
nIntFaceTot = (homogenization_RGC_Ngrains(1,homID)-1)*homogenization_RGC_Ngrains(2,homID)*homogenization_RGC_Ngrains(3,homID) + &
homogenization_RGC_Ngrains(1,homID)*(homogenization_RGC_Ngrains(2,homID)-1)*homogenization_RGC_Ngrains(3,homID) + &
homogenization_RGC_Ngrains(1,homID)*homogenization_RGC_Ngrains(2,homID)*(homogenization_RGC_Ngrains(3,homID)-1)
nIntFaceTot=(homogenization_RGC_Ngrains(1,homID)-1_pInt)*homogenization_RGC_Ngrains(2,homID)*homogenization_RGC_Ngrains(3,homID)&
+ homogenization_RGC_Ngrains(1,homID)*(homogenization_RGC_Ngrains(2,homID)-1_pInt)*homogenization_RGC_Ngrains(3,homID)&
+ homogenization_RGC_Ngrains(1,homID)*homogenization_RGC_Ngrains(2,homID)*(homogenization_RGC_Ngrains(3,homID)-1_pInt)
c = 0_pInt
homogenization_RGC_postResults = 0.0_pReal
do o = 1,homogenization_Noutput(mesh_element(3,el))
do o = 1_pInt,homogenization_Noutput(mesh_element(3,el))
select case(homogenization_RGC_output(o,homID))
case('constitutivework')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1)
c = c + 1
c = c + 1_pInt
case('magnitudemismatch')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2)
homogenization_RGC_postResults(c+2) = state%p(3*nIntFaceTot+3)
homogenization_RGC_postResults(c+3) = state%p(3*nIntFaceTot+4)
c = c + 3
c = c + 3_pInt
case('penaltyenergy')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+5)
c = c + 1
c = c + 1_pInt
case('volumediscrepancy')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+6)
c = c + 1
c = c + 1_pInt
case('averagerelaxrate')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+7)
c = c + 1
c = c + 1_pInt
case('maximumrelaxrate')
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+8)
c = c + 1
c = c + 1_pInt
end select
enddo
@ -960,24 +960,24 @@ subroutine homogenization_RGC_stressPenalty(&
!* -------------------------------------------------------------------------------------------------------------
!*** Computing the mismatch and penalty stress tensor of all grains
do iGrain = 1,homogenization_Ngrains(mesh_element(3,el))
do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el))
Gmoduli = homogenization_RGC_equivalentModuli(iGrain,ip,el)
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) ! get the grain ID in local 3-dimensional index (x,y,z)-position
!* Looping over all six interfaces of each grain
do iFace = 1,nFace
do iFace = 1_pInt,nFace
intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the interface normal
iGNghb3 = iGrain3 ! identify the neighboring grain across the interface
iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(dble(intFace(1))/dble(abs(intFace(1))))
iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal),pInt)
if (iGNghb3(1) < 1) iGNghb3(1) = nGDim(1) ! with periodicity along e1 direction
if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1
if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1_pInt
if (iGNghb3(2) < 1) iGNghb3(2) = nGDim(2) ! with periodicity along e2 direction
if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1
if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1_pInt
if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) ! with periodicity along e3 direction
if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1
if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1_pInt
iGNghb = homogenization_RGC_grain3to1(iGNghb3,homID) ! get the ID of the neighboring grain
Gmoduli = homogenization_RGC_equivalentModuli(iGNghb,ip,el) ! collecting the shear modulus and Burgers vector of the neighbor
muGNghb = Gmoduli(1)
@ -987,10 +987,10 @@ subroutine homogenization_RGC_stressPenalty(&
!* Compute the mismatch tensor of all interfaces
nDefNorm = 0.0_pReal
nDef = 0.0_pReal
do i = 1,3
do j = 1,3
do k = 1,3
do l = 1,3
do i = 1_pInt,3_pInt
do j = 1_pInt,3_pInt
do k = 1_pInt,3_pInt
do l = 1_pInt,3_pInt
nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l)! compute the interface mismatch tensor from the jump of deformation gradient
enddo
enddo
@ -1011,10 +1011,10 @@ subroutine homogenization_RGC_stressPenalty(&
! endif
!* Compute the stress penalty of all interfaces
do i = 1,3
do j = 1,3
do k = 1,3
do l = 1,3
do i = 1_pInt,3_pInt
do j = 1_pInt,3_pInt
do k = 1_pInt,3_pInt
do l = 1_pInt,3_pInt
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*homogenization_RGC_xiAlpha(homID) &
*surfCorr(abs(intFace(1)))/homogenization_RGC_dAlpha(abs(intFace(1)),homID) &
*cosh(homogenization_RGC_ciAlpha(homID)*nDefNorm) &
@ -1073,16 +1073,16 @@ subroutine homogenization_RGC_volumePenalty(&
!* Compute the volumes of grains and of cluster
vDiscrep = math_det33(fAvg) ! compute the volume of the cluster
do iGrain = 1,nGrain
do iGrain = 1_pInt,nGrain
gVol(iGrain) = math_det33(fDef(:,:,iGrain)) ! compute the volume of individual grains
vDiscrep = vDiscrep - gVol(iGrain)/dble(nGrain) ! calculate the difference/dicrepancy between
vDiscrep = vDiscrep - gVol(iGrain)/real(nGrain,pReal) ! calculate the difference/dicrepancy between
! the volume of the cluster and the the total volume of grains
enddo
!* Calculate the stress and penalty due to volume discrepancy
vPen = 0.0_pReal
do iGrain = 1,nGrain
vPen(:,:,iGrain) = -1.0_pReal/dble(nGrain)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
do iGrain = 1_pInt,nGrain
vPen(:,:,iGrain) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* &
gVol(iGrain)*transpose(math_inv33(fDef(:,:,iGrain)))
@ -1128,11 +1128,11 @@ function homogenization_RGC_surfaceCorrection(&
invC = 0.0_pReal
call math_invert33(avgC,invC,detF,error)
homogenization_RGC_surfaceCorrection = 0.0_pReal
do iBase = 1,3
do iBase = 1_pInt,3_pInt
intFace = (/iBase,1_pInt,1_pInt,1_pInt/)
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of the interface
do i = 1,3
do j = 1,3
do i = 1_pInt,3_pInt
do j = 1_pInt,3_pInt
homogenization_RGC_surfaceCorrection(iBase) = & ! compute the component of (the inverse of) the stretch in the direction of the normal
homogenization_RGC_surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j)
enddo
@ -1154,8 +1154,6 @@ function homogenization_RGC_equivalentModuli(&
use prec, only: pReal,pInt,p_vec
use constitutive, only: constitutive_homogenizedC,constitutive_averageBurgers
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use material, only: homogenization_typeInstance
implicit none
!* Definition of variables
@ -1188,8 +1186,6 @@ function homogenization_RGC_relaxationVector(&
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use material, only: homogenization_typeInstance
implicit none
!* Definition of variables
@ -1219,7 +1215,6 @@ function homogenization_RGC_interfaceNormal(&
use prec, only: pReal,pInt,p_vec
use math, only: math_mul33x3
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
implicit none
!* Definition of variables
@ -1263,7 +1258,7 @@ function homogenization_RGC_getInterface(&
integer(pInt) iDir
!* Direction of interface normal
iDir = (int(dble(iFace-1)/2.0_pReal)+1)*(-1_pInt)**iFace
iDir = (int(real(iFace-1_pInt,pReal)/2.0_pReal,pInt)+1_pInt)*(-1_pInt)**iFace
homogenization_RGC_getInterface(1) = iDir
!* Identify the interface position by the direction of its normal
@ -1281,8 +1276,7 @@ function homogenization_RGC_grain1to3(&
homID & ! homogenization ID
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use prec, only: pInt,p_vec
implicit none
!* Definition of variables
@ -1306,8 +1300,7 @@ function homogenization_RGC_grain3to1(&
homID & ! homogenization ID
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use prec, only: pInt,p_vec
implicit none
!* Definition of variables
@ -1318,7 +1311,7 @@ function homogenization_RGC_grain3to1(&
!* Get the grain ID
nGDim = homogenization_RGC_Ngrains(:,homID)
homogenization_RGC_grain3to1 = grain3(1) + nGDim(1)*(grain3(2)-1) + nGDim(1)*nGDim(2)*(grain3(3)-1)
homogenization_RGC_grain3to1 = grain3(1) + nGDim(1)*(grain3(2)-1_pInt) + nGDim(1)*nGDim(2)*(grain3(3)-1_pInt)
endfunction
@ -1330,8 +1323,7 @@ function homogenization_RGC_interface4to1(&
homID & ! homogenization ID
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use prec, only: pInt,p_vec
implicit none
!* Definition of variables
@ -1342,22 +1334,22 @@ function homogenization_RGC_interface4to1(&
nGDim = homogenization_RGC_Ngrains(:,homID)
!* Compute the total number of interfaces, which ...
nIntFace(1) = (nGDim(1)-1)*nGDim(2)*nGDim(3) ! ... normal //e1
nIntFace(2) = nGDim(1)*(nGDim(2)-1)*nGDim(3) ! ... normal //e2
nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1) ! ... normal //e3
nIntFace(1) = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! ... normal //e1
nIntFace(2) = nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! ... normal //e2
nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) ! ... normal //e3
!* Get the corresponding interface ID in 1D global array
if (abs(iFace4D(1)) == 1_pInt) then ! ... interface with normal //e1
homogenization_RGC_interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1) &
+ nGDim(2)*nGDim(3)*(iFace4D(2)-1)
homogenization_RGC_interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1_pInt) &
+ nGDim(2)*nGDim(3)*(iFace4D(2)-1_pInt)
if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) homogenization_RGC_interface4to1 = 0_pInt
elseif (abs(iFace4D(1)) == 2_pInt) then ! ... interface with normal //e2
homogenization_RGC_interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1) &
+ nGDim(3)*nGDim(1)*(iFace4D(3)-1) + nIntFace(1)
homogenization_RGC_interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1_pInt) &
+ nGDim(3)*nGDim(1)*(iFace4D(3)-1_pInt) + nIntFace(1)
if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) homogenization_RGC_interface4to1 = 0_pInt
elseif (abs(iFace4D(1)) == 3_pInt) then ! ... interface with normal //e3
homogenization_RGC_interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1) &
+ nGDim(1)*nGDim(2)*(iFace4D(4)-1) + nIntFace(1) + nIntFace(2)
homogenization_RGC_interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1_pInt) &
+ nGDim(1)*nGDim(2)*(iFace4D(4)-1_pInt) + nIntFace(1) + nIntFace(2)
if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) homogenization_RGC_interface4to1 = 0_pInt
endif
@ -1372,7 +1364,6 @@ function homogenization_RGC_interface1to4(&
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
implicit none
!* Definition of variables
@ -1383,30 +1374,53 @@ function homogenization_RGC_interface1to4(&
nGDim = homogenization_RGC_Ngrains(:,homID)
!* Compute the total number of interfaces, which ...
nIntFace(1) = (nGDim(1)-1)*nGDim(2)*nGDim(3) ! ... normal //e1
nIntFace(2) = nGDim(1)*(nGDim(2)-1)*nGDim(3) ! ... normal //e2
nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1) ! ... normal //e3
nIntFace(1) = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! ... normal //e1
nIntFace(2) = nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! ... normal //e2
nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) ! ... normal //e3
!* Get the corresponding interface ID in 4D (normal and local position)
if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! ... interface with normal //e1
homogenization_RGC_interface1to4(1) = 1
if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! ... interface with normal //e1
homogenization_RGC_interface1to4(1) = 1_pInt
homogenization_RGC_interface1to4(3) = mod((iFace1D-1_pInt),nGDim(2))+1_pInt
homogenization_RGC_interface1to4(4) = mod(int(real(iFace1D-1_pInt,pReal)/real(nGDim(2),pReal),pInt),nGDim(3))+1
homogenization_RGC_interface1to4(2) = int(real(iFace1D-1,pReal)/real(nGDim(2),pReal)/real(nGDim(3),pReal),pInt)+1
elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then
! ... interface with normal //e2
homogenization_RGC_interface1to4(1) = 2
homogenization_RGC_interface1to4(4) = mod((iFace1D-nIntFace(1)-1),nGDim(3))+1
homogenization_RGC_interface1to4(2) = mod(int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)),nGDim(1))+1
homogenization_RGC_interface1to4(3) = int(real(iFace1D-nIntFace(1)-1,pReal)/real(nGDim(3),pReal)/real(nGDim(1),pReal))+1
elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then
! ... interface with normal //e3
homogenization_RGC_interface1to4(1) = 3
homogenization_RGC_interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1),nGDim(1))+1
homogenization_RGC_interface1to4(3) = mod(int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/&
real(nGDim(1),pReal),pInt),nGDim(2))+1
homogenization_RGC_interface1to4(4) = int(real(iFace1D-nIntFace(2)-nIntFace(1)-1,pReal)/&
real(nGDim(1),pReal)/real(nGDim(2),pReal),pInt)+1
homogenization_RGC_interface1to4(4) = mod(&
int(&
real(iFace1D-1_pInt,pReal)/&
real(nGDim(2),pReal)&
,pInt)&
,nGDim(3))+1_pInt
homogenization_RGC_interface1to4(2) = int(&
real(iFace1D-1_pInt,pReal)/&
real(nGDim(2),pReal)/&
real(nGDim(3),pReal)&
,pInt)+1_pInt
elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! ... interface with normal //e2
homogenization_RGC_interface1to4(1) = 2_pInt
homogenization_RGC_interface1to4(4) = mod((iFace1D-nIntFace(1)-1_pInt),nGDim(3))+1_pInt
homogenization_RGC_interface1to4(2) = mod(&
int(&
real(iFace1D-nIntFace(1)-1_pInt,pReal)/&
real(nGDim(3),pReal)&
,pInt)&
,nGDim(1))+1_pInt
homogenization_RGC_interface1to4(3) = int(&
real(iFace1D-nIntFace(1)-1_pInt,pReal)/&
real(nGDim(3),pReal)/&
real(nGDim(1),pReal)&
,pInt)+1_pInt
elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! ... interface with normal //e3
homogenization_RGC_interface1to4(1) = 3_pInt
homogenization_RGC_interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1_pInt),nGDim(1))+1_pInt
homogenization_RGC_interface1to4(3) = mod(&
int(&
real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/&
real(nGDim(1),pReal)&
,pInt)&
,nGDim(2))+1_pInt
homogenization_RGC_interface1to4(4) = int(&
real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/&
real(nGDim(1),pReal)/&
real(nGDim(2),pReal)&
,pInt)+1_pInt
endif
endfunction
@ -1442,18 +1456,18 @@ subroutine homogenization_RGC_grainDeformation(&
integer(pInt), dimension (3) :: iGrain3
integer(pInt) homID, iGrain,iFace,i,j
!
integer(pInt), parameter :: nFace = 6
integer(pInt), parameter :: nFace = 6_pInt
!* Compute the deformation gradient of individual grains due to relaxations
homID = homogenization_typeInstance(mesh_element(3,el))
F = 0.0_pReal
do iGrain = 1,homogenization_Ngrains(mesh_element(3,el))
do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el))
iGrain3 = homogenization_RGC_grain1to3(iGrain,homID)
do iFace = 1,nFace
do iFace = 1_pInt,nFace
intFace = homogenization_RGC_getInterface(iFace,iGrain3)
aVect = homogenization_RGC_relaxationVector(intFace,state,homID)
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el)
forall (i=1:3,j=1:3) &
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations
enddo
F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! relaxed deformation gradient

54
code/homogenization_isostrain.f90
View File

@ -59,17 +59,18 @@ CONTAINS
!* Module initialization *
!**************************************
subroutine homogenization_isostrain_init(&
file & ! file pointer to material configuration
myFile & ! file pointer to material configuration
)
use, intrinsic :: iso_fortran_env
use prec, only: pInt, pReal
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 math, only: math_Mandel3333to66, math_Voigt66to3333
use IO
use material
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j, output, mySize
integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: maxNchunks = 2_pInt
integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions
integer(pInt) section, i, j, output, mySize
integer :: maxNinstance, k !no pInt (stores a system dependen value from 'count'
character(len=64) tag
character(len=1024) line
@ -81,7 +82,7 @@ subroutine homogenization_isostrain_init(&
!$OMP END CRITICAL (write2out)
maxNinstance = count(homogenization_type == homogenization_isostrain_label)
if (maxNinstance == 0_pInt) return
if (maxNinstance == 0) return
allocate(homogenization_isostrain_sizeState(maxNinstance)) ; homogenization_isostrain_sizeState = 0_pInt
allocate(homogenization_isostrain_sizePostResults(maxNinstance)); homogenization_isostrain_sizePostResults = 0_pInt
@ -91,16 +92,16 @@ subroutine homogenization_isostrain_init(&
allocate(homogenization_isostrain_output(maxval(homogenization_Noutput), &
maxNinstance)) ; homogenization_isostrain_output = ''
rewind(file)
rewind(myFile)
line = ''
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= material_partHomogenization) ! wind forward to <homogenization>
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
@ -121,18 +122,18 @@ subroutine homogenization_isostrain_init(&
endif
enddo
100 do i = 1_pInt,maxNinstance ! sanity checks
100 do k = 1,maxNinstance ! sanity checks
enddo
do i = 1,maxNinstance
do k = 1,maxNinstance
homogenization_isostrain_sizeState(i) = 0_pInt
do j = 1,maxval(homogenization_Noutput)
do j = 1_pInt,maxval(homogenization_Noutput)
select case(homogenization_isostrain_output(j,i))
case('ngrains')
mySize = 1
mySize = 1_pInt
case default
mySize = 0
mySize = 0_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
@ -180,7 +181,7 @@ subroutine homogenization_isostrain_partitionDeformation(&
el & ! my element
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains,homogenization_Ngrains
implicit none
@ -193,7 +194,7 @@ subroutine homogenization_isostrain_partitionDeformation(&
integer(pInt) i
! homID = homogenization_typeInstance(mesh_element(3,el))
forall (i = 1:homogenization_Ngrains(mesh_element(3,el))) &
forall (i = 1_pInt:homogenization_Ngrains(mesh_element(3,el))) &
F(1:3,1:3,i) = avgF
return
@ -215,7 +216,6 @@ function homogenization_isostrain_updateState(&
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains
implicit none
@ -249,7 +249,7 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(&
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains, homogenization_Ngrains
implicit none
@ -263,8 +263,8 @@ subroutine homogenization_isostrain_averageStressAndItsTangent(&
! homID = homogenization_typeInstance(mesh_element(3,el))
Ngrains = homogenization_Ngrains(mesh_element(3,el))
avgP = sum(P,3)/dble(Ngrains)
dAvgPdAvgF = sum(dPdF,5)/dble(Ngrains)
avgP = sum(P,3)/real(Ngrains,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal)
return
@ -281,7 +281,7 @@ function homogenization_isostrain_averageTemperature(&
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use mesh, only: mesh_element
use material, only: homogenization_maxNgrains, homogenization_Ngrains
implicit none
@ -293,7 +293,7 @@ function homogenization_isostrain_averageTemperature(&
! homID = homogenization_typeInstance(mesh_element(3,el))
Ngrains = homogenization_Ngrains(mesh_element(3,el))
homogenization_isostrain_averageTemperature = sum(Temperature(1:Ngrains))/dble(Ngrains)
homogenization_isostrain_averageTemperature = sum(Temperature(1:Ngrains))/real(Ngrains,pReal)
return
@ -325,11 +325,11 @@ pure function homogenization_isostrain_postResults(&
c = 0_pInt
homogenization_isostrain_postResults = 0.0_pReal
do o = 1,homogenization_Noutput(mesh_element(3,el))
do o = 1_pInt,homogenization_Noutput(mesh_element(3,el))
select case(homogenization_isostrain_output(o,homID))
case ('ngrains')
homogenization_isostrain_postResults(c+1) = real(homogenization_isostrain_Ngrains(homID),pReal)
c = c + 1
homogenization_isostrain_postResults(c+1_pInt) = real(homogenization_isostrain_Ngrains(homID),pReal)
c = c + 1_pInt
end select
enddo