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

polishing, adding _pInt etc. where applicable 

post_results now handels zero increment different (like FEM, it is always there even if it is not counted)
This commit is contained in:
Martin Diehl 2012-02-13 14:18:07 +00:00
parent f03e7c459c
commit 156ec4582a
13 changed files with 643 additions and 615 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
code/DAMASK_marc.f90
View File

@ -1,4 +1,4 @@
! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
!
! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit.

2
code/DAMASK_spectral.f90
View File

@ -132,7 +132,7 @@ program DAMASK_spectral
0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal],&
[ 6_pInt, 6_pInt])
real(pReal), dimension(3,3,3,3) :: temp_3333 = 0.0_pReal
integer(pInt) :: size_reduced = 0.0_pReal ! number of stress BCs
integer(pInt) :: size_reduced = 0_pInt ! number of stress BCs
!--------------------------------------------------------------------------------------------------
! pointwise data

2
code/DAMASK_spectral_interface.f90
View File

@ -35,7 +35,7 @@ CONTAINS
!
!********************************************************************
subroutine DAMASK_interface_init()
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
implicit none
character(len=1024) commandLine, hostName, userName

6
code/constitutive.f90
View File

@ -101,7 +101,7 @@ logical knownConstitution
! --- PARSE CONSTITUTIONS FROM CONFIG FILE ---
if (.not. IO_open_jobFile(fileunit,material_localFileExt)) then ! no local material configuration present...
if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100) ! ...and cannot open material.config file
if (.not. IO_open_file(fileunit,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file
endif
call constitutive_j2_init(fileunit)
call constitutive_phenopowerlaw_init(fileunit)
@ -114,7 +114,7 @@ close(fileunit)
! --- WRITE DESCRIPTION FILE FOR CONSTITUTIVE PHASE OUTPUT ---
if(.not. IO_write_jobFile(fileunit,'outputConstitutive')) then ! problems in writing file
call IO_error (50)
call IO_error (50_pInt)
endif
do p = 1,material_Nphase
i = phase_constitutionInstance(p) ! which instance of a constitution is present phase
@ -323,7 +323,7 @@ endif
constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(myInstance)
case default
call IO_error(200,material_phase(g,i,e)) ! unknown constitution
call IO_error(200_pInt,material_phase(g,i,e)) ! unknown constitution
end select
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p

263
code/constitutive_dislotwin.f90
View File

@ -287,148 +287,148 @@ do ! read thru sections of
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
section = section + 1_pInt
output = 0_pInt ! reset output counter
endif
if (section > 0 .and. phase_constitution(section) == constitutive_dislotwin_label) then ! one of my sections
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
select case(tag)
case ('(output)')
output = output + 1
constitutive_dislotwin_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
output = output + 1_pInt
constitutive_dislotwin_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('lattice_structure')
constitutive_dislotwin_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
constitutive_dislotwin_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('covera_ratio')
constitutive_dislotwin_CoverA(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_CoverA(i) = IO_floatValue(line,positions,2_pInt)
case ('c11')
constitutive_dislotwin_C11(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_C11(i) = IO_floatValue(line,positions,2_pInt)
case ('c12')
constitutive_dislotwin_C12(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_C12(i) = IO_floatValue(line,positions,2_pInt)
case ('c13')
constitutive_dislotwin_C13(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_C13(i) = IO_floatValue(line,positions,2_pInt)
case ('c33')
constitutive_dislotwin_C33(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_C33(i) = IO_floatValue(line,positions,2_pInt)
case ('c44')
constitutive_dislotwin_C44(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_C44(i) = IO_floatValue(line,positions,2_pInt)
case ('nslip')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_dislotwin_Nslip(j,i) = IO_intValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_dislotwin_Nslip(j,i) = IO_intValue(line,positions,1_pInt+j)
case ('ntwin')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_dislotwin_Ntwin(j,i) = IO_intValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_dislotwin_Ntwin(j,i) = IO_intValue(line,positions,1_pInt+j)
case ('rhoedge0')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_dislotwin_rhoEdge0(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_dislotwin_rhoEdge0(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('rhoedgedip0')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_dislotwin_rhoEdgeDip0(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_dislotwin_rhoEdgeDip0(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('slipburgers')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_dislotwin_burgersPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_dislotwin_burgersPerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('twinburgers')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_dislotwin_burgersPerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_dislotwin_burgersPerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('qedge')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_dislotwin_QedgePerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_dislotwin_QedgePerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('v0')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_dislotwin_v0PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_dislotwin_v0PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('ndot0')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_dislotwin_Ndot0PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_dislotwin_Ndot0PerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('twinsize')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_dislotwin_twinsizePerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_dislotwin_twinsizePerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('clambdaslip')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_dislotwin_CLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_dislotwin_CLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('grainsize')
constitutive_dislotwin_GrainSize(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_GrainSize(i) = IO_floatValue(line,positions,2_pInt)
case ('maxtwinfraction')
constitutive_dislotwin_MaxTwinFraction(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_MaxTwinFraction(i) = IO_floatValue(line,positions,2_pInt)
case ('pexponent')
constitutive_dislotwin_p(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_p(i) = IO_floatValue(line,positions,2_pInt)
case ('qexponent')
constitutive_dislotwin_q(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_q(i) = IO_floatValue(line,positions,2_pInt)
case ('rexponent')
constitutive_dislotwin_r(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_r(i) = IO_floatValue(line,positions,2_pInt)
case ('d0')
constitutive_dislotwin_D0(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_D0(i) = IO_floatValue(line,positions,2_pInt)
case ('qsd')
constitutive_dislotwin_Qsd(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_Qsd(i) = IO_floatValue(line,positions,2_pInt)
case ('atol_rho')
constitutive_dislotwin_aTolRho(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_aTolRho(i) = IO_floatValue(line,positions,2_pInt)
case ('cmfptwin')
constitutive_dislotwin_Cmfptwin(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_Cmfptwin(i) = IO_floatValue(line,positions,2_pInt)
case ('cthresholdtwin')
constitutive_dislotwin_Cthresholdtwin(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_Cthresholdtwin(i) = IO_floatValue(line,positions,2_pInt)
case ('solidsolutionstrength')
constitutive_dislotwin_SolidSolutionStrength(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_SolidSolutionStrength(i) = IO_floatValue(line,positions,2_pInt)
case ('l0')
constitutive_dislotwin_L0(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_L0(i) = IO_floatValue(line,positions,2_pInt)
case ('cedgedipmindistance')
constitutive_dislotwin_CEdgeDipMinDistance(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_CEdgeDipMinDistance(i) = IO_floatValue(line,positions,2_pInt)
case ('catomicvolume')
constitutive_dislotwin_CAtomicVolume(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_CAtomicVolume(i) = IO_floatValue(line,positions,2_pInt)
case ('interactionslipslip')
forall (j = 1:lattice_maxNinteraction) &
constitutive_dislotwin_interactionSlipSlip(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_dislotwin_interactionSlipSlip(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('interactionsliptwin')
forall (j = 1:lattice_maxNinteraction) &
constitutive_dislotwin_interactionSlipTwin(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_dislotwin_interactionSlipTwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('interactiontwinslip')
forall (j = 1:lattice_maxNinteraction) &
constitutive_dislotwin_interactionTwinSlip(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_dislotwin_interactionTwinSlip(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('interactiontwintwin')
forall (j = 1:lattice_maxNinteraction) &
constitutive_dislotwin_interactionTwinTwin(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_dislotwin_interactionTwinTwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('sfe_0k')
constitutive_dislotwin_SFE_0K(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_SFE_0K(i) = IO_floatValue(line,positions,2_pInt)
case ('dsfe_dt')
constitutive_dislotwin_dSFE_dT(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_dSFE_dT(i) = IO_floatValue(line,positions,2_pInt)
case ('shearbandresistance')
constitutive_dislotwin_sbResistance(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_sbResistance(i) = IO_floatValue(line,positions,2_pInt)
case ('shearbandvelocity')
constitutive_dislotwin_sbVelocity(i) = IO_floatValue(line,positions,2)
constitutive_dislotwin_sbVelocity(i) = IO_floatValue(line,positions,2_pInt)
end select
endif
enddo
100 do i = 1,maxNinstance
100 do i = 1_pInt,maxNinstance
constitutive_dislotwin_structure(i) = &
lattice_initializeStructure(constitutive_dislotwin_structureName(i),constitutive_dislotwin_CoverA(i))
myStructure = constitutive_dislotwin_structure(i)
!* Sanity checks
if (myStructure < 1 .or. myStructure > 3) call IO_error(205)
if (sum(constitutive_dislotwin_Nslip(:,i)) <= 0_pInt) call IO_error(225)
if (sum(constitutive_dislotwin_Ntwin(:,i)) < 0_pInt) call IO_error(225) !***
if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt)
if (sum(constitutive_dislotwin_Nslip(:,i)) <= 0_pInt) call IO_error(225_pInt)
if (sum(constitutive_dislotwin_Ntwin(:,i)) < 0_pInt) call IO_error(225_pInt) !***
do f = 1,lattice_maxNslipFamily
if (constitutive_dislotwin_Nslip(f,i) > 0_pInt) then
if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(220)
if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(220)
if (constitutive_dislotwin_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(221)
if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226)
if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(220_pInt)
if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(220_pInt)
if (constitutive_dislotwin_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(221_pInt)
if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt)
endif
enddo
do f = 1,lattice_maxNtwinFamily
if (constitutive_dislotwin_Ntwin(f,i) > 0_pInt) then
if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221) !***
if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(226) !***
if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221_pInt) !***
if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(226_pInt) !***
endif
enddo
! if (any(constitutive_dislotwin_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) call IO_error(229)
if (constitutive_dislotwin_CAtomicVolume(i) <= 0.0_pReal) call IO_error(230)
if (constitutive_dislotwin_D0(i) <= 0.0_pReal) call IO_error(231)
if (constitutive_dislotwin_Qsd(i) <= 0.0_pReal) call IO_error(232)
if (constitutive_dislotwin_aTolRho(i) <= 0.0_pReal) call IO_error(233)
if (constitutive_dislotwin_sbResistance(i) <= 0.0_pReal) call IO_error(234)
if (constitutive_dislotwin_sbVelocity(i) < 0.0_pReal) call IO_error(235)
if (constitutive_dislotwin_SFE_0K(i) == 0.0_pReal .AND. &
constitutive_dislotwin_dSFE_dT(i) == 0.0_pReal) call IO_error(223)
if (constitutive_dislotwin_CAtomicVolume(i) <= 0.0_pReal) call IO_error(230_pInt)
if (constitutive_dislotwin_D0(i) <= 0.0_pReal) call IO_error(231_pInt)
if (constitutive_dislotwin_Qsd(i) <= 0.0_pReal) call IO_error(232_pInt)
if (constitutive_dislotwin_aTolRho(i) <= 0.0_pReal) call IO_error(233_pInt)
if (constitutive_dislotwin_sbResistance(i) <= 0.0_pReal) call IO_error(234_pInt)
if (constitutive_dislotwin_sbVelocity(i) < 0.0_pReal) call IO_error(235_pInt)
if (constitutive_dislotwin_SFE_0K(i) == 0.0_pReal .and. &
constitutive_dislotwin_dSFE_dT(i) == 0.0_pReal) call IO_error(223_pInt)
!* Determine total number of active slip or twin systems
constitutive_dislotwin_Nslip(:,i) = min(lattice_NslipSystem(:,myStructure),constitutive_dislotwin_Nslip(:,i))
@ -473,25 +473,25 @@ 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
do i = 1,maxNinstance
do i = 1_pInt,maxNinstance
myStructure = constitutive_dislotwin_structure(i)
!* Inverse lookup of my slip system family
l = 0_pInt
do f = 1,lattice_maxNslipFamily
do k = 1,constitutive_dislotwin_Nslip(f,i)
l = l + 1
do f = 1_pInt,lattice_maxNslipFamily
do k = 1_pInt,constitutive_dislotwin_Nslip(f,i)
l = l + 1_pInt
constitutive_dislotwin_slipFamily(l,i) = f
constitutive_dislotwin_slipSystemLattice(l,i) = sum(lattice_NslipSystem(1:f-1,myStructure)) + k
constitutive_dislotwin_slipSystemLattice(l,i) = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) + k
enddo; enddo
!* Inverse lookup of my twin system family
l = 0_pInt
do f = 1,lattice_maxNtwinFamily
do k = 1,constitutive_dislotwin_Ntwin(f,i)
l = l + 1
do f = 1_pInt,lattice_maxNtwinFamily
do k = 1_pInt,constitutive_dislotwin_Ntwin(f,i)
l = l + 1_pInt
constitutive_dislotwin_twinFamily(l,i) = f
constitutive_dislotwin_twinSystemLattice(l,i) = sum(lattice_NtwinSystem(1:f-1,myStructure)) + k
constitutive_dislotwin_twinSystemLattice(l,i) = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) + k
enddo; enddo
!* Determine size of state array
@ -504,7 +504,7 @@ do i = 1,maxNinstance
size(constitutive_dislotwin_listDependentSlipStates)*ns+size(constitutive_dislotwin_listDependentTwinStates)*nt
!* Determine size of postResults array
do o = 1,maxval(phase_Noutput)
do o = 1_pInt,maxval(phase_Noutput)
select case(constitutive_dislotwin_output(o,i))
case('edge_density', &
'dipole_density', &
@ -545,14 +545,14 @@ do i = 1,maxNinstance
!* Elasticity matrix and shear modulus according to material.config
select case (myStructure)
case(1:2) ! cubic(s)
forall(k=1:3)
forall(j=1:3) &
case(1_pInt:2_pInt) ! cubic(s)
forall(k=1_pInt:3_pInt)
forall(j=1_pInt:3_pInt) &
constitutive_dislotwin_Cslip_66(k,j,i) = constitutive_dislotwin_C12(i)
constitutive_dislotwin_Cslip_66(k,k,i) = constitutive_dislotwin_C11(i)
constitutive_dislotwin_Cslip_66(k+3,k+3,i) = constitutive_dislotwin_C44(i)
constitutive_dislotwin_Cslip_66(k+3_pInt,k+3_pInt,i) = constitutive_dislotwin_C44(i)
end forall
case(3:) ! all hex
case(3_pInt:) ! all hex
constitutive_dislotwin_Cslip_66(1,1,i) = constitutive_dislotwin_C11(i)
constitutive_dislotwin_Cslip_66(2,2,i) = constitutive_dislotwin_C11(i)
constitutive_dislotwin_Cslip_66(3,3,i) = constitutive_dislotwin_C33(i)
@ -572,23 +572,24 @@ do i = 1,maxNinstance
0.2_pReal*(constitutive_dislotwin_C11(i)-constitutive_dislotwin_C12(i))+0.3_pReal*constitutive_dislotwin_C44(i)
!* Construction of the twin elasticity matrices
do j=1,lattice_maxNtwinFamily
do k=1,constitutive_dislotwin_Ntwin(j,i)
do l=1,3 ; do m=1,3 ; do n=1,3 ; do o=1,3 ; do p=1,3 ; do q=1,3 ; do r=1,3 ; do s=1,3
constitutive_dislotwin_Ctwin_3333(l,m,n,o,sum(constitutive_dislotwin_Nslip(1:j-1,i))+k,i) = &
constitutive_dislotwin_Ctwin_3333(l,m,n,o,sum(constitutive_dislotwin_Nslip(1:j-1,i))+k,i) + &
do j=1_pInt,lattice_maxNtwinFamily
do k=1_pInt,constitutive_dislotwin_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_dislotwin_Ctwin_3333(l,m,n,o,sum(constitutive_dislotwin_Nslip(1:j-1_pInt,i))+k,i) = &
constitutive_dislotwin_Ctwin_3333(l,m,n,o,sum(constitutive_dislotwin_Nslip(1:j-1_pInt,i))+k,i) + &
constitutive_dislotwin_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_dislotwin_Ctwin_66(:,:,k,i) = math_Mandel3333to66(constitutive_dislotwin_Ctwin_3333(:,:,:,:,k,i))
enddo
enddo
!* Burgers vector, dislocation velocity prefactor, mean free path prefactor and minimum dipole distance for each slip system
do s = 1,constitutive_dislotwin_totalNslip(i)
do s = 1_pInt,constitutive_dislotwin_totalNslip(i)
f = constitutive_dislotwin_slipFamily(s,i)
constitutive_dislotwin_burgersPerSlipSystem(s,i) = constitutive_dislotwin_burgersPerSlipFamily(f,i)
constitutive_dislotwin_QedgePerSlipSystem(s,i) = constitutive_dislotwin_QedgePerSlipFamily(f,i)
@ -597,7 +598,7 @@ do i = 1,maxNinstance
enddo
!* Burgers vector, nucleation rate prefactor and twin size for each twin system
do s = 1,constitutive_dislotwin_totalNtwin(i)
do s = 1_pInt,constitutive_dislotwin_totalNtwin(i)
f = constitutive_dislotwin_twinFamily(s,i)
constitutive_dislotwin_burgersPerTwinSystem(s,i) = constitutive_dislotwin_burgersPerTwinFamily(f,i)
constitutive_dislotwin_Ndot0PerTwinSystem(s,i) = constitutive_dislotwin_Ndot0PerTwinFamily(f,i)
@ -605,41 +606,41 @@ do i = 1,maxNinstance
enddo
!* Construction of interaction matrices
do s1 = 1,constitutive_dislotwin_totalNslip(i)
do s2 = 1,constitutive_dislotwin_totalNslip(i)
do s1 = 1_pInt,constitutive_dislotwin_totalNslip(i)
do s2 = 1_pInt,constitutive_dislotwin_totalNslip(i)
constitutive_dislotwin_interactionMatrixSlipSlip(s1,s2,i) = &
constitutive_dislotwin_interactionSlipSlip(lattice_interactionSlipSlip(constitutive_dislotwin_slipSystemLattice(s1,i), &
constitutive_dislotwin_slipSystemLattice(s2,i), &
myStructure),i)
enddo; enddo
do s1 = 1,constitutive_dislotwin_totalNslip(i)
do t2 = 1,constitutive_dislotwin_totalNtwin(i)
do s1 = 1_pInt,constitutive_dislotwin_totalNslip(i)
do t2 = 1_pInt,constitutive_dislotwin_totalNtwin(i)
constitutive_dislotwin_interactionMatrixSlipTwin(s1,t2,i) = &
constitutive_dislotwin_interactionSlipTwin(lattice_interactionSlipTwin(constitutive_dislotwin_slipSystemLattice(s1,i), &
constitutive_dislotwin_twinSystemLattice(t2,i), &
myStructure),i)
constitutive_dislotwin_interactionSlipTwin(&
lattice_interactionSlipTwin(constitutive_dislotwin_slipSystemLattice(s1,i), &
constitutive_dislotwin_twinSystemLattice(t2,i),myStructure),i)
enddo; enddo
do t1 = 1,constitutive_dislotwin_totalNtwin(i)
do s2 = 1,constitutive_dislotwin_totalNslip(i)
do t1 = 1_pInt,constitutive_dislotwin_totalNtwin(i)
do s2 = 1_pInt,constitutive_dislotwin_totalNslip(i)
constitutive_dislotwin_interactionMatrixTwinSlip(t1,s2,i) = &
constitutive_dislotwin_interactionTwinSlip(lattice_interactionTwinSlip(constitutive_dislotwin_twinSystemLattice(t1,i), &
constitutive_dislotwin_slipSystemLattice(s2,i), &
myStructure),i)
constitutive_dislotwin_interactionTwinSlip(lattice_interactionTwinSlip(&
constitutive_dislotwin_twinSystemLattice(t1,i), &
constitutive_dislotwin_slipSystemLattice(s2,i), myStructure),i)
enddo; enddo
do t1 = 1,constitutive_dislotwin_totalNtwin(i)
do t2 = 1,constitutive_dislotwin_totalNtwin(i)
do t1 = 1_pInt,constitutive_dislotwin_totalNtwin(i)
do t2 = 1_pInt,constitutive_dislotwin_totalNtwin(i)
constitutive_dislotwin_interactionMatrixTwinTwin(t1,t2,i) = &
constitutive_dislotwin_interactionTwinTwin(lattice_interactionTwinTwin(constitutive_dislotwin_twinSystemLattice(t1,i), &
constitutive_dislotwin_twinSystemLattice(t2,i), &
myStructure),i)
constitutive_dislotwin_interactionTwinTwin(&
lattice_interactionTwinTwin(constitutive_dislotwin_twinSystemLattice(t1,i), &
constitutive_dislotwin_twinSystemLattice(t2,i), myStructure),i)
enddo; enddo
!* Calculation of forest projections for edge dislocations
do s1 = 1,constitutive_dislotwin_totalNslip(i)
do s2 = 1,constitutive_dislotwin_totalNslip(i)
do s1 = 1_pInt,constitutive_dislotwin_totalNslip(i)
do s2 = 1_pInt,constitutive_dislotwin_totalNslip(i)
constitutive_dislotwin_forestProjectionEdge(s1,s2,i) = &
abs(math_mul3x3(lattice_sn(:,constitutive_dislotwin_slipSystemLattice(s1,i),myStructure), &
lattice_st(:,constitutive_dislotwin_slipSystemLattice(s2,i),myStructure)))
@ -678,7 +679,7 @@ constitutive_dislotwin_stateInit = 0.0_pReal
!* Initialize basic slip state variables
s1 = 0_pInt
do f = 1,lattice_maxNslipFamily
do f = 1_pInt,lattice_maxNslipFamily
s0 = s1 + 1_pInt
s1 = s0 + constitutive_dislotwin_Nslip(f,myInstance) - 1_pInt
do s = s0,s1
@ -687,24 +688,24 @@ do f = 1,lattice_maxNslipFamily
enddo
enddo
constitutive_dislotwin_stateInit(1:ns) = rhoEdge0
constitutive_dislotwin_stateInit(ns+1:2*ns) = rhoEdgeDip0
constitutive_dislotwin_stateInit(ns+1:2_pInt*ns) = rhoEdgeDip0
!* Initialize dependent slip microstructural variables
forall (s = 1:ns) &
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*ns+nt+1:3*ns+nt) = invLambdaSlip0
constitutive_dislotwin_stateInit(2_pInt*ns+nt+1:3_pInt*ns+nt) = invLambdaSlip0
forall (s = 1:ns) &
forall (s = 1_pInt:ns) &
MeanFreePathSlip0(s) = &
constitutive_dislotwin_GrainSize(myInstance)/(1.0_pReal+invLambdaSlip0(s)*constitutive_dislotwin_GrainSize(myInstance))
constitutive_dislotwin_stateInit(4*ns+2*nt+1:5*ns+2*nt) = MeanFreePathSlip0
constitutive_dislotwin_stateInit(4_pInt*ns+2_pInt*nt+1:5_pInt*ns+2_pInt*nt) = MeanFreePathSlip0
forall (s = 1:ns) &
forall (s = 1_pInt:ns) &
tauSlipThreshold0(s) = constitutive_dislotwin_SolidSolutionStrength(myInstance)+ &
constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)* &
sqrt(dot_product((rhoEdge0+rhoEdgeDip0),constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance)))
constitutive_dislotwin_stateInit(5*ns+3*nt+1:6*ns+3*nt) = tauSlipThreshold0
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) &

79
code/constitutive_j2.f90
View File

@ -83,14 +83,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+2_pInt*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k, output, mySize
character(len=64) tag
character(len=1024) line
@ -103,9 +104,9 @@ subroutine constitutive_j2_init(file)
!$OMP END CRITICAL (write2out)
maxNinstance = count(phase_constitution == constitutive_j2_label)
if (maxNinstance == 0) return
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,*)
@ -131,7 +132,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 <phase>
read(file,'(a1024)',END=100) line
@ -142,53 +143,53 @@ subroutine constitutive_j2_init(file)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
section = section + 1_pInt
output = 0_pInt ! reset output counter
endif
if (section > 0 .and. phase_constitution(section) == constitutive_j2_label) then ! one of my sections
if (section > 0_pInt .and. phase_constitution(section) == constitutive_j2_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
select case(tag)
case ('(output)')
output = output + 1
constitutive_j2_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
output = output + 1_pInt
constitutive_j2_output(output,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)
end select
endif
enddo
100 do i = 1,maxNinstance ! sanity checks
if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(210)
if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211)
if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(212)
if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(213)
if (constitutive_j2_a(i) <= 0.0_pReal) call IO_error(241)
if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(240)
if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(242)
100 do i = 1_pInt,maxNinstance ! sanity checks
if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(210_pInt)
if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211_pInt)
if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(212_pInt)
if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(213_pInt)
if (constitutive_j2_a(i) <= 0.0_pReal) call IO_error(241_pInt)
if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(240_pInt)
if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(242_pInt)
enddo
do i = 1,maxNinstance
do j = 1,maxval(phase_Noutput)
do i = 1_pInt,maxNinstance
do j = 1_pInt,maxval(phase_Noutput)
select case(constitutive_j2_output(j,i))
case('flowstress')
mySize = 1_pInt
@ -205,22 +206,24 @@ 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) &
forall(k=1_pInt:3_pInt)
forall(j=1_pInt:3_pInt) &
constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i)
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))
constitutive_j2_Cslip_66(k+3_pInt,k+3_pInt,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i))
end forall
constitutive_j2_Cslip_66(1:6,1:6,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(1:6,1:6,i)))
enddo
return
endsubroutine
@ -241,6 +244,7 @@ pure function constitutive_j2_stateInit(myInstance)
endfunction
!*********************************************************************
!* relevant microstructural state *
!*********************************************************************
@ -261,6 +265,7 @@ real(pReal), dimension(constitutive_j2_sizeState(myInstance)) :: &
constitutive_j2_aTolState = constitutive_j2_aTolResistance(myInstance)
endfunction
@ -288,6 +293,7 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el)
return
endfunction
@ -310,6 +316,7 @@ subroutine constitutive_j2_microstructure(Temperature,state,ipc,ip,el)
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
endsubroutine
@ -335,6 +342,7 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v,
material_phase, &
phase_constitutionInstance
implicit none
!*** input variables ***!
@ -458,6 +466,7 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el)
endfunction
!****************************************************************
!* calculates the rate of change of temperature *
!****************************************************************

186
code/constitutive_nonlocal.f90
View File

@ -140,6 +140,7 @@ CONTAINS
!**************************************
subroutine constitutive_nonlocal_init(file)
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, &
@ -328,7 +329,7 @@ constitutive_nonlocal_peierlsStressPerSlipFamily = 0.0_pReal
rewind(file)
line = ''
section = 0
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
@ -339,98 +340,108 @@ do
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
section = section + 1_pInt
output = 0_pInt ! reset output counter
cycle
endif
if (section > 0 .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections
if (section > 0_pInt .and. phase_constitution(section) == constitutive_nonlocal_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','/nonlocal/')
cycle
case ('(output)')
output = output + 1
constitutive_nonlocal_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
output = output + 1_pInt
constitutive_nonlocal_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('lattice_structure')
constitutive_nonlocal_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
constitutive_nonlocal_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('c/a_ratio','covera_ratio')
constitutive_nonlocal_CoverA(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_CoverA(i) = IO_floatValue(line,positions,2_pInt)
case ('c11')
constitutive_nonlocal_C11(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_C11(i) = IO_floatValue(line,positions,2_pInt)
case ('c12')
constitutive_nonlocal_C12(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_C12(i) = IO_floatValue(line,positions,2_pInt)
case ('c13')
constitutive_nonlocal_C13(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_C13(i) = IO_floatValue(line,positions,2_pInt)
case ('c33')
constitutive_nonlocal_C33(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_C33(i) = IO_floatValue(line,positions,2_pInt)
case ('c44')
constitutive_nonlocal_C44(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_C44(i) = IO_floatValue(line,positions,2_pInt)
case ('nslip')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_Nslip(f,i) = IO_intValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_Nslip(f,i)&
= IO_intValue(line,positions,1_pInt+f)
case ('rhosgledgepos0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgePos0(f,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgePos0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
case ('rhosgledgeneg0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgeNeg0(f,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgeNeg0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
case ('rhosglscrewpos0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewPos0(f,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewPos0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
case ('rhosglscrewneg0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewNeg0(f,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewNeg0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
case ('rhodipedge0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipEdge0(f,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipEdge0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
case ('rhodipscrew0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipScrew0(f,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipScrew0(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
case ('lambda0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_lambda0PerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_lambda0PerSlipFamily(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
case ('burgers')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_burgersPerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_burgersPerSlipFamily(f,i)&
= IO_floatValue(line,positions,1_pInt+f)
case('cutoffradius','r')
constitutive_nonlocal_R(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_R(i) = IO_floatValue(line,positions,2_pInt)
case('minimumdipoleheightedge','ddipminedge')
forall (f = 1:lattice_maxNslipFamily) &
constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,1,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,1_pInt,i)&
= IO_floatValue(line,positions,1_pInt+f)
case('minimumdipoleheightscrew','ddipminscrew')
forall (f = 1:lattice_maxNslipFamily) &
constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,2,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_minimumDipoleHeightPerSlipFamily(f,2_pInt,i)&
= IO_floatValue(line,positions,1_pInt+f)
case('atomicvolume')
constitutive_nonlocal_atomicVolume(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_atomicVolume(i) = IO_floatValue(line,positions,2_pInt)
case('selfdiffusionprefactor','dsd0')
constitutive_nonlocal_Dsd0(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_Dsd0(i) = IO_floatValue(line,positions,2_pInt)
case('selfdiffusionenergy','qsd')
constitutive_nonlocal_Qsd(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_Qsd(i) = IO_floatValue(line,positions,2_pInt)
case('atol_rho')
constitutive_nonlocal_aTolRho(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_aTolRho(i) = IO_floatValue(line,positions,2_pInt)
case ('interaction_slipslip')
forall (it = 1:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(it,i) = IO_floatValue(line,positions,1+it)
forall (it = 1_pInt:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(it,i)&
= IO_floatValue(line,positions,1_pInt+it)
case('peierlsstressedge')
forall (f = 1:lattice_maxNslipFamily) &
constitutive_nonlocal_peierlsStressPerSlipFamily(f,1,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_peierlsStressPerSlipFamily(f,1_pInt,i)&
= IO_floatValue(line,positions,1_pInt+f)
case('peierlsstressscrew')
forall (f = 1:lattice_maxNslipFamily) &
constitutive_nonlocal_peierlsStressPerSlipFamily(f,2,i) = IO_floatValue(line,positions,1+f)
forall (f = 1_pInt:lattice_maxNslipFamily) constitutive_nonlocal_peierlsStressPerSlipFamily(f,2_pInt,i)&
= IO_floatValue(line,positions,1_pInt+f)
case('doublekinkwidth')
constitutive_nonlocal_doublekinkwidth(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_doublekinkwidth(i) = IO_floatValue(line,positions,2_pInt)
case('solidsolutionenergy')
constitutive_nonlocal_solidSolutionEnergy(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_solidSolutionEnergy(i) = IO_floatValue(line,positions,2_pInt)
case('solidsolutionsize')
constitutive_nonlocal_solidSolutionSize(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_solidSolutionSize(i) = IO_floatValue(line,positions,2_pInt)
case('solidsolutionconcentration')
constitutive_nonlocal_solidSolutionConcentration(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_solidSolutionConcentration(i) = IO_floatValue(line,positions,2_pInt)
case('p')
constitutive_nonlocal_p(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_p(i) = IO_floatValue(line,positions,2_pInt)
case('q')
constitutive_nonlocal_q(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_q(i) = IO_floatValue(line,positions,2_pInt)
case('viscosity','glideviscosity')
constitutive_nonlocal_viscosity(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_viscosity(i) = IO_floatValue(line,positions,2_pInt)
case('attackfrequency','fattack')
constitutive_nonlocal_fattack(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_fattack(i) = IO_floatValue(line,positions,2_pInt)
case('rhosglscatter')
constitutive_nonlocal_rhoSglScatter(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_rhoSglScatter(i) = IO_floatValue(line,positions,2_pInt)
case('surfacetransmissivity')
constitutive_nonlocal_surfaceTransmissivity(i) = IO_floatValue(line,positions,2)
constitutive_nonlocal_surfaceTransmissivity(i) = IO_floatValue(line,positions,2_pInt)
case default
call IO_error(236,ext_msg=tag)
call IO_error(236_pInt,ext_msg=tag)
end select
endif
enddo
@ -651,13 +662,13 @@ do i = 1,maxNinstance
!*** elasticity matrix and shear modulus according to material.config
select case (myStructure)
case(1:2) ! cubic(s)
forall(k=1:3)
forall(j=1:3) constitutive_nonlocal_Cslip_66(k,j,i) = constitutive_nonlocal_C12(i)
case(1_pInt:2_pInt) ! cubic(s)
forall(k=1_pInt:3_pInt)
forall(j=1_pInt:3_pInt) constitutive_nonlocal_Cslip_66(k,j,i) = constitutive_nonlocal_C12(i)
constitutive_nonlocal_Cslip_66(k,k,i) = constitutive_nonlocal_C11(i)
constitutive_nonlocal_Cslip_66(k+3,k+3,i) = constitutive_nonlocal_C44(i)
constitutive_nonlocal_Cslip_66(k+3_pInt,k+3_pInt,i) = constitutive_nonlocal_C44(i)
end forall
case(3:) ! all hex
case(3_pInt:) ! all hex
constitutive_nonlocal_Cslip_66(1,1,i) = constitutive_nonlocal_C11(i)
constitutive_nonlocal_Cslip_66(2,2,i) = constitutive_nonlocal_C11(i)
constitutive_nonlocal_Cslip_66(3,3,i) = constitutive_nonlocal_C33(i)
@ -681,7 +692,7 @@ do i = 1,maxNinstance
/ ( 4.0_pReal*constitutive_nonlocal_C11(i) + 6.0_pReal*constitutive_nonlocal_C12(i) &
+ 2.0_pReal*constitutive_nonlocal_C44(i) ) ! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
do s1 = 1,ns
do s1 = 1_pInt,ns
f = constitutive_nonlocal_slipFamily(s1,i)
!*** burgers vector, mean free path prefactor and minimum dipole distance for each slip system
@ -708,17 +719,16 @@ do i = 1,maxNinstance
constitutive_nonlocal_interactionMatrixSlipSlip(s1,s2,i) &
= constitutive_nonlocal_interactionSlipSlip(lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), &
constitutive_nonlocal_slipSystemLattice(s2,i), &
myStructure), &
i)
myStructure), i)
enddo
!*** rotation matrix from lattice configuration to slip system
constitutive_nonlocal_lattice2slip(1:3,1:3,s1,i) &
= math_transpose33( reshape((/ lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
= math_transpose33( reshape([ lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
-lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure)/), (/3,3/)))
lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure)], [3,3]))
enddo
enddo
@ -1058,9 +1068,9 @@ if (.not. phase_localConstitution(phase)) then
if (neighboring_ns == ns) then
if (neighboring_el /= el .or. neighboring_ip /= ip) then
connection_latticeConf(1:3,n) = math_mul33x3(invFe, neighboring_ipCoords - ipCoords)
forall (s = 1:ns, c = 1:2) &
neighboring_rhoExcess(c,s,n) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*ns+s) & ! positive mobiles
- state(g,neighboring_ip,neighboring_el)%p((2*c-1)*ns+s) ! negative mobiles
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
neighboring_rhoExcess(c,s,n) = state(g,neighboring_ip,neighboring_el)%p((2_pInt*c-2_pInt)*ns+s) & ! positive mobiles
- state(g,neighboring_ip,neighboring_el)%p((2_pInt*c-1_pInt)*ns+s) ! negative mobiles
else
! thats myself! probably using periodic images
connection_latticeConf(1:3,n) = areaNormal_latticeConf(1:3,n) * FVsize
@ -1068,7 +1078,7 @@ if (.not. phase_localConstitution(phase)) then
endif
else
! different number of active slip systems
call IO_error(-1,ext_msg='different number of active slip systems in neighboring IPs of same crystal structure')
call IO_error(-1_pInt,ext_msg='different number of active slip systems in neighboring IPs of same crystal structure')
endif
else
! local neighbor or different lattice structure or different constitution instance
@ -1449,15 +1459,15 @@ enddo
!*** get dislocation velocity and its tangent and store the velocity in the state array
if (myStructure == 1_pInt) then ! for fcc all velcities are equal
call constitutive_nonlocal_kinetics(v(1:ns,1), tau, 1, Temperature, state, g, ip, el, dv_dtau(1:ns,1))
do t = 1,4
call constitutive_nonlocal_kinetics(v(1:ns,1), tau, 1_pInt, Temperature, state, g, ip, el, dv_dtau(1:ns,1))
do t = 1_pInt,4_pInt
v(1:ns,t) = v(1:ns,1)
dv_dtau(1:ns,t) = dv_dtau(1:ns,1)
state%p((12+t)*ns+1:(13+t)*ns) = v(1:ns,1)
state%p((12_pInt+t)*ns+1:(13_pInt+t)*ns) = v(1:ns,1)
enddo
else ! for all other lattice structures the velcities may vary with character and sign
do t = 1,4
c = (t-1)/2+1
do t = 1_pInt,4_pInt
c = (t-1_pInt)/2_pInt+1_pInt
call constitutive_nonlocal_kinetics(v(1:ns,t), tau, c, Temperature, state, g, ip, el, dv_dtau(1:ns,t))
state%p((12+t)*ns+1:(13+t)*ns) = v(1:ns,t)
enddo
@ -1466,8 +1476,8 @@ endif
!*** Bauschinger effect
forall (s = 1:ns, t = 5:8, state%p((t-1)*ns+s) * v(s,t-4) < 0.0_pReal) &
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(state%p((t-1)*ns+s))
forall (s = 1_pInt:ns, t = 5_pInt:8_pInt, state%p((t-1)*ns+s) * v(s,t-4_pInt) < 0.0_pReal) &
rhoSgl(s,t-4_pInt) = rhoSgl(s,t-4_pInt) + abs(state%p((t-1_pInt)*ns+s))
!*** Calculation of gdot and its tangent
@ -1636,7 +1646,7 @@ logical considerEnteringFlux, &
considerLeavingFlux
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g
write(6,*)
@ -1644,10 +1654,10 @@ logical considerEnteringFlux, &
#endif
select case(mesh_element(2,el))
case (1,6,7,8,9)
case (1_pInt,6_pInt,7_pInt,8_pInt,9_pInt)
! all fine
case default
call IO_error(-1,el,ip,g,'element type not supported for nonlocal constitution')
call IO_error(-1_pInt,el,ip,g,'element type not supported for nonlocal constitution')
end select
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
@ -1662,17 +1672,17 @@ dUpper = 0.0_pReal
!*** shortcut to state variables
forall (s = 1:ns, t = 1:4) &
rhoSgl(s,t) = max(state(g,ip,el)%p((t-1)*ns+s), 0.0_pReal)
forall (s = 1:ns, t = 5:8) &
rhoSgl(s,t) = state(g,ip,el)%p((t-1)*ns+s)
forall (s = 1:ns, c = 1:2) &
rhoDip(s,c) = max(state(g,ip,el)%p((7+c)*ns+s), 0.0_pReal)
rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
tauBack = state(g,ip,el)%p(12*ns+1:13*ns)
forall (t = 1:4) &
v(1:ns,t) = state(g,ip,el)%p((12+t)*ns+1:(13+t)*ns)
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) &
rhoSgl(s,t) = max(state(g,ip,el)%p((t-1_pInt)*ns+s), 0.0_pReal)
forall (s = 1_pInt:ns, t = 5_pInt:8_pInt) &
rhoSgl(s,t) = state(g,ip,el)%p((t-1_pInt)*ns+s)
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
rhoDip(s,c) = max(state(g,ip,el)%p((7_pInt+c)*ns+s), 0.0_pReal)
rhoForest = state(g,ip,el)%p(10_pInt*ns+1:11_pInt*ns)
tauThreshold = state(g,ip,el)%p(11_pInt*ns+1_pInt:12_pInt*ns)
tauBack = state(g,ip,el)%p(12_pInt*ns+1:13_pInt*ns)
forall (t = 1_pInt:4_pInt) &
v(1_pInt:ns,t) = state(g,ip,el)%p((12_pInt+t)*ns+1_pInt:(13_pInt+t)*ns)
!*** sanity check for timestep
@ -1687,13 +1697,13 @@ endif
!****************************************************************************
!*** Calculate shear rate
forall (t = 1:4) &
gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgers(1:ns,myInstance) * v(1:ns,t)
forall (s = 1:ns, t = 1:4, rhoSgl(s,t+4) * v(s,t) < 0.0_pReal) & ! contribution of used rho for changing sign of v
forall (t = 1_pInt:4_pInt) &
gdot(1_pInt:ns,t) = rhoSgl(1_pInt:ns,t) * constitutive_nonlocal_burgers(1:ns,myInstance) * v(1:ns,t)
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt, rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pReal) & ! contribution of used rho for changing sign of v
gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgers(s,myInstance) * v(s,t)
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
if (debug_verbosity > 6_pInt .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot
endif
@ -1801,7 +1811,7 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then
enddo
endif
opposite_n = n + mod(n,2) - mod(n+1,2)
opposite_n = n + mod(n,2_pInt) - mod(n+1_pInt,2_pInt)
opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el)
opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el)

173
code/constitutive_phenopowerlaw.f90
View File

@ -263,88 +263,93 @@ subroutine constitutive_phenopowerlaw_init(file)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
section = section + 1_pInt
output = 0_pInt ! reset output counter
endif
if (section > 0 .and. phase_constitution(section) == constitutive_phenopowerlaw_label) then ! one of my sections
if (section > 0_pInt .and. phase_constitution(section) == constitutive_phenopowerlaw_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 ('(output)')
output = output + 1
constitutive_phenopowerlaw_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
output = output + 1_pInt
constitutive_phenopowerlaw_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('lattice_structure')
constitutive_phenopowerlaw_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
constitutive_phenopowerlaw_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('covera_ratio')
constitutive_phenopowerlaw_CoverA(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_CoverA(i) = IO_floatValue(line,positions,2_pInt)
case ('c11')
constitutive_phenopowerlaw_C11(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_C11(i) = IO_floatValue(line,positions,2_pInt)
case ('c12')
constitutive_phenopowerlaw_C12(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_C12(i) = IO_floatValue(line,positions,2_pInt)
case ('c13')
constitutive_phenopowerlaw_C13(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_C13(i) = IO_floatValue(line,positions,2_pInt)
case ('c33')
constitutive_phenopowerlaw_C33(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_C33(i) = IO_floatValue(line,positions,2_pInt)
case ('c44')
constitutive_phenopowerlaw_C44(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_C44(i) = IO_floatValue(line,positions,2_pInt)
case ('nslip')
forall (j = 1:lattice_maxNslipFamily) constitutive_phenopowerlaw_Nslip(j,i) = IO_intValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily)&
constitutive_phenopowerlaw_Nslip(j,i) = IO_intValue(line,positions,1_pInt+j)
case ('gdot0_slip')
constitutive_phenopowerlaw_gdot0_slip(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_gdot0_slip(i) = IO_floatValue(line,positions,2_pInt)
case ('n_slip')
constitutive_phenopowerlaw_n_slip(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_n_slip(i) = IO_floatValue(line,positions,2_pInt)
case ('tau0_slip')
forall (j = 1:lattice_maxNslipFamily) constitutive_phenopowerlaw_tau0_slip(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily)&
constitutive_phenopowerlaw_tau0_slip(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('tausat_slip')
forall (j = 1:lattice_maxNslipFamily) constitutive_phenopowerlaw_tausat_slip(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily)&
constitutive_phenopowerlaw_tausat_slip(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('a_slip', 'w0_slip')
constitutive_phenopowerlaw_a_slip(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_a_slip(i) = IO_floatValue(line,positions,2_pInt)
case ('ntwin')
forall (j = 1:lattice_maxNtwinFamily) constitutive_phenopowerlaw_Ntwin(j,i) = IO_intValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily)&
constitutive_phenopowerlaw_Ntwin(j,i) = IO_intValue(line,positions,1_pInt+j)
case ('gdot0_twin')
constitutive_phenopowerlaw_gdot0_twin(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_gdot0_twin(i) = IO_floatValue(line,positions,2_pInt)
case ('n_twin')
constitutive_phenopowerlaw_n_twin(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_n_twin(i) = IO_floatValue(line,positions,2_pInt)
case ('tau0_twin')
forall (j = 1:lattice_maxNtwinFamily) constitutive_phenopowerlaw_tau0_twin(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily)&
constitutive_phenopowerlaw_tau0_twin(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('s_pr')
constitutive_phenopowerlaw_spr(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_spr(i) = IO_floatValue(line,positions,2_pInt)
case ('twin_b')
constitutive_phenopowerlaw_twinB(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_twinB(i) = IO_floatValue(line,positions,2_pInt)
case ('twin_c')
constitutive_phenopowerlaw_twinC(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_twinC(i) = IO_floatValue(line,positions,2_pInt)
case ('twin_d')
constitutive_phenopowerlaw_twinD(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_twinD(i) = IO_floatValue(line,positions,2_pInt)
case ('twin_e')
constitutive_phenopowerlaw_twinE(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_twinE(i) = IO_floatValue(line,positions,2_pInt)
case ('h0_slipslip')
constitutive_phenopowerlaw_h0_slipslip(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_h0_slipslip(i) = IO_floatValue(line,positions,2_pInt)
case ('h0_sliptwin')
constitutive_phenopowerlaw_h0_sliptwin(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_h0_sliptwin(i) = IO_floatValue(line,positions,2_pInt)
case ('h0_twinslip')
constitutive_phenopowerlaw_h0_twinslip(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_h0_twinslip(i) = IO_floatValue(line,positions,2_pInt)
case ('h0_twintwin')
constitutive_phenopowerlaw_h0_twintwin(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_h0_twintwin(i) = IO_floatValue(line,positions,2_pInt)
case ('atol_resistance')
constitutive_phenopowerlaw_aTolResistance(i) = IO_floatValue(line,positions,2)
constitutive_phenopowerlaw_aTolResistance(i) = IO_floatValue(line,positions,2_pInt)
case ('interaction_slipslip')
forall (j = 1:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_slipslip(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_slipslip(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('interaction_sliptwin')
forall (j = 1:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_sliptwin(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_sliptwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('interaction_twinslip')
forall (j = 1:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_twinslip(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_twinslip(j,i) = IO_floatValue(line,positions,1_pInt+j)
case ('interaction_twintwin')
forall (j = 1:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_twintwin(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_twintwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
end select
endif
enddo
100 do i = 1,maxNinstance
100 do i = 1_pInt,maxNinstance
constitutive_phenopowerlaw_structure(i) = lattice_initializeStructure(constitutive_phenopowerlaw_structureName(i), & ! get structure
constitutive_phenopowerlaw_CoverA(i))
@ -357,21 +362,21 @@ subroutine constitutive_phenopowerlaw_init(file)
constitutive_phenopowerlaw_totalNslip(i) = sum(constitutive_phenopowerlaw_Nslip(:,i)) ! how many slip systems altogether
constitutive_phenopowerlaw_totalNtwin(i) = sum(constitutive_phenopowerlaw_Ntwin(:,i)) ! how many twin systems altogether
if (constitutive_phenopowerlaw_structure(i) < 1 ) call IO_error(205,i)
if (constitutive_phenopowerlaw_structure(i) < 1 ) call IO_error(205_pInt,i)
if (any(constitutive_phenopowerlaw_tau0_slip(:,i) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(210,i)
if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(211,i)
if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(212,i)
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(210_pInt,i)
if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(211_pInt,i)
if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(212_pInt,i)
if (any(constitutive_phenopowerlaw_tausat_slip(:,i) <= 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(213,i)
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(213_pInt,i)
if (any(constitutive_phenopowerlaw_a_slip(i) == 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(214,i)
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(214_pInt,i)
if (any(constitutive_phenopowerlaw_tau0_twin(:,i) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(210,i)
constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(210_pInt,i)
if ( constitutive_phenopowerlaw_gdot0_twin(i) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211,i)
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211_pInt,i)
if ( constitutive_phenopowerlaw_n_twin(i) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(212,i)
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(212_pInt,i)
if (constitutive_phenopowerlaw_aTolResistance(i) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolResistance(i) = 1.0_pReal ! default absolute tolerance 1 Pa
@ -395,8 +400,8 @@ subroutine constitutive_phenopowerlaw_init(file)
constitutive_phenopowerlaw_hardeningMatrix_twintwin = 0.0_pReal
do i = 1,maxNinstance
do j = 1,maxval(phase_Noutput)
do i = 1_pInt,maxNinstance
do j = 1_pInt,maxval(phase_Noutput)
select case(constitutive_phenopowerlaw_output(j,i))
case('resistance_slip', &
'shearrate_slip', &
@ -424,21 +429,21 @@ subroutine constitutive_phenopowerlaw_init(file)
enddo
constitutive_phenopowerlaw_sizeDotState(i) = constitutive_phenopowerlaw_totalNslip(i)+ &
constitutive_phenopowerlaw_totalNtwin(i)+ 2 ! s_slip, s_twin, sum(gamma), sum(f)
constitutive_phenopowerlaw_totalNtwin(i)+ 2_pInt ! s_slip, s_twin, sum(gamma), sum(f)
constitutive_phenopowerlaw_sizeState(i) = constitutive_phenopowerlaw_totalNslip(i)+ &
constitutive_phenopowerlaw_totalNtwin(i)+ 2 ! s_slip, s_twin, sum(gamma), sum(f)
constitutive_phenopowerlaw_totalNtwin(i)+ 2_pInt ! s_slip, s_twin, sum(gamma), sum(f)
myStructure = constitutive_phenopowerlaw_structure(i)
select case (lattice_symmetryType(myStructure)) ! assign elasticity tensor
case(1) ! cubic(s)
forall(k=1:3)
forall(j=1:3) &
case(1_pInt) ! cubic(s)
forall(k=1_pInt:3_pInt)
forall(j=1_pInt:3_pInt) &
constitutive_phenopowerlaw_Cslip_66(k,j,i) = constitutive_phenopowerlaw_C12(i)
constitutive_phenopowerlaw_Cslip_66(k,k,i) = constitutive_phenopowerlaw_C11(i)
constitutive_phenopowerlaw_Cslip_66(k+3,k+3,i) = constitutive_phenopowerlaw_C44(i)
constitutive_phenopowerlaw_Cslip_66(k+3_pInt,k+3_pInt,i) = constitutive_phenopowerlaw_C44(i)
end forall
case(2) ! hex
case(2_pInt) ! hex
constitutive_phenopowerlaw_Cslip_66(1,1,i) = constitutive_phenopowerlaw_C11(i)
constitutive_phenopowerlaw_Cslip_66(2,2,i) = constitutive_phenopowerlaw_C11(i)
constitutive_phenopowerlaw_Cslip_66(3,3,i) = constitutive_phenopowerlaw_C33(i)
@ -456,12 +461,12 @@ subroutine constitutive_phenopowerlaw_init(file)
constitutive_phenopowerlaw_Cslip_66(:,:,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_phenopowerlaw_Cslip_66(:,:,i)))
do f = 1,lattice_maxNslipFamily ! >>> interaction slip -- X
index_myFamily = sum(constitutive_phenopowerlaw_Nslip(1:f-1,i))
do j = 1,constitutive_phenopowerlaw_Nslip(f,i) ! loop over (active) systems in my family (slip)
do o = 1,lattice_maxNslipFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Nslip(1:o-1,i))
do k = 1,constitutive_phenopowerlaw_Nslip(o,i) ! loop over (active) systems in other family (slip)
do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X
index_myFamily = sum(constitutive_phenopowerlaw_Nslip(1:f-1_pInt,i))
do j = 1_pInt,constitutive_phenopowerlaw_Nslip(f,i) ! loop over (active) systems in my family (slip)
do o = 1_pInt,lattice_maxNslipFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Nslip(1:o-1_pInt,i))
do k = 1_pInt,constitutive_phenopowerlaw_Nslip(o,i) ! loop over (active) systems in other family (slip)
constitutive_phenopowerlaw_hardeningMatrix_slipslip(index_otherFamily+k,index_myFamily+j,i) = &
constitutive_phenopowerlaw_interaction_slipslip(lattice_interactionSlipSlip( &
sum(lattice_NslipSystem(1:o-1,myStructure))+k, &
@ -469,39 +474,39 @@ subroutine constitutive_phenopowerlaw_init(file)
myStructure), i )
enddo; enddo
do o = 1,lattice_maxNtwinFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Ntwin(1:o-1,i))
do k = 1,constitutive_phenopowerlaw_Ntwin(o,i) ! loop over (active) systems in other family (twin)
do o = 1_pInt,lattice_maxNtwinFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Ntwin(1:o-1_pInt,i))
do k = 1_pInt,constitutive_phenopowerlaw_Ntwin(o,i) ! loop over (active) systems in other family (twin)
constitutive_phenopowerlaw_hardeningMatrix_sliptwin(index_otherFamily+k,index_myFamily+j,i) = &
constitutive_phenopowerlaw_interaction_sliptwin(lattice_interactionSlipTwin( &
sum(lattice_NtwinSystem(1:o-1,myStructure))+k, &
sum(lattice_NslipSystem(1:f-1,myStructure))+j, &
sum(lattice_NtwinSystem(1:o-1_pInt,myStructure))+k, &
sum(lattice_NslipSystem(1:f-1_pInt,myStructure))+j, &
myStructure), i )
enddo; enddo
enddo; enddo
do f = 1,lattice_maxNtwinFamily ! >>> interaction twin -- X
index_myFamily = sum(constitutive_phenopowerlaw_Ntwin(1:f-1,i))
do j = 1,constitutive_phenopowerlaw_Ntwin(f,i) ! loop over (active) systems in my family (twin)
do f = 1_pInt,lattice_maxNtwinFamily ! >>> interaction twin -- X
index_myFamily = sum(constitutive_phenopowerlaw_Ntwin(1:f-1_pInt,i))
do j = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,i) ! loop over (active) systems in my family (twin)
do o = 1,lattice_maxNslipFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Nslip(1:o-1,i))
do k = 1,constitutive_phenopowerlaw_Nslip(o,i) ! loop over (active) systems in other family (slip)
do o = 1_pInt,lattice_maxNslipFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Nslip(1:o-1_pInt,i))
do k = 1_pInt,constitutive_phenopowerlaw_Nslip(o,i) ! loop over (active) systems in other family (slip)
constitutive_phenopowerlaw_hardeningMatrix_twinslip(index_otherFamily+k,index_myFamily+j,i) = &
constitutive_phenopowerlaw_interaction_twinslip(lattice_interactionTwinSlip( &
sum(lattice_NslipSystem(1:o-1,myStructure))+k, &
sum(lattice_NtwinSystem(1:f-1,myStructure))+j, &
sum(lattice_NslipSystem(1:o-1_pInt,myStructure))+k, &
sum(lattice_NtwinSystem(1:f-1_pInt,myStructure))+j, &
myStructure), i )
enddo; enddo
do o = 1,lattice_maxNtwinFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Ntwin(1:o-1,i))
do k = 1,constitutive_phenopowerlaw_Ntwin(o,i) ! loop over (active) systems in other family (twin)
do o = 1_pInt,lattice_maxNtwinFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Ntwin(1:o-1_pInt,i))
do k = 1_pInt,constitutive_phenopowerlaw_Ntwin(o,i) ! loop over (active) systems in other family (twin)
constitutive_phenopowerlaw_hardeningMatrix_twintwin(index_otherFamily+k,index_myFamily+j,i) = &
constitutive_phenopowerlaw_interaction_twintwin(lattice_interactionTwinTwin( &
sum(lattice_NtwinSystem(1:o-1,myStructure))+k, &
sum(lattice_NtwinSystem(1:f-1,myStructure))+j, &
sum(lattice_NtwinSystem(1:o-1_pInt,myStructure))+k, &
sum(lattice_NtwinSystem(1:f-1_pInt,myStructure))+j, &
myStructure), i )
enddo; enddo

259
code/constitutive_titanmod.f90
View File

@ -389,199 +389,199 @@ do ! read thru sections of
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
section = section + 1_pInt
output = 0_pInt ! reset output counter
endif
if (section > 0 .and. phase_constitution(section) == constitutive_titanmod_label) then ! one of my sections
if (section > 0_pInt .and. phase_constitution(section) == constitutive_titanmod_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 ('(output)')
output = output + 1
constitutive_titanmod_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
output = output + 1_pInt
constitutive_titanmod_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
write(6,*) tag,constitutive_titanmod_output(output,i)
case ('lattice_structure')
constitutive_titanmod_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
constitutive_titanmod_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt))
write(6,*) tag,constitutive_titanmod_structureName(i)
case ('covera_ratio')
constitutive_titanmod_CoverA(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_CoverA(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag,constitutive_titanmod_CoverA(i)
case ('c11')
constitutive_titanmod_C11(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_C11(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag,constitutive_titanmod_C11(i)
case ('c12')
constitutive_titanmod_C12(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_C12(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag,constitutive_titanmod_C12(i)
case ('c13')
constitutive_titanmod_C13(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_C13(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag,constitutive_titanmod_C13(i)
case ('c33')
constitutive_titanmod_C33(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_C33(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag,constitutive_titanmod_C33(i)
case ('c44')
constitutive_titanmod_C44(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_C44(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag,constitutive_titanmod_C44(i)
case ('debyefrequency')
constitutive_titanmod_debyefrequency(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_debyefrequency(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag,constitutive_titanmod_debyefrequency(i)
case ('kinkf0')
constitutive_titanmod_kinkf0(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_kinkf0(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag,constitutive_titanmod_kinkf0(i)
case ('nslip')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_Nslip(j,i) = IO_intValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_Nslip(j,i) = IO_intValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_Nslip(1:4,i)
case ('ntwin')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_Ntwin(j,i) = IO_intValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_Ntwin(j,i) = IO_intValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_Ntwin(1:4,i)
case ('rho_edge0')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_rho_edge0(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
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) &
constitutive_titanmod_rho_screw0(j,i) = IO_floatValue(line,positions,1+j)
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')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_burgersPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_burgersPerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_burgersPerSlipFamily(1:4,i)
case ('twinburgers')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_burgersPerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_burgersPerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_burgersPerTwinFamily(1:4,i)
case ('f0')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_f0_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_f0_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_f0_PerSlipFamily(1:4,i)
case ('twinf0')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_twinf0_PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_twinf0_PerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_twinf0_PerTwinFamily(1:4,i)
case ('tau0e')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_tau0e_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_tau0e_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_tau0e_PerSlipFamily(1:4,i)
case ('twintau0')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_twintau0_PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_twintau0_PerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_twintau0_PerTwinFamily(1:4,i)
case ('tau0s')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_tau0s_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_tau0s_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_tau0s_PerSlipFamily(1:4,i)
case ('capre')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_capre_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_capre_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_capre_PerSlipFamily(1:4,i)
case ('caprs')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_caprs_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_caprs_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_caprs_PerSlipFamily(1:4,i)
case ('v0e')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_v0e_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_v0e_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_v0e_PerSlipFamily(1:4,i)
case ('twingamma0')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_twingamma0_PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_twingamma0_PerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_twingamma0_PerTwinFamily(1:4,i)
case ('v0s')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_v0s_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_v0s_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_v0s_PerSlipFamily(1:4,i)
case ('kinkcriticallength')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_kinkcriticallength_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_kinkcriticallength_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_kinkcriticallength_PerSlipFamily(1:4,i)
case ('twinsize')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_twinsizePerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_twinsizePerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case ('celambdaslip')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_CeLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_CeLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case ('twinlambdaslip')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_twinlambdaslipPerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_twinlambdaslipPerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_twinlambdaslipPerTwinFamily(1:4,i)
case ('cslambdaslip')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_CsLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_CsLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case ('grainsize')
constitutive_titanmod_GrainSize(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_GrainSize(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag
case ('maxtwinfraction')
constitutive_titanmod_MaxTwinFraction(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_MaxTwinFraction(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag
case ('pe')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_pe_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_pe_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_pe_PerSlipFamily(1:4,i)
case ('twinp')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_twinp_PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_twinp_PerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_twinp_PerTwinFamily(1:4,i)
case ('ps')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_ps_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_ps_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_ps_PerSlipFamily(1:4,i)
case ('qe')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_qe_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_qe_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_qe_PerSlipFamily(1:4,i)
case ('twinq')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_twinq_PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_twinq_PerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_twinq_PerTwinFamily(1:4,i)
case ('qs')
forall (j = 1:lattice_maxNslipFamily) &
constitutive_titanmod_qs_PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNslipFamily) &
constitutive_titanmod_qs_PerSlipFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_qs_PerSlipFamily(1:4,i)
case ('twinshearconstant')
forall (j = 1:lattice_maxNtwinFamily) &
constitutive_titanmod_twinshearconstant_PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNtwinFamily) &
constitutive_titanmod_twinshearconstant_PerTwinFamily(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag,constitutive_titanmod_twinshearconstant_PerTwinFamily(1:4,i)
case ('dc')
constitutive_titanmod_dc(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_dc(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag
case ('twinhpconstant')
constitutive_titanmod_twinhpconstant(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_twinhpconstant(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag
case ('atol_rho')
constitutive_titanmod_aTolRho(i) = IO_floatValue(line,positions,2)
constitutive_titanmod_aTolRho(i) = IO_floatValue(line,positions,2_pInt)
write(6,*) tag
case ('interactionslipslip')
forall (j = 1:lattice_maxNinteraction) &
constitutive_titanmod_interactionSlipSlip(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_titanmod_interactionSlipSlip(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case ('interactionee')
forall (j = 1:lattice_maxNinteraction) &
constitutive_titanmod_interaction_ee(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_titanmod_interaction_ee(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case ('interactionss')
forall (j = 1:lattice_maxNinteraction) &
constitutive_titanmod_interaction_ss(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_titanmod_interaction_ss(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case ('interactiones')
forall (j = 1:lattice_maxNinteraction) &
constitutive_titanmod_interaction_es(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_titanmod_interaction_es(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case ('interactionsliptwin')
forall (j = 1:lattice_maxNinteraction) &
constitutive_titanmod_interactionSlipTwin(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_titanmod_interactionSlipTwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case ('interactiontwinslip')
forall (j = 1:lattice_maxNinteraction) &
constitutive_titanmod_interactionTwinSlip(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_titanmod_interactionTwinSlip(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
case ('interactiontwintwin')
forall (j = 1:lattice_maxNinteraction) &
constitutive_titanmod_interactionTwinTwin(j,i) = IO_floatValue(line,positions,1+j)
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_titanmod_interactionTwinTwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
write(6,*) tag
end select
endif
@ -589,43 +589,43 @@ enddo
write(6,*) 'Material Property reading done'
100 do i = 1,maxNinstance
100 do i = 1_pInt,maxNinstance
constitutive_titanmod_structure(i) = &
lattice_initializeStructure(constitutive_titanmod_structureName(i),constitutive_titanmod_CoverA(i))
myStructure = constitutive_titanmod_structure(i)
!* Sanity checks
if (myStructure < 1 .or. myStructure > 3) call IO_error(205)
if (sum(constitutive_titanmod_Nslip(:,i)) <= 0_pInt) call IO_error(207)
if (sum(constitutive_titanmod_Ntwin(:,i)) < 0_pInt) call IO_error(208) !***
if (myStructure < 1_pInt .or. myStructure > 3_pInt) call IO_error(205_pInt)
if (sum(constitutive_titanmod_Nslip(:,i)) <= 0_pInt) call IO_error(207_pInt)
if (sum(constitutive_titanmod_Ntwin(:,i)) < 0_pInt) call IO_error(208_pInt) !***
do f = 1,lattice_maxNslipFamily
if (constitutive_titanmod_Nslip(f,i) > 0_pInt) then
if (constitutive_titanmod_rho_edge0(f,i) < 0.0_pReal) call IO_error(209)
if (constitutive_titanmod_rho_screw0(f,i) < 0.0_pReal) call IO_error(210)
if (constitutive_titanmod_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(211)
if (constitutive_titanmod_f0_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(212)
if (constitutive_titanmod_tau0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(229)
if (constitutive_titanmod_tau0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(233)
if (constitutive_titanmod_capre_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(234)
if (constitutive_titanmod_caprs_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235)
if (constitutive_titanmod_v0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226)
if (constitutive_titanmod_v0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226)
if (constitutive_titanmod_kinkcriticallength_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(238)
if (constitutive_titanmod_rho_edge0(f,i) < 0.0_pReal) call IO_error(209_pInt)
if (constitutive_titanmod_rho_screw0(f,i) < 0.0_pReal) call IO_error(210_pInt)
if (constitutive_titanmod_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(211_pInt)
if (constitutive_titanmod_f0_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(212_pInt)
if (constitutive_titanmod_tau0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(229_pInt)
if (constitutive_titanmod_tau0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(233_pInt)
if (constitutive_titanmod_capre_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(234_pInt)
if (constitutive_titanmod_caprs_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235_pInt)
if (constitutive_titanmod_v0e_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt)
if (constitutive_titanmod_v0s_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt)
if (constitutive_titanmod_kinkcriticallength_PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(238_pInt)
endif
enddo
do f = 1,lattice_maxNtwinFamily
if (constitutive_titanmod_Ntwin(f,i) > 0_pInt) then
if (constitutive_titanmod_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221) !***
if (constitutive_titanmod_twinf0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(228)
if (constitutive_titanmod_twinshearconstant_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(228)
if (constitutive_titanmod_twintau0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(229)
if (constitutive_titanmod_twingamma0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(226)
if (constitutive_titanmod_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221_pInt) !***
if (constitutive_titanmod_twinf0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(228_pInt)
if (constitutive_titanmod_twinshearconstant_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(228_pInt)
if (constitutive_titanmod_twintau0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(229_pInt)
if (constitutive_titanmod_twingamma0_PerTwinFamily(f,i) <= 0.0_pReal) call IO_error(226_pInt)
endif
enddo
! if (any(constitutive_titanmod_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) call IO_error(229)
if (constitutive_titanmod_dc(i) <= 0.0_pReal) call IO_error(231)
if (constitutive_titanmod_twinhpconstant(i) <= 0.0_pReal) call IO_error(232)
if (constitutive_titanmod_aTolRho(i) <= 0.0_pReal) call IO_error(233)
if (constitutive_titanmod_dc(i) <= 0.0_pReal) call IO_error(231_pInt)
if (constitutive_titanmod_twinhpconstant(i) <= 0.0_pReal) call IO_error(232_pInt)
if (constitutive_titanmod_aTolRho(i) <= 0.0_pReal) call IO_error(233_pInt)
!* Determine total number of active slip or twin systems
constitutive_titanmod_Nslip(:,i) = min(lattice_NslipSystem(:,myStructure),constitutive_titanmod_Nslip(:,i))
@ -724,25 +724,25 @@ constitutive_titanmod_Ctwin_66 = 0.0_pReal
constitutive_titanmod_Ctwin_3333 = 0.0_pReal
write(6,*) 'Allocated slip system variables'
do i = 1,maxNinstance
do i = 1_pInt,maxNinstance
myStructure = constitutive_titanmod_structure(i)
!* Inverse lookup of slip system family
l = 0_pInt
do f = 1,lattice_maxNslipFamily
do k = 1,constitutive_titanmod_Nslip(f,i)
l = l + 1
do f = 1_pInt,lattice_maxNslipFamily
do k = 1_pInt,constitutive_titanmod_Nslip(f,i)
l = l + 1_pInt
constitutive_titanmod_slipFamily(l,i) = f
constitutive_titanmod_slipSystemLattice(l,i) = sum(lattice_NslipSystem(1:f-1,myStructure)) + k
constitutive_titanmod_slipSystemLattice(l,i) = sum(lattice_NslipSystem(1:f-1_pInt,myStructure)) + k
enddo; enddo
!* Inverse lookup of twin system family
l = 0_pInt
do f = 1,lattice_maxNtwinFamily
do k = 1,constitutive_titanmod_Ntwin(f,i)
l = l + 1
do f = 1_pInt,lattice_maxNtwinFamily
do k = 1_pInt,constitutive_titanmod_Ntwin(f,i)
l = l + 1_pInt
constitutive_titanmod_twinFamily(l,i) = f
constitutive_titanmod_twinSystemLattice(l,i) = sum(lattice_NtwinSystem(1:f-1,myStructure)) + k
constitutive_titanmod_twinSystemLattice(l,i) = sum(lattice_NtwinSystem(1:f-1_pInt,myStructure)) + k
enddo; enddo
!* Determine size of state array
@ -809,14 +809,14 @@ write(6,*) 'Determining elasticity matrix'
!* Elasticity matrix and shear modulus according to material.config
select case (myStructure)
case(1:2) ! cubic(s)
forall(k=1:3)
forall(j=1:3) &
case(1_pInt:2_pInt) ! cubic(s)
forall(k=1_pInt:3_pInt)
forall(j=1_pInt:3_pInt) &
constitutive_titanmod_Cslip_66(k,j,i) = constitutive_titanmod_C12(i)
constitutive_titanmod_Cslip_66(k,k,i) = constitutive_titanmod_C11(i)
constitutive_titanmod_Cslip_66(k+3,k+3,i) = constitutive_titanmod_C44(i)
constitutive_titanmod_Cslip_66(k+3_pInt,k+3_pInt,i) = constitutive_titanmod_C44(i)
end forall
case(3:) ! all hex
case(3_pInt:) ! all hex
constitutive_titanmod_Cslip_66(1,1,i) = constitutive_titanmod_C11(i)
constitutive_titanmod_Cslip_66(2,2,i) = constitutive_titanmod_C11(i)
constitutive_titanmod_Cslip_66(3,3,i) = constitutive_titanmod_C33(i)
@ -836,9 +836,10 @@ write(6,*) 'Determining elasticity matrix'
0.2_pReal*(constitutive_titanmod_C11(i)-constitutive_titanmod_C12(i))+0.3_pReal*constitutive_titanmod_C44(i)
!* Construction of the twin elasticity matrices
do j=1,lattice_maxNtwinFamily
do k=1,constitutive_titanmod_Ntwin(j,i)
do l=1,3 ; do m=1,3 ; do n=1,3 ; do o=1,3 ; do p=1,3 ; do q=1,3 ; do r=1,3 ; do s=1,3
do j=1_pInt,lattice_maxNtwinFamily
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_Cslip_3333(p,q,r,s,i)*&

116
code/crystallite.f90
View File

@ -103,6 +103,7 @@ 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 prec, only: pInt, &
pReal
use debug, only: debug_info, &
@ -221,8 +222,7 @@ allocate(crystallite_orientation(4,gMax,iMax,eMax)); crystalli
allocate(crystallite_orientation0(4,gMax,iMax,eMax)); crystallite_orientation0 = 0.0_pReal
allocate(crystallite_rotation(4,gMax,iMax,eMax)); crystallite_rotation = 0.0_pReal
allocate(crystallite_disorientation(4,nMax,gMax,iMax,eMax)); crystallite_disorientation = 0.0_pReal
allocate(crystallite_symmetryID(gMax,iMax,eMax)); crystallite_symmetryID = 0_pInt
allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true.
allocate(crystallite_symmetryID(gMax,iMax,eMax)); crystallite_symmetryID = 0_pIntallocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true.
allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false.
allocate(crystallite_todo(gMax,iMax,eMax)); crystallite_todo = .false.
allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true.
@ -234,10 +234,10 @@ allocate(crystallite_sizePostResult(maxval(crystallite_Noutput), &
if (.not. IO_open_jobFile(file,material_localFileExt)) then ! no local material configuration present...
if (.not. IO_open_file(file,material_configFile)) call IO_error(100) ! ...and cannot open material.config file
if (.not. IO_open_file(file,material_configFile)) call IO_error(100_pInt) ! ...and cannot open material.config file
endif
line = ''
section = 0
section = 0_pInt
do while (IO_lc(IO_getTag(line,'<','>')) /= material_partCrystallite) ! wind forward to <crystallite>
read(file,'(a1024)',END=100) line
@ -248,38 +248,38 @@ do ! read thru sections of
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
section = section + 1_pInt
output = 0_pInt ! reset output counter
endif
if (section > 0) then
if (section > 0_pInt) then
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 ('(output)')
output = output + 1
crystallite_output(output,section) = IO_lc(IO_stringValue(line,positions,2))
output = output + 1_pInt
crystallite_output(output,section) = IO_lc(IO_stringValue(line,positions,2_pInt))
end select
endif
enddo
100 close(file)
do i = 1,material_Ncrystallite ! sanity checks
do i = 1_pInt,material_Ncrystallite ! sanity checks
enddo
do i = 1,material_Ncrystallite
do j = 1,crystallite_Noutput(i)
do i = 1_pInt,material_Ncrystallite
do j = 1_pInt,crystallite_Noutput(i)
select case(crystallite_output(j,i))
case('phase','texture','volume')
mySize = 1
mySize = 1_pInt
case('orientation','grainrotation') ! orientation as quaternion, or deviation from initial grain orientation in axis-angle form (angle in degrees)
mySize = 4
mySize = 4_pInt
case('eulerangles') ! Bunge (3-1-3) Euler angles
mySize = 3
mySize = 3_pInt
case('defgrad','f','fe','fp','lp','e','ee','p','firstpiola','1stpiola','s','tstar','secondpiola','2ndpiola')
mySize = 9
mySize = 9_pInt
case default
mySize = 0
mySize = 0_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
@ -290,7 +290,7 @@ do i = 1,material_Ncrystallite
enddo
crystallite_maxSizePostResults = 0_pInt
do j = 1,material_Nmicrostructure
do j = 1_pInt,material_Nmicrostructure
if (microstructure_active(j)) &
crystallite_maxSizePostResults = max(crystallite_maxSizePostResults,&
crystallite_sizePostResults(microstructure_crystallite(j)))
@ -298,9 +298,9 @@ enddo
! write description file for crystallite output
if(.not. IO_write_jobFile(file,'outputCrystallite')) call IO_error (50) ! problems in writing file
if(.not. IO_write_jobFile(file,'outputCrystallite')) call IO_error (50_pInt) ! problems in writing file
do p = 1,material_Ncrystallite
do p = 1_pInt,material_Ncrystallite
write(file,*)
write(file,'(a)') '['//trim(crystallite_name(p))//']'
write(file,*)
@ -619,10 +619,10 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
crystallite_subTstar0_v(1:6,g,i,e) = crystallite_Tstar_v(1:6,g,i,e) ! ...2nd PK stress
!$OMP FLUSH(crystallite_subF0)
elseif (formerSubStep > subStepMinCryst) then ! this crystallite just converged
if (debug_verbosity > 4) then
if (debug_verbosity > 4_pInt) then
!$OMP CRITICAL (distributionCrystallite)
debug_CrystalliteLoopDistribution(min(nCryst+1,NiterationCrystallite)) = &
debug_CrystalliteLoopDistribution(min(nCryst+1,NiterationCrystallite)) + 1
debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) = &
debug_CrystalliteLoopDistribution(min(nCryst+1_pInt,NiterationCrystallite)) + 1_pInt
!$OMP END CRITICAL (distributionCrystallite)
endif
endif
@ -640,7 +640,7 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
! cant restore dotState here, since not yet calculated in first cutback after initialization
!$OMP FLUSH(crystallite_invFp)
#ifndef _OPENMP
if (debug_verbosity > 4 &
if (debug_verbosity > 4_pInt &
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,f10.8)') '<< CRYST >> cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
crystallite_subStep(g,i,e)
@ -672,20 +672,20 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
if (any(crystallite_todo)) then
select case(numerics_integrator(numerics_integrationMode))
case(1)
case(1_pInt)
call crystallite_integrateStateFPI()
case(2)
case(2_pInt)
call crystallite_integrateStateEuler()
case(3)
case(3_pInt)
call crystallite_integrateStateAdaptiveEuler()
case(4)
case(4_pInt)
call crystallite_integrateStateRK4()
case(5)
case(5_pInt)
call crystallite_integrateStateRKCK45()
endselect
endif
NiterationCrystallite = NiterationCrystallite + 1
NiterationCrystallite = NiterationCrystallite + 1_pInt
enddo ! cutback loop
@ -774,12 +774,12 @@ if(updateJaco) then
! --- INITIALIZE UNPERTURBED STATE ---
select case(numerics_integrator(numerics_integrationMode))
case(1) ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state
case(1_pInt) ! Fix-point method: restore to last converged state at end of subinc, since this is probably closest to perturbed state
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState)
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
mySizeState = constitutive_sizeState(g,i,e)
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState)
@ -792,13 +792,13 @@ if(updateJaco) then
crystallite_Fe = Fe_backup
crystallite_Lp = Lp_backup
crystallite_Tstar_v = Tstar_v_backup
case(2,3) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step
case(4,5) ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress
case(2_pInt,3_pInt) ! explicit Euler methods: nothing to restore (except for F), since we are only doing a stress integration step
case(4_pInt,5_pInt) ! explicit Runge-Kutta methods: restore to start of subinc, since we are doing a full integration of state and stress
!$OMP PARALLEL DO PRIVATE(myNgrains,mySizeState,mySizeDotState)
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
mySizeState = constitutive_sizeState(g,i,e)
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_subState0(g,i,e)%p(1:mySizeState)
@ -822,15 +822,15 @@ if(updateJaco) then
where (crystallite_todo) crystallite_converged = .false. ! start out non-converged
select case(numerics_integrator(numerics_integrationMode))
case(1)
case(1_pInt)
call crystallite_integrateStateFPI()
case(2)
case(2_pInt)
call crystallite_integrateStateEuler()
case(3)
case(3_pInt)
call crystallite_integrateStateAdaptiveEuler()
case(4)
case(4_pInt)
call crystallite_integrateStateRK4()
case(5)
case(5_pInt)
call crystallite_integrateStateRKCK45()
end select
@ -838,12 +838,12 @@ if(updateJaco) then
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
if (crystallite_requested(g,i,e) .and. crystallite_converged(g,i,e)) then ! converged state warrants stiffness update
select case(perturbation)
case(1)
case(1_pInt)
dPdF_perturbation1(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl
case(2)
case(2_pInt)
dPdF_perturbation2(1:3,1:3,k,l,g,i,e) = (crystallite_P(1:3,1:3,g,i,e) - P_backup(1:3,1:3,g,i,e)) / myPert ! tangent dP_ij/dFg_kl
end select
endif
@ -861,14 +861,14 @@ if(updateJaco) then
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
if (crystallite_requested(g,i,e) .and. convergenceFlag_backup(g,i,e)) then ! central solution converged
select case(pert_method)
case(1)
case(1_pInt)
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e)
case(2)
case(2_pInt)
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e)
case(3)
case(3_pInt)
crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.5_pReal* ( dPdF_perturbation1(1:3,1:3,1:3,1:3,g,i,e) &
+ dPdF_perturbation2(1:3,1:3,1:3,1:3,g,i,e))
end select
@ -887,7 +887,7 @@ if(updateJaco) then
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
mySizeState = constitutive_sizeState(g,i,e)
mySizeDotState = constitutive_sizeDotState(g,i,e)
constitutive_state(g,i,e)%p(1:mySizeState) = constitutive_state_backup(g,i,e)%p(1:mySizeState)
@ -977,7 +977,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
@ -1724,7 +1724,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
@ -1996,7 +1996,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
@ -2890,14 +2890,14 @@ LpLoop: do
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
dR_dLp = math_identity2nd(9) - math_mul99x99(dLp_dT_constitutive, dT_dLp)
dR_dLp = math_identity2nd(9_pInt) - math_mul99x99(dLp_dT_constitutive, dT_dLp)
inv_dR_dLp = 0.0_pReal
call math_invert(9,dR_dLp,inv_dR_dLp,dummy,error) ! invert dR/dLp --> dLp/dR
call math_invert(9_pInt,dR_dLp,inv_dR_dLp,dummy,error) ! invert dR/dLp --> dLp/dR
if (error) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
if (debug_verbosity > 4_pInt) then
write(6,'(a,i8,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on dR/dLp inversion at el ip g ',e,i,g
if (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,*)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dR_dLp)
@ -3057,8 +3057,8 @@ logical error
call math_pDecomposition(crystallite_Fe(1:3,1:3,g,i,e), U, R, error) ! polar decomposition of Fe
if (error) then
call IO_warning(650, e, i, g)
orientation = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! fake orientation
call IO_warning(650_pInt, e, i, g)
orientation = [1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] ! fake orientation
else
orientation = math_RtoQuaternion(transpose(R))
endif

140
code/homogenization_RGC.f90
View File

@ -117,44 +117,44 @@ subroutine homogenization_RGC_init(&
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
section = section + 1_pInt
output = 0_pInt ! reset output counter
endif
if (section > 0 .and. homogenization_type(section) == homogenization_RGC_label) then ! one of my sections
if (section > 0_pInt .and. homogenization_type(section) == homogenization_RGC_label) then ! one of my sections
i = homogenization_typeInstance(section) ! which instance of my type is present homogenization
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 ('(output)')
output = output + 1
homogenization_RGC_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
output = output + 1_pInt
homogenization_RGC_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('clustersize')
homogenization_RGC_Ngrains(1,i) = IO_intValue(line,positions,2)
homogenization_RGC_Ngrains(2,i) = IO_intValue(line,positions,3)
homogenization_RGC_Ngrains(3,i) = IO_intValue(line,positions,4)
homogenization_RGC_Ngrains(1,i) = IO_intValue(line,positions,2_pInt)
homogenization_RGC_Ngrains(2,i) = IO_intValue(line,positions,3_pInt)
homogenization_RGC_Ngrains(3,i) = IO_intValue(line,positions,4_pInt)
case ('scalingparameter')
homogenization_RGC_xiAlpha(i) = IO_floatValue(line,positions,2)
homogenization_RGC_xiAlpha(i) = IO_floatValue(line,positions,2_pInt)
case ('overproportionality')
homogenization_RGC_ciAlpha(i) = IO_floatValue(line,positions,2)
homogenization_RGC_ciAlpha(i) = IO_floatValue(line,positions,2_pInt)
case ('grainsize')
homogenization_RGC_dAlpha(1,i) = IO_floatValue(line,positions,2)
homogenization_RGC_dAlpha(2,i) = IO_floatValue(line,positions,3)
homogenization_RGC_dAlpha(3,i) = IO_floatValue(line,positions,4)
homogenization_RGC_dAlpha(1,i) = IO_floatValue(line,positions,2_pInt)
homogenization_RGC_dAlpha(2,i) = IO_floatValue(line,positions,3_pInt)
homogenization_RGC_dAlpha(3,i) = IO_floatValue(line,positions,4_pInt)
case ('clusterorientation')
homogenization_RGC_angles(1,i) = IO_floatValue(line,positions,2)
homogenization_RGC_angles(2,i) = IO_floatValue(line,positions,3)
homogenization_RGC_angles(3,i) = IO_floatValue(line,positions,4)
homogenization_RGC_angles(1,i) = IO_floatValue(line,positions,2_pInt)
homogenization_RGC_angles(2,i) = IO_floatValue(line,positions,3_pInt)
homogenization_RGC_angles(3,i) = IO_floatValue(line,positions,4_pInt)
end select
endif
enddo
!*** assigning cluster orientations
do e = 1,mesh_NcpElems
do e = 1_pInt,mesh_NcpElems
if (homogenization_type(mesh_element(3,e)) == homogenization_RGC_label) then
myInstance = homogenization_typeInstance(mesh_element(3,e))
if (all (homogenization_RGC_angles(:,myInstance) >= 399.9_pReal)) then
homogenization_RGC_orientation(:,:,1,e) = math_EulerToR(math_sampleRandomOri())
do i = 1,FE_Nips(mesh_element(2,e))
do i = 1_pInt,FE_Nips(mesh_element(2,e))
if (microstructure_elemhomo(mesh_element(4,e))) then
homogenization_RGC_orientation(:,:,i,e) = homogenization_RGC_orientation(:,:,1,e)
else
@ -162,44 +162,44 @@ subroutine homogenization_RGC_init(&
endif
enddo
else
do i = 1,FE_Nips(mesh_element(2,e))
do i = 1_pInt,FE_Nips(mesh_element(2,e))
homogenization_RGC_orientation(:,:,i,e) = math_EulerToR(homogenization_RGC_angles(:,myInstance)*inRad)
enddo
endif
endif
enddo
100 if (debug_verbosity == 4) then
100 if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
do i = 1,maxNinstance
do i = 1_pInt,maxNinstance
write(6,'(a15,1x,i4)') 'instance: ', i
write(6,*)
write(6,'(a25,3(1x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1,3)
write(6,'(a25,3(1x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1_pInt,3_pInt)
write(6,'(a25,1x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i)
write(6,'(a25,1x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i)
write(6,'(a25,3(1x,e10.3))') 'grain size: ',(homogenization_RGC_dAlpha(j,i),j=1,3)
write(6,'(a25,3(1x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1,3)
write(6,'(a25,3(1x,e10.3))') 'grain size: ',(homogenization_RGC_dAlpha(j,i),j=1_pInt,3_pInt)
write(6,'(a25,3(1x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1_pInt,3_pInt)
enddo
!$OMP END CRITICAL (write2out)
endif
do i = 1,maxNinstance
do j = 1,maxval(homogenization_Noutput)
do i = 1_pInt,maxNinstance
do j = 1_pInt,maxval(homogenization_Noutput)
select case(homogenization_RGC_output(j,i))
case('constitutivework')
mySize = 1
mySize = 1_pInt
case('magnitudemismatch')
mySize = 3
mySize = 3_pInt
case('penaltyenergy')
mySize = 1
mySize = 1_pInt
case('volumediscrepancy')
mySize = 1
mySize = 1_pInt
case('averagerelaxrate')
mySize = 1
mySize = 1_pInt
case('maximumrelaxrate')
mySize = 1
mySize = 1_pInt
case default
mySize = 0
mySize = 0_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
@ -212,9 +212,9 @@ subroutine homogenization_RGC_init(&
homogenization_RGC_sizeState(i) &
= 3*(homogenization_RGC_Ngrains(1,i)-1)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) &
+ 3*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1)*homogenization_RGC_Ngrains(3,i) &
+ 3*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1) &
= 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) &
+ 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
@ -272,16 +272,16 @@ subroutine homogenization_RGC_partitionDeformation(&
integer(pInt), dimension (3) :: iGrain3
integer(pInt) homID, iGrain,iFace,i,j
!
integer(pInt), parameter :: nFace = 6
integer(pInt), parameter :: nFace = 6_pInt
!* Debugging the overall deformation gradient
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' =========='
write(6,'(1x,a32)')'Overall deformation gradient: '
do i = 1,3
write(6,'(1x,3(e14.8,1x))')(avgF(i,j), j = 1,3)
do i = 1_pInt,3_pInt
write(6,'(1x,3(e14.8,1x))')(avgF(i,j), j = 1_pInt,3_pInt)
enddo
write(6,*)' '
call flush(6)
@ -293,21 +293,21 @@ subroutine homogenization_RGC_partitionDeformation(&
F = 0.0_pReal
do iGrain = 1,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) ! identifying 6 interfaces of each grain
aVect = homogenization_RGC_relaxationVector(intFace,state,homID)! get the relaxation vectors for each interface from global relaxation vector array
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of each interface
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) ! calculating deformation relaxations due to interface relaxation
enddo
F(:,:,iGrain) = F(:,:,iGrain) + avgF(:,:) ! resulting relaxed deformation gradient
!* Debugging the grain deformation gradients
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain
do i = 1,3
write(6,'(1x,3(e14.8,1x))')(F(i,j,iGrain), j = 1,3)
do i = 1_pInt,3_pInt
write(6,'(1x,3(e14.8,1x))')(F(i,j,iGrain), j = 1_pInt,3_pInt)
enddo
write(6,*)' '
call flush(6)
@ -368,7 +368,7 @@ function homogenization_RGC_updateState(&
real(pReal) residMax,stresMax,constitutiveWork,penaltyEnergy,volDiscrep
logical error
!
integer(pInt), parameter :: nFace = 6
integer(pInt), parameter :: nFace = 6_pInt
!
real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
@ -383,11 +383,11 @@ function homogenization_RGC_updateState(&
+ nGDim(1)*nGDim(2)*(nGDim(3)-1)
!* Allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
allocate(resid(3*nIntFaceTot)); resid = 0.0_pReal
allocate(resid(3_pInt*nIntFaceTot)); resid = 0.0_pReal
allocate(tract(nIntFaceTot,3)); tract = 0.0_pReal
allocate(relax(3*nIntFaceTot)); relax = state%p(1:3*nIntFaceTot)
allocate(drelax(3*nIntFaceTot))
drelax = state%p(1:3*nIntFaceTot) - state0%p(1:3*nIntFaceTot)
allocate(relax(3_pInt*nIntFaceTot)); relax = state%p(1:3_pInt*nIntFaceTot)
allocate(drelax(3_pInt*nIntFaceTot))
drelax = state%p(1:3_pInt*nIntFaceTot) - state0%p(1:3_pInt*nIntFaceTot)
!* Debugging the obtained state
if (debug_verbosity == 4) then
@ -432,19 +432,19 @@ function homogenization_RGC_updateState(&
!* 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)
intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N)
normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the 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)
intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P)
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
tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1)))/(refRelaxRate_RGC*dt))**viscPower_RGC, &
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
tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) &
@ -720,15 +720,15 @@ function homogenization_RGC_updateState(&
!* -------------------------------------------------------------------------------------------------------------
!*** Computing the update of the state variable (relaxation vectors) using the Jacobian matrix
allocate(jnverse(3*nIntFaceTot,3*nIntFaceTot)); jnverse = 0.0_pReal
call math_invert(3*nIntFaceTot,jmatrix,jnverse,ival,error) ! Compute the inverse of the overall Jacobian matrix
allocate(jnverse(3_pInt*nIntFaceTot,3_pInt*nIntFaceTot)); jnverse = 0.0_pReal
call math_invert(3_pInt*nIntFaceTot,jmatrix,jnverse,ival,error) ! Compute the inverse of the overall Jacobian matrix
!* Debugging the inverse Jacobian matrix
if (debug_verbosity == 4) then
if (debug_verbosity == 4_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian inverse'
do i = 1,3*nIntFaceTot
write(6,'(1x,100(e10.4,1x))')(jnverse(i,j), j = 1,3*nIntFaceTot)
write(6,'(1x,100(e10.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
call flush(6)
@ -1292,9 +1292,9 @@ function homogenization_RGC_grain1to3(&
!* Get the grain position
nGDim = homogenization_RGC_Ngrains(:,homID)
homogenization_RGC_grain1to3(3) = 1+(grain1-1)/(nGDim(1)*nGDim(2))
homogenization_RGC_grain1to3(2) = 1+mod((grain1-1)/nGDim(1),nGDim(2))
homogenization_RGC_grain1to3(1) = 1+mod((grain1-1),nGDim(1))
homogenization_RGC_grain1to3(3) = 1_pInt+(grain1-1_pInt)/(nGDim(1)*nGDim(2))
homogenization_RGC_grain1to3(2) = 1_pInt+mod((grain1-1_pInt)/nGDim(1),nGDim(2))
homogenization_RGC_grain1to3(1) = 1_pInt+mod((grain1-1_pInt),nGDim(1))
endfunction
@ -1390,21 +1390,23 @@ function homogenization_RGC_interface1to4(&
!* 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
homogenization_RGC_interface1to4(3) = mod((iFace1D-1),nGDim(2))+1
homogenization_RGC_interface1to4(4) = mod(int(dble(iFace1D-1)/dble(nGDim(2))),nGDim(3))+1
homogenization_RGC_interface1to4(2) = int(dble(iFace1D-1)/dble(nGDim(2))/dble(nGDim(3)))+1
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(dble(iFace1D-nIntFace(1)-1)/dble(nGDim(3))),nGDim(1))+1
homogenization_RGC_interface1to4(3) = int(dble(iFace1D-nIntFace(1)-1)/dble(nGDim(3))/dble(nGDim(1)))+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(dble(iFace1D-nIntFace(2)-nIntFace(1)-1)/dble(nGDim(1))),nGDim(2))+1
homogenization_RGC_interface1to4(4) = int(dble(iFace1D-nIntFace(2)-nIntFace(1)-1)/dble(nGDim(1))/dble(nGDim(2)))+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
endif
endfunction

18
code/homogenization_isostrain.f90
View File

@ -81,7 +81,7 @@ subroutine homogenization_isostrain_init(&
!$OMP END CRITICAL (write2out)
maxNinstance = count(homogenization_type == homogenization_isostrain_label)
if (maxNinstance == 0) return
if (maxNinstance == 0_pInt) return
allocate(homogenization_isostrain_sizeState(maxNinstance)) ; homogenization_isostrain_sizeState = 0_pInt
allocate(homogenization_isostrain_sizePostResults(maxNinstance)); homogenization_isostrain_sizePostResults = 0_pInt
@ -93,7 +93,7 @@ subroutine homogenization_isostrain_init(&
rewind(file)
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
@ -104,24 +104,24 @@ subroutine homogenization_isostrain_init(&
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
section = section + 1_pInt
output = 0_pInt ! reset output counter
endif
if (section > 0 .and. homogenization_type(section) == homogenization_isostrain_label) then ! one of my sections
i = homogenization_typeInstance(section) ! which instance of my type is present homogenization
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 ('(output)')
output = output + 1
homogenization_isostrain_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
output = output + 1_pInt
homogenization_isostrain_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('ngrains')
homogenization_isostrain_Ngrains(i) = IO_intValue(line,positions,2)
homogenization_isostrain_Ngrains(i) = IO_intValue(line,positions,2_pInt)
end select
endif
enddo
100 do i = 1,maxNinstance ! sanity checks
100 do i = 1_pInt,maxNinstance ! sanity checks
enddo
do i = 1,maxNinstance

2
processing/post/postResults.py
View File

@ -93,7 +93,7 @@ class MPIEspectral_result: # mimic py_post result object
self.N_loadcases = self._keyedPackedArray('loadcases',count=1,type='i',default=1)[0]
self._frequencies = self._keyedPackedArray('frequencies',count=self.N_loadcases,type='i',default=1)
self._increments = self._keyedPackedArray('increments',count=self.N_loadcases,type='i')
self._increments[0] -= 1 # delete zero'th entry
#self._increments[0] -= 1 # delete zero'th entry (might be needed for older spectralOuts
self._times = self._keyedPackedArray('times',count=self.N_loadcases,type='d',default=0.0)
self._logscales = self._keyedPackedArray('logscales',count=self.N_loadcases,type='i',default=0)
self.dimension = self._keyedPackedArray('dimension',count=3,type='d')