Merge branch 'development' into release_bugfix

This commit is contained in:
Martin Diehl 2018-05-26 16:37:39 +02:00
commit 466cb57d7c
11 changed files with 89 additions and 16798 deletions

View File

@ -1 +0,0 @@
env/DAMASK.csh

View File

@ -1 +0,0 @@
env/DAMASK.sh

View File

@ -1 +0,0 @@
env/DAMASK.zsh

@ -1 +1 @@
Subproject commit b7d1d309146e017caa5744333c2e4a4532a6fc20 Subproject commit cd02f6c1a481491eb4517651516b8311348b4777

View File

@ -1 +1 @@
v2.0.2 v2.0.2-4-g35ba260

File diff suppressed because one or more lines are too long

File diff suppressed because it is too large Load Diff

File diff suppressed because one or more lines are too long

View File

@ -108,7 +108,7 @@ module numerics
character(len=64), private :: & character(len=64), private :: &
fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag
character(len=64), protected, public :: & character(len=64), protected, public :: &
spectral_solver = 'basicpetsc' , & !< spectral solution method spectral_solver = 'basic', & !< spectral solution method
spectral_derivative = 'continuous' !< spectral spatial derivative method spectral_derivative = 'continuous' !< spectral spatial derivative method
character(len=1024), protected, public :: & character(len=1024), protected, public :: &
petsc_defaultOptions = '-mech_snes_type ngmres & petsc_defaultOptions = '-mech_snes_type ngmres &

View File

@ -53,7 +53,7 @@ module plastic_isotropic
dilatation = .false. dilatation = .false.
end type end type
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable, target, private :: param !< containers of constitutive parameters (len Ninstance)
type, private :: tIsotropicState !< internal state aliases type, private :: tIsotropicState !< internal state aliases
real(pReal), pointer, dimension(:) :: & ! scalars along NipcMyInstance real(pReal), pointer, dimension(:) :: & ! scalars along NipcMyInstance
@ -61,20 +61,10 @@ module plastic_isotropic
accumulatedShear accumulatedShear
end type end type
type, private :: tIsotropicAbsTol !< internal alias for abs tolerance in state
real(pReal), pointer :: & ! scalars
flowstress, &
accumulatedShear
end type
type(tIsotropicState), allocatable, dimension(:), private :: & !< state aliases per instance type(tIsotropicState), allocatable, dimension(:), private :: & !< state aliases per instance
state, & state, &
state0, &
dotState dotState
type(tIsotropicAbsTol), allocatable, dimension(:), private :: & !< state aliases per instance
stateAbsTol
public :: & public :: &
plastic_isotropic_init, & plastic_isotropic_init, &
plastic_isotropic_LpAndItsTangent, & plastic_isotropic_LpAndItsTangent, &
@ -130,6 +120,7 @@ subroutine plastic_isotropic_init(fileUnit)
implicit none implicit none
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
type(tParameters), pointer :: p
integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt) :: & integer(pInt) :: &
@ -184,14 +175,11 @@ subroutine plastic_isotropic_init(fileUnit)
endif endif
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
phase = phase + 1_pInt ! advance section counter phase = phase + 1_pInt ! advance section counter
if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then
instance = phase_plasticityInstance(phase) ! count instances of my constitutive law
allocate(param(instance)%outputID(phase_Noutput(phase))) ! allocate space for IDs of every requested output
endif
cycle ! skip to next line cycle ! skip to next line
endif endif
if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
p => param(instance) ! shorthand pointer to parameter object of my constitutive law
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key
@ -201,58 +189,58 @@ subroutine plastic_isotropic_init(fileUnit)
select case(outputtag) select case(outputtag)
case ('flowstress') case ('flowstress')
plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
param(instance)%outputID (plastic_isotropic_Noutput(instance)) = flowstress_ID
plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag
p%outputID = [p%outputID,flowstress_ID]
case ('strainrate') case ('strainrate')
plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
param(instance)%outputID (plastic_isotropic_Noutput(instance)) = strainrate_ID
plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag
p%outputID = [p%outputID,strainrate_ID]
end select end select
case ('/dilatation/') case ('/dilatation/')
param(instance)%dilatation = .true. p%dilatation = .true.
case ('tau0') case ('tau0')
param(instance)%tau0 = IO_floatValue(line,chunkPos,2_pInt) p%tau0 = IO_floatValue(line,chunkPos,2_pInt)
case ('gdot0') case ('gdot0')
param(instance)%gdot0 = IO_floatValue(line,chunkPos,2_pInt) p%gdot0 = IO_floatValue(line,chunkPos,2_pInt)
case ('n') case ('n')
param(instance)%n = IO_floatValue(line,chunkPos,2_pInt) p%n = IO_floatValue(line,chunkPos,2_pInt)
case ('h0') case ('h0')
param(instance)%h0 = IO_floatValue(line,chunkPos,2_pInt) p%h0 = IO_floatValue(line,chunkPos,2_pInt)
case ('h0_slope','slopelnrate') case ('h0_slope','slopelnrate')
param(instance)%h0_slopeLnRate = IO_floatValue(line,chunkPos,2_pInt) p%h0_slopeLnRate = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat') case ('tausat')
param(instance)%tausat = IO_floatValue(line,chunkPos,2_pInt) p%tausat = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfita') case ('tausat_sinhfita')
param(instance)%tausat_SinhFitA = IO_floatValue(line,chunkPos,2_pInt) p%tausat_SinhFitA = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitb') case ('tausat_sinhfitb')
param(instance)%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt) p%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitc') case ('tausat_sinhfitc')
param(instance)%tausat_SinhFitC = IO_floatValue(line,chunkPos,2_pInt) p%tausat_SinhFitC = IO_floatValue(line,chunkPos,2_pInt)
case ('tausat_sinhfitd') case ('tausat_sinhfitd')
param(instance)%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt) p%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt)
case ('a', 'w0') case ('a', 'w0')
param(instance)%a = IO_floatValue(line,chunkPos,2_pInt) p%a = IO_floatValue(line,chunkPos,2_pInt)
case ('taylorfactor') case ('taylorfactor')
param(instance)%fTaylor = IO_floatValue(line,chunkPos,2_pInt) p%fTaylor = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_flowstress') case ('atol_flowstress')
param(instance)%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt) p%aTolFlowstress = IO_floatValue(line,chunkPos,2_pInt)
case ('atol_shear') case ('atol_shear')
param(instance)%aTolShear = IO_floatValue(line,chunkPos,2_pInt) p%aTolShear = IO_floatValue(line,chunkPos,2_pInt)
case default case default
@ -261,25 +249,24 @@ subroutine plastic_isotropic_init(fileUnit)
enddo parsingFile enddo parsingFile
allocate(state(maxNinstance)) ! internal state aliases allocate(state(maxNinstance)) ! internal state aliases
allocate(state0(maxNinstance))
allocate(dotState(maxNinstance)) allocate(dotState(maxNinstance))
allocate(stateAbsTol(maxNinstance))
initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop over every plasticity initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop over every plasticity
myPhase: if (phase_plasticity(phase) == PLASTICITY_isotropic_ID) then ! isolate instances of own constitutive description myPhase: if (phase_plasticity(phase) == PLASTICITY_isotropic_ID) then ! isolate instances of own constitutive description
NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc) NipcMyPhase = count(material_phase == phase) ! number of own material points (including point components ipc)
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
p => param(instance)
extmsg = '' extmsg = ''
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity checks
if (param(instance)%aTolShear <= 0.0_pReal) param(instance)%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6 if (p%aTolShear <= 0.0_pReal) p%aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6
if (param(instance)%tau0 < 0.0_pReal) extmsg = trim(extmsg)//' tau0' if (p%tau0 < 0.0_pReal) extmsg = trim(extmsg)//' tau0'
if (param(instance)%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' if (p%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0'
if (param(instance)%n <= 0.0_pReal) extmsg = trim(extmsg)//' n' if (p%n <= 0.0_pReal) extmsg = trim(extmsg)//' n'
if (param(instance)%tausat <= 0.0_pReal) extmsg = trim(extmsg)//' tausat' if (p%tausat <= 0.0_pReal) extmsg = trim(extmsg)//' tausat'
if (param(instance)%a <= 0.0_pReal) extmsg = trim(extmsg)//' a' if (p%a <= 0.0_pReal) extmsg = trim(extmsg)//' a'
if (param(instance)%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//' taylorfactor' if (p%fTaylor <= 0.0_pReal) extmsg = trim(extmsg)//' taylorfactor'
if (param(instance)%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_flowstress' if (p%aTolFlowstress <= 0.0_pReal) extmsg = trim(extmsg)//' atol_flowstress'
if (extmsg /= '') then if (extmsg /= '') then
extmsg = trim(extmsg)//' ('//PLASTICITY_ISOTROPIC_label//')' ! prepare error message identifier extmsg = trim(extmsg)//' ('//PLASTICITY_ISOTROPIC_label//')' ! prepare error message identifier
call IO_error(211_pInt,ip=instance,ext_msg=extmsg) call IO_error(211_pInt,ip=instance,ext_msg=extmsg)
@ -287,7 +274,7 @@ subroutine plastic_isotropic_init(fileUnit)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Determine size of postResults array ! Determine size of postResults array
outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance) outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
select case(param(instance)%outputID(o)) select case(p%outputID(o))
case(flowstress_ID,strainrate_ID) case(flowstress_ID,strainrate_ID)
mySize = 1_pInt mySize = 1_pInt
case default case default
@ -331,31 +318,20 @@ subroutine plastic_isotropic_init(fileUnit)
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal) allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! globally required state aliases ! locally defined state aliases and initialization of state0 and aTolState
plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase)
plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases
state(instance)%flowstress => plasticState(phase)%state (1,1:NipcMyPhase) state(instance)%flowstress => plasticState(phase)%state (1,1:NipcMyPhase)
state0(instance)%flowstress => plasticState(phase)%state0 (1,1:NipcMyPhase)
dotState(instance)%flowstress => plasticState(phase)%dotState (1,1:NipcMyPhase) dotState(instance)%flowstress => plasticState(phase)%dotState (1,1:NipcMyPhase)
stateAbsTol(instance)%flowstress => plasticState(phase)%aTolState(1) plasticState(phase)%state0(1,1:NipcMyPhase) = p%tau0
plasticState(phase)%aTolState(1) = p%aTolFlowstress
state(instance)%accumulatedShear => plasticState(phase)%state (2,1:NipcMyPhase) state(instance)%accumulatedShear => plasticState(phase)%state (2,1:NipcMyPhase)
state0(instance)%accumulatedShear => plasticState(phase)%state0 (2,1:NipcMyPhase)
dotState(instance)%accumulatedShear => plasticState(phase)%dotState (2,1:NipcMyPhase) dotState(instance)%accumulatedShear => plasticState(phase)%dotState (2,1:NipcMyPhase)
stateAbsTol(instance)%accumulatedShear => plasticState(phase)%aTolState(2) plasticState(phase)%state0 (2,1:NipcMyPhase) = 0.0_pReal
plasticState(phase)%aTolState(2) = p%aTolShear
!-------------------------------------------------------------------------------------------------- ! global alias
! init state plasticState(phase)%slipRate => plasticState(phase)%dotState(2:2,1:NipcMyPhase)
state0(instance)%flowstress = param(instance)%tau0 plasticState(phase)%accumulatedSlip => plasticState(phase)%state (2:2,1:NipcMyPhase)
state0(instance)%accumulatedShear = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! init absolute state tolerances
stateAbsTol(instance)%flowstress = param(instance)%aTolFlowstress
stateAbsTol(instance)%accumulatedShear = param(instance)%aTolShear
endif myPhase endif myPhase
enddo initializeInstances enddo initializeInstances
@ -400,6 +376,8 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(tParameters), pointer :: p
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
@ -414,6 +392,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33) squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
@ -423,11 +402,11 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dTstar99 = 0.0_pReal dLp_dTstar99 = 0.0_pReal
else else
gamma_dot = param(instance)%gdot0 & gamma_dot = p%gdot0 &
* ( sqrt(1.5_pReal) * norm_Tstar_dev / param(instance)%fTaylor / state(instance)%flowstress(of) ) & * ( sqrt(1.5_pReal) * norm_Tstar_dev / p%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n **p%n
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/param(instance)%fTaylor Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/p%fTaylor
if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) & .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
@ -441,13 +420,13 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp ! Calculation of the tangent of Lp
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * & dLp_dTstar_3333(k,l,m,n) = (p%n-1.0_pReal) * &
Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
dLp_dTstar99 = math_Plain3333to99(gamma_dot / param(instance)%fTaylor * & dLp_dTstar99 = math_Plain3333to99(gamma_dot / p%fTaylor * &
dLp_dTstar_3333 / norm_Tstar_dev) dLp_dTstar_3333 / norm_Tstar_dev)
end if end if
end subroutine plastic_isotropic_LpAndItsTangent end subroutine plastic_isotropic_LpAndItsTangent
@ -479,9 +458,11 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,e
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(tParameters), pointer :: p
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal) :: & real(pReal) :: &
gamma_dot, & !< strainrate gamma_dot, & !< strainrate
norm_Tstar_sph, & !< euclidean norm of Tstar_sph norm_Tstar_sph, & !< euclidean norm of Tstar_sph
squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph
@ -491,34 +472,34 @@ real(pReal) :: &
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v)) ! spherical part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33) squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33)
norm_Tstar_sph = sqrt(squarenorm_Tstar_sph) norm_Tstar_sph = sqrt(squarenorm_Tstar_sph)
if (param(instance)%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero if (p%dilatation .and. norm_Tstar_sph > 0.0_pReal) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero
gamma_dot = param(instance)%gdot0 & gamma_dot = p%gdot0 &
* (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) & * (sqrt(1.5_pReal) * norm_Tstar_sph / p%fTaylor / state(instance)%flowstress(of) ) &
**param(instance)%n **p%n
Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/p%fTaylor
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Li ! Calculation of the tangent of Li
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * & dLi_dTstar_3333(k,l,m,n) = (p%n-1.0_pReal) * &
Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * & dLi_dTstar_3333 = gamma_dot / p%fTaylor * &
dLi_dTstar_3333 / norm_Tstar_sph dLi_dTstar_3333 / norm_Tstar_sph
else else
Li = 0.0_pReal Li = 0.0_pReal
dLi_dTstar_3333 = 0.0_pReal dLi_dTstar_3333 = 0.0_pReal
endif endif
end subroutine plastic_isotropic_LiAndItsTangent
end subroutine plastic_isotropic_LiAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -541,6 +522,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(tParameters), pointer :: p
real(pReal), dimension(6) :: & real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: & real(pReal) :: &
@ -554,10 +536,11 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress ! norm of (deviatoric) 2nd Piola-Kirchhoff stress
if (param(instance)%dilatation) then if (p%dilatation) then
norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v)) norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
else else
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
@ -566,26 +549,26 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
end if end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! strain rate ! strain rate
gamma_dot = param(instance)%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & gamma_dot = p%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
/ &!----------------------------------------------------------------------------------- / &!-----------------------------------------------------------------------------------
(param(instance)%fTaylor*state(instance)%flowstress(of) ))**param(instance)%n (p%fTaylor*state(instance)%flowstress(of) ))**p%n
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! hardening coefficient ! hardening coefficient
if (abs(gamma_dot) > 1e-12_pReal) then if (abs(gamma_dot) > 1e-12_pReal) then
if (dEq0(param(instance)%tausat_SinhFitA)) then if (dEq0(p%tausat_SinhFitA)) then
saturation = param(instance)%tausat saturation = p%tausat
else else
saturation = param(instance)%tausat & saturation = p%tausat &
+ asinh( (gamma_dot / param(instance)%tausat_SinhFitA& + asinh( (gamma_dot / p%tausat_SinhFitA&
)**(1.0_pReal / param(instance)%tausat_SinhFitD)& )**(1.0_pReal / p%tausat_SinhFitD)&
)**(1.0_pReal / param(instance)%tausat_SinhFitC) & )**(1.0_pReal / p%tausat_SinhFitC) &
/ ( param(instance)%tausat_SinhFitB & / ( p%tausat_SinhFitB &
* (gamma_dot / param(instance)%gdot0)**(1.0_pReal / param(instance)%n) & * (gamma_dot / p%gdot0)**(1.0_pReal / p%n) &
) )
endif endif
hardening = ( param(instance)%h0 + param(instance)%h0_slopeLnRate * log(gamma_dot) ) & hardening = ( p%h0 + p%h0_slopeLnRate * log(gamma_dot) ) &
* abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**param(instance)%a & * abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**p%a &
* sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation) * sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation)
else else
hardening = 0.0_pReal hardening = 0.0_pReal
@ -614,6 +597,9 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(tParameters), pointer :: p
real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
plastic_isotropic_postResults plastic_isotropic_postResults
@ -629,10 +615,11 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember of = phasememberAt(ipc,ip,el) ! phasememberAt should be tackled by material and be renamed to material_phasemember
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
p => param(instance)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress ! norm of (deviatoric) 2nd Piola-Kirchhoff stress
if (param(instance)%dilatation) then if (p%dilatation) then
norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v)) norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
else else
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
@ -644,15 +631,15 @@ function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
plastic_isotropic_postResults = 0.0_pReal plastic_isotropic_postResults = 0.0_pReal
outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance) outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
select case(param(instance)%outputID(o)) select case(p%outputID(o))
case (flowstress_ID) case (flowstress_ID)
plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of) plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of)
c = c + 1_pInt c = c + 1_pInt
case (strainrate_ID) case (strainrate_ID)
plastic_isotropic_postResults(c+1_pInt) = & plastic_isotropic_postResults(c+1_pInt) = &
param(instance)%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & p%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v &
/ &!---------------------------------------------------------------------------------- / &!----------------------------------------------------------------------------------
(param(instance)%fTaylor * state(instance)%flowstress(of)) ) ** param(instance)%n (p%fTaylor * state(instance)%flowstress(of)) ) ** p%n
c = c + 1_pInt c = c + 1_pInt
end select end select
enddo outputsLoop enddo outputsLoop

View File

@ -22,7 +22,7 @@ module spectral_mech_basic
private private
character (len=*), parameter, public :: & character (len=*), parameter, public :: &
DAMASK_spectral_SolverBasicPETSC_label = 'basicpetsc' DAMASK_spectral_SolverBasicPETSC_label = 'basic'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types