Merge branch 'development' of magit1.mpie.de:damask/DAMASK into development
This commit is contained in:
commit
7c39d01b35
|
@ -251,8 +251,6 @@ subroutine CPFEM_general(age, dt)
|
||||||
crystallite_Tstar0_v, &
|
crystallite_Tstar0_v, &
|
||||||
crystallite_Tstar_v
|
crystallite_Tstar_v
|
||||||
use homogenization, only: &
|
use homogenization, only: &
|
||||||
materialpoint_F, &
|
|
||||||
materialpoint_F0, &
|
|
||||||
materialpoint_stressAndItsTangent, &
|
materialpoint_stressAndItsTangent, &
|
||||||
materialpoint_postResults
|
materialpoint_postResults
|
||||||
use IO, only: &
|
use IO, only: &
|
||||||
|
|
|
@ -546,7 +546,7 @@ function IO_hybridIA(Nast,ODFfileName)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! math module is not available
|
! math module is not available
|
||||||
real(pReal), parameter :: PI = 3.14159265358979323846264338327950288419716939937510_pReal
|
real(pReal), parameter :: PI = 3.141592653589793_pReal
|
||||||
real(pReal), parameter :: INRAD = PI/180.0_pReal
|
real(pReal), parameter :: INRAD = PI/180.0_pReal
|
||||||
|
|
||||||
integer(pInt) :: i,j,bin,NnonZero,Nset,Nreps,reps,phi1,Phi,phi2
|
integer(pInt) :: i,j,bin,NnonZero,Nset,Nreps,reps,phi1,Phi,phi2
|
||||||
|
@ -666,7 +666,7 @@ function IO_hybridIA(Nast,ODFfileName)
|
||||||
else
|
else
|
||||||
prob = 0.0_pReal
|
prob = 0.0_pReal
|
||||||
endif
|
endif
|
||||||
dV_V(phi2,Phi,phi1) = prob*dg_0*sin((Phi-1.0_pReal+center)*deltas(2))
|
dV_V(phi2,Phi,phi1) = prob*dg_0*sin((real(Phi-1_pInt,pReal)+center)*deltas(2))
|
||||||
enddo; enddo; enddo
|
enddo; enddo; enddo
|
||||||
close(FILEUNIT)
|
close(FILEUNIT)
|
||||||
dV_V = dV_V/sum_dV_V ! normalize to 1
|
dV_V = dV_V/sum_dV_V ! normalize to 1
|
||||||
|
@ -713,7 +713,7 @@ function IO_hybridIA(Nast,ODFfileName)
|
||||||
do i=1_pInt,Nast
|
do i=1_pInt,Nast
|
||||||
if (i < Nast) then
|
if (i < Nast) then
|
||||||
call random_number(rnd)
|
call random_number(rnd)
|
||||||
j = nint(rnd*(Nreps-i)+i+0.5_pReal,pInt)
|
j = nint(rnd*real(Nreps-i,pReal)+real(i,pReal)+0.5_pReal,pInt)
|
||||||
else
|
else
|
||||||
j = i
|
j = i
|
||||||
endif
|
endif
|
||||||
|
|
|
@ -604,8 +604,6 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
|
||||||
use material, only: &
|
use material, only: &
|
||||||
phase_plasticity, &
|
phase_plasticity, &
|
||||||
material_phase, &
|
material_phase, &
|
||||||
material_homog, &
|
|
||||||
phaseAt, phasememberAt, &
|
|
||||||
phase_kinematics, &
|
phase_kinematics, &
|
||||||
phase_Nkinematics, &
|
phase_Nkinematics, &
|
||||||
PLASTICITY_isotropic_ID, &
|
PLASTICITY_isotropic_ID, &
|
||||||
|
|
|
@ -197,8 +197,6 @@ subroutine crystallite_init
|
||||||
nMax, & !< maximum number of ip neighbors
|
nMax, & !< maximum number of ip neighbors
|
||||||
myNcomponents, & !< number of components at current IP
|
myNcomponents, & !< number of components at current IP
|
||||||
section = 0_pInt, &
|
section = 0_pInt, &
|
||||||
j, &
|
|
||||||
p, &
|
|
||||||
mySize
|
mySize
|
||||||
|
|
||||||
character(len=65536) :: &
|
character(len=65536) :: &
|
||||||
|
@ -511,7 +509,8 @@ end subroutine crystallite_init
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine crystallite_stressAndItsTangent(updateJaco)
|
subroutine crystallite_stressAndItsTangent(updateJaco)
|
||||||
use prec, only: &
|
use prec, only: &
|
||||||
tol_math_check
|
tol_math_check, &
|
||||||
|
dNeq
|
||||||
use numerics, only: &
|
use numerics, only: &
|
||||||
subStepMinCryst, &
|
subStepMinCryst, &
|
||||||
subStepSizeCryst, &
|
subStepSizeCryst, &
|
||||||
|
@ -803,7 +802,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
||||||
endif
|
endif
|
||||||
else
|
else
|
||||||
subFracIntermediate = maxval(crystallite_subFrac, mask=.not.crystallite_localPlasticity)
|
subFracIntermediate = maxval(crystallite_subFrac, mask=.not.crystallite_localPlasticity)
|
||||||
if (abs(subFracIntermediate) > tiny(0.0_pReal)) then
|
if (dNeq(subFracIntermediate,0.0_pReal)) then
|
||||||
crystallite_neighborEnforcedCutback = .false. ! look for ips that require a cutback because of a nonconverged neighbor
|
crystallite_neighborEnforcedCutback = .false. ! look for ips that require a cutback because of a nonconverged neighbor
|
||||||
!$OMP PARALLEL
|
!$OMP PARALLEL
|
||||||
!$OMP DO PRIVATE(neighboring_e,neighboring_i)
|
!$OMP DO PRIVATE(neighboring_e,neighboring_i)
|
||||||
|
@ -844,7 +843,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
|
||||||
!$OMP DO PRIVATE(neighboring_e,neighboring_i)
|
!$OMP DO PRIVATE(neighboring_e,neighboring_i)
|
||||||
do e = FEsolving_execElem(1),FEsolving_execElem(2)
|
do e = FEsolving_execElem(1),FEsolving_execElem(2)
|
||||||
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
|
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
|
||||||
if (.not. crystallite_localPlasticity(1,i,e) .and. abs(crystallite_subFrac(1,i,e)) > tiny(0.0_pReal)) then
|
if (.not. crystallite_localPlasticity(1,i,e) .and. dNeq(crystallite_subFrac(1,i,e),0.0_pReal)) then
|
||||||
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))
|
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))
|
||||||
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
|
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
|
||||||
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
|
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
|
||||||
|
@ -3352,7 +3351,8 @@ end subroutine crystallite_integrateStateFPI
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
logical function crystallite_stateJump(ipc,ip,el)
|
logical function crystallite_stateJump(ipc,ip,el)
|
||||||
use prec, only: &
|
use prec, only: &
|
||||||
prec_isNaN
|
prec_isNaN, &
|
||||||
|
dNeq
|
||||||
use debug, only: &
|
use debug, only: &
|
||||||
debug_level, &
|
debug_level, &
|
||||||
debug_crystallite, &
|
debug_crystallite, &
|
||||||
|
@ -3404,7 +3404,7 @@ logical function crystallite_stateJump(ipc,ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
#ifndef _OPENMP
|
#ifndef _OPENMP
|
||||||
if (any(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c) /= 0.0_pReal) &
|
if (any(dNeq(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c),0.0_pReal)) &
|
||||||
.and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
|
.and. iand(debug_level(debug_crystallite), 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) &
|
||||||
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
|
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
|
||||||
|
@ -3459,7 +3459,8 @@ logical function crystallite_integrateStress(&
|
||||||
)
|
)
|
||||||
use prec, only: pLongInt, &
|
use prec, only: pLongInt, &
|
||||||
tol_math_check, &
|
tol_math_check, &
|
||||||
prec_isNaN
|
prec_isNaN, &
|
||||||
|
dEq
|
||||||
use numerics, only: nStress, &
|
use numerics, only: nStress, &
|
||||||
aTol_crystalliteStress, &
|
aTol_crystalliteStress, &
|
||||||
rTol_crystalliteStress, &
|
rTol_crystalliteStress, &
|
||||||
|
@ -3607,7 +3608,7 @@ logical function crystallite_integrateStress(&
|
||||||
!* inversion of Fp_current...
|
!* inversion of Fp_current...
|
||||||
|
|
||||||
invFp_current = math_inv33(Fp_current)
|
invFp_current = math_inv33(Fp_current)
|
||||||
if (all(abs(invFp_current) <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
|
failedInversionFp: if (all(dEq(invFp_current,0.0_pReal))) then
|
||||||
#ifndef _OPENMP
|
#ifndef _OPENMP
|
||||||
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
|
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
|
||||||
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',&
|
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',&
|
||||||
|
@ -3617,13 +3618,13 @@ logical function crystallite_integrateStress(&
|
||||||
endif
|
endif
|
||||||
#endif
|
#endif
|
||||||
return
|
return
|
||||||
endif
|
endif failedInversionFp
|
||||||
A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
|
A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
|
||||||
|
|
||||||
!* inversion of Fi_current...
|
!* inversion of Fi_current...
|
||||||
|
|
||||||
invFi_current = math_inv33(Fi_current)
|
invFi_current = math_inv33(Fi_current)
|
||||||
if (all(abs(invFi_current) <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
|
failedInversionFi: if (all(dEq(invFi_current,0.0_pReal))) then
|
||||||
#ifndef _OPENMP
|
#ifndef _OPENMP
|
||||||
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
|
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
|
||||||
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip ipc ',&
|
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip ipc ',&
|
||||||
|
@ -3633,7 +3634,7 @@ logical function crystallite_integrateStress(&
|
||||||
endif
|
endif
|
||||||
#endif
|
#endif
|
||||||
return
|
return
|
||||||
endif
|
endif failedInversionFi
|
||||||
|
|
||||||
!* start LpLoop with normal step length
|
!* start LpLoop with normal step length
|
||||||
|
|
||||||
|
@ -3883,7 +3884,7 @@ logical function crystallite_integrateStress(&
|
||||||
invFp_new = math_mul33x33(invFp_current,B)
|
invFp_new = math_mul33x33(invFp_current,B)
|
||||||
invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
|
invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
|
||||||
Fp_new = math_inv33(invFp_new)
|
Fp_new = math_inv33(invFp_new)
|
||||||
if (all(abs(Fp_new)<= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
|
failedInversionInvFp: if (all(dEq(Fp_new,0.0_pReal))) then
|
||||||
#ifndef _OPENMP
|
#ifndef _OPENMP
|
||||||
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
|
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
|
||||||
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip ipc ',&
|
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip ipc ',&
|
||||||
|
@ -3895,7 +3896,7 @@ logical function crystallite_integrateStress(&
|
||||||
endif
|
endif
|
||||||
#endif
|
#endif
|
||||||
return
|
return
|
||||||
endif
|
endif failedInversionInvFp
|
||||||
Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi_new) ! calc resulting Fe
|
Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi_new) ! calc resulting Fe
|
||||||
|
|
||||||
!* calculate 1st Piola-Kirchhoff stress
|
!* calculate 1st Piola-Kirchhoff stress
|
||||||
|
@ -4116,7 +4117,7 @@ function crystallite_postResults(ipc, ip, el)
|
||||||
mySize = 1_pInt
|
mySize = 1_pInt
|
||||||
detF = math_det33(crystallite_partionedF(1:3,1:3,ipc,ip,el)) ! V_current = det(F) * V_reference
|
detF = math_det33(crystallite_partionedF(1:3,1:3,ipc,ip,el)) ! V_current = det(F) * V_reference
|
||||||
crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) &
|
crystallite_postResults(c+1) = detF * mesh_ipVolume(ip,el) &
|
||||||
/ homogenization_Ngrains(mesh_element(3,el)) ! grain volume (not fraction but absolute)
|
/ real(homogenization_Ngrains(mesh_element(3,el)),pReal) ! grain volume (not fraction but absolute)
|
||||||
case (orientation_ID)
|
case (orientation_ID)
|
||||||
mySize = 4_pInt
|
mySize = 4_pInt
|
||||||
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,ipc,ip,el) ! grain orientation as quaternion
|
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,ipc,ip,el) ! grain orientation as quaternion
|
||||||
|
|
|
@ -223,7 +223,7 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
|
||||||
use material, only: &
|
use material, only: &
|
||||||
homogenization_Ngrains, &
|
homogenization_Ngrains, &
|
||||||
mappingHomogenization, &
|
mappingHomogenization, &
|
||||||
phaseAt, phasememberAt, &
|
phaseAt, &
|
||||||
phase_source, &
|
phase_source, &
|
||||||
phase_Nsources, &
|
phase_Nsources, &
|
||||||
SOURCE_damage_isoBrittle_ID, &
|
SOURCE_damage_isoBrittle_ID, &
|
||||||
|
@ -280,8 +280,8 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
phiDot = phiDot/homogenization_Ngrains(mappingHomogenization(2,ip,el))
|
phiDot = phiDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
|
||||||
dPhiDot_dPhi = dPhiDot_dPhi/homogenization_Ngrains(mappingHomogenization(2,ip,el))
|
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
|
||||||
|
|
||||||
end subroutine damage_local_getSourceAndItsTangent
|
end subroutine damage_local_getSourceAndItsTangent
|
||||||
|
|
||||||
|
|
|
@ -184,7 +184,7 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
|
||||||
use material, only: &
|
use material, only: &
|
||||||
homogenization_Ngrains, &
|
homogenization_Ngrains, &
|
||||||
mappingHomogenization, &
|
mappingHomogenization, &
|
||||||
phaseAt, phasememberAt, &
|
phaseAt, &
|
||||||
phase_source, &
|
phase_source, &
|
||||||
phase_Nsources, &
|
phase_Nsources, &
|
||||||
SOURCE_damage_isoBrittle_ID, &
|
SOURCE_damage_isoBrittle_ID, &
|
||||||
|
@ -241,8 +241,8 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
phiDot = phiDot/homogenization_Ngrains(mappingHomogenization(2,ip,el))
|
phiDot = phiDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
|
||||||
dPhiDot_dPhi = dPhiDot_dPhi/homogenization_Ngrains(mappingHomogenization(2,ip,el))
|
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
|
||||||
|
|
||||||
end subroutine damage_nonlocal_getSourceAndItsTangent
|
end subroutine damage_nonlocal_getSourceAndItsTangent
|
||||||
|
|
||||||
|
@ -279,9 +279,7 @@ function damage_nonlocal_getDiffusion33(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
damage_nonlocal_getDiffusion33 = &
|
damage_nonlocal_getDiffusion33 = &
|
||||||
charLength*charLength* &
|
charLength**2_pInt*damage_nonlocal_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
|
||||||
damage_nonlocal_getDiffusion33/ &
|
|
||||||
homogenization_Ngrains(homog)
|
|
||||||
|
|
||||||
end function damage_nonlocal_getDiffusion33
|
end function damage_nonlocal_getDiffusion33
|
||||||
|
|
||||||
|
@ -310,7 +308,8 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
|
||||||
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el))
|
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_DamageMobility(material_phase(ipc,ip,el))
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
damage_nonlocal_getMobility = damage_nonlocal_getMobility /homogenization_Ngrains(mesh_element(3,el))
|
damage_nonlocal_getMobility = damage_nonlocal_getMobility/&
|
||||||
|
real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
|
|
||||||
end function damage_nonlocal_getMobility
|
end function damage_nonlocal_getMobility
|
||||||
|
|
||||||
|
|
|
@ -46,10 +46,10 @@ module debug
|
||||||
integer(pInt),protected, dimension(debug_maxNtype+2_pInt), public :: & ! specific ones, and 2 for "all" and "other"
|
integer(pInt),protected, dimension(debug_maxNtype+2_pInt), public :: & ! specific ones, and 2 for "all" and "other"
|
||||||
debug_level = 0_pInt
|
debug_level = 0_pInt
|
||||||
|
|
||||||
integer(pInt), public :: &
|
integer(pLongInt), public :: &
|
||||||
debug_cumLpCalls = 0_pInt, & !< total number of calls to LpAndItsTangent
|
debug_cumLpCalls = 0_pLongInt, & !< total number of calls to LpAndItsTangent
|
||||||
debug_cumDeltaStateCalls = 0_pInt, & !< total number of calls to deltaState
|
debug_cumDeltaStateCalls = 0_pLongInt, & !< total number of calls to deltaState
|
||||||
debug_cumDotStateCalls = 0_pInt !< total number of calls to dotState
|
debug_cumDotStateCalls = 0_pLongInt !< total number of calls to dotState
|
||||||
|
|
||||||
integer(pInt), protected, public :: &
|
integer(pInt), protected, public :: &
|
||||||
debug_e = 1_pInt, &
|
debug_e = 1_pInt, &
|
||||||
|
@ -67,6 +67,7 @@ module debug
|
||||||
debug_jacobianMaxLocation = 0_pInt, &
|
debug_jacobianMaxLocation = 0_pInt, &
|
||||||
debug_jacobianMinLocation = 0_pInt
|
debug_jacobianMinLocation = 0_pInt
|
||||||
|
|
||||||
|
|
||||||
integer(pInt), dimension(:), allocatable, public :: &
|
integer(pInt), dimension(:), allocatable, public :: &
|
||||||
debug_CrystalliteLoopDistribution, & !< distribution of crystallite cutbacks
|
debug_CrystalliteLoopDistribution, & !< distribution of crystallite cutbacks
|
||||||
debug_MaterialpointStateLoopDistribution, &
|
debug_MaterialpointStateLoopDistribution, &
|
||||||
|
|
|
@ -386,6 +386,8 @@ end subroutine homogenization_RGC_partitionDeformation
|
||||||
! "happy" with result
|
! "happy" with result
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
|
function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
|
||||||
|
use prec, only: &
|
||||||
|
dEq
|
||||||
use debug, only: &
|
use debug, only: &
|
||||||
debug_level, &
|
debug_level, &
|
||||||
debug_homogenization,&
|
debug_homogenization,&
|
||||||
|
@ -441,10 +443,10 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
|
||||||
real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
|
real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
|
||||||
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
|
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
|
||||||
|
|
||||||
if(abs(dt) < tiny(0.0_pReal)) then ! zero time step
|
zeroTimeStep: if(dEq(dt,0.0_pReal)) then
|
||||||
homogenization_RGC_updateState = .true. ! pretend everything is fine and return
|
homogenization_RGC_updateState = .true. ! pretend everything is fine and return
|
||||||
return
|
return
|
||||||
endif
|
endif zeroTimeStep
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! get the dimension of the cluster (grains and interfaces)
|
! get the dimension of the cluster (grains and interfaces)
|
||||||
|
|
|
@ -77,7 +77,6 @@ subroutine kinematics_thermal_expansion_init(fileUnit)
|
||||||
integer(pInt) :: maxNinstance,phase,instance,kinematics
|
integer(pInt) :: maxNinstance,phase,instance,kinematics
|
||||||
character(len=65536) :: &
|
character(len=65536) :: &
|
||||||
tag = '', &
|
tag = '', &
|
||||||
output = '', &
|
|
||||||
line = ''
|
line = ''
|
||||||
|
|
||||||
mainProcess: if (worldrank == 0) then
|
mainProcess: if (worldrank == 0) then
|
||||||
|
|
|
@ -1232,6 +1232,8 @@ end subroutine material_parseTexture
|
||||||
!! calculates the volume of the grains and deals with texture components and hybridIA
|
!! calculates the volume of the grains and deals with texture components and hybridIA
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine material_populateGrains
|
subroutine material_populateGrains
|
||||||
|
use prec, only: &
|
||||||
|
dEq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
math_RtoEuler, &
|
math_RtoEuler, &
|
||||||
math_EulerToR, &
|
math_EulerToR, &
|
||||||
|
@ -1371,7 +1373,7 @@ subroutine material_populateGrains
|
||||||
else
|
else
|
||||||
forall (i = 1_pInt:FE_Nips(t)) & ! loop over IPs
|
forall (i = 1_pInt:FE_Nips(t)) & ! loop over IPs
|
||||||
volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = &
|
volumeOfGrain(grain+(i-1)*dGrains+1_pInt:grain+i*dGrains) = &
|
||||||
mesh_ipVolume(i,e)/dGrains ! assign IPvolume/Ngrains@IP to all grains of IP
|
mesh_ipVolume(i,e)/real(dGrains,pReal) ! assign IPvolume/Ngrains@IP to all grains of IP
|
||||||
grain = grain + FE_Nips(t) * dGrains ! wind forward by Nips*Ngrains@IP
|
grain = grain + FE_Nips(t) * dGrains ! wind forward by Nips*Ngrains@IP
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
|
@ -1393,7 +1395,7 @@ subroutine material_populateGrains
|
||||||
|
|
||||||
NgrainsOfConstituent = 0_pInt ! reset counter of grains per constituent
|
NgrainsOfConstituent = 0_pInt ! reset counter of grains per constituent
|
||||||
forall (i = 1_pInt:myNconstituents) &
|
forall (i = 1_pInt:myNconstituents) &
|
||||||
NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro) * myNgrains, pInt) ! do rounding integer conversion
|
NgrainsOfConstituent(i) = nint(microstructure_fraction(i,micro)*real(myNgrains,pReal),pInt)! do rounding integer conversion
|
||||||
do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong?
|
do while (sum(NgrainsOfConstituent) /= myNgrains) ! total grain count over constituents wrong?
|
||||||
sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change
|
sgn = sign(1_pInt, myNgrains - sum(NgrainsOfConstituent)) ! direction of required change
|
||||||
extreme = 0.0_pReal
|
extreme = 0.0_pReal
|
||||||
|
@ -1434,24 +1436,24 @@ subroutine material_populateGrains
|
||||||
! ...has texture components
|
! ...has texture components
|
||||||
if (texture_ODFfile(textureID) == '') then
|
if (texture_ODFfile(textureID) == '') then
|
||||||
gauss: do t = 1_pInt,texture_Ngauss(textureID) ! loop over Gauss components
|
gauss: do t = 1_pInt,texture_Ngauss(textureID) ! loop over Gauss components
|
||||||
do g = 1_pInt,int(myNorientations*texture_Gauss(5,t,textureID),pInt) ! loop over required grain count
|
do g = 1_pInt,int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID),pInt) ! loop over required grain count
|
||||||
orientationOfGrain(:,grain+constituentGrain+g) = &
|
orientationOfGrain(:,grain+constituentGrain+g) = &
|
||||||
math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
|
math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
|
||||||
texture_Gauss( 4,t,textureID))
|
texture_Gauss( 4,t,textureID))
|
||||||
enddo
|
enddo
|
||||||
constituentGrain = &
|
constituentGrain = &
|
||||||
constituentGrain + int(myNorientations*texture_Gauss(5,t,textureID)) ! advance counter for grains of current constituent
|
constituentGrain + int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID)) ! advance counter for grains of current constituent
|
||||||
enddo gauss
|
enddo gauss
|
||||||
|
|
||||||
fiber: do t = 1_pInt,texture_Nfiber(textureID) ! loop over fiber components
|
fiber: do t = 1_pInt,texture_Nfiber(textureID) ! loop over fiber components
|
||||||
do g = 1_pInt,int(myNorientations*texture_Fiber(6,t,textureID),pInt) ! loop over required grain count
|
do g = 1_pInt,int(real(myNorientations,pReal)*texture_Fiber(6,t,textureID),pInt) ! loop over required grain count
|
||||||
orientationOfGrain(:,grain+constituentGrain+g) = &
|
orientationOfGrain(:,grain+constituentGrain+g) = &
|
||||||
math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
|
math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
|
||||||
texture_Fiber(3:4,t,textureID),&
|
texture_Fiber(3:4,t,textureID),&
|
||||||
texture_Fiber( 5,t,textureID))
|
texture_Fiber( 5,t,textureID))
|
||||||
enddo
|
enddo
|
||||||
constituentGrain = &
|
constituentGrain = &
|
||||||
constituentGrain + int(myNorientations*texture_fiber(6,t,textureID),pInt) ! advance counter for grains of current constituent
|
constituentGrain + int(real(myNorientations,pReal)*texture_fiber(6,t,textureID),pInt) ! advance counter for grains of current constituent
|
||||||
enddo fiber
|
enddo fiber
|
||||||
|
|
||||||
random: do constituentGrain = constituentGrain+1_pInt,myNorientations ! fill remainder with random
|
random: do constituentGrain = constituentGrain+1_pInt,myNorientations ! fill remainder with random
|
||||||
|
@ -1462,7 +1464,7 @@ subroutine material_populateGrains
|
||||||
else
|
else
|
||||||
orientationOfGrain(1:3,grain+1_pInt:grain+myNorientations) = &
|
orientationOfGrain(1:3,grain+1_pInt:grain+myNorientations) = &
|
||||||
IO_hybridIA(myNorientations,texture_ODFfile(textureID))
|
IO_hybridIA(myNorientations,texture_ODFfile(textureID))
|
||||||
if (all(orientationOfGrain(1:3,grain+1_pInt) == -1.0_pReal)) call IO_error(156_pInt)
|
if (all(dEq(orientationOfGrain(1:3,grain+1_pInt),-1.0_pReal))) call IO_error(156_pInt)
|
||||||
endif
|
endif
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
@ -1501,7 +1503,7 @@ subroutine material_populateGrains
|
||||||
|
|
||||||
do j = 1_pInt,NgrainsOfConstituent(i)-1_pInt ! walk thru grains of current constituent
|
do j = 1_pInt,NgrainsOfConstituent(i)-1_pInt ! walk thru grains of current constituent
|
||||||
call random_number(rnd)
|
call random_number(rnd)
|
||||||
t = nint(rnd*(NgrainsOfConstituent(i)-j)+j+0.5_pReal,pInt) ! select a grain in remaining list
|
t = nint(rnd*real(NgrainsOfConstituent(i)-j,pReal)+real(j,pReal)+0.5_pReal,pInt) ! select a grain in remaining list
|
||||||
m = phaseOfGrain(grain+t) ! exchange current with random
|
m = phaseOfGrain(grain+t) ! exchange current with random
|
||||||
phaseOfGrain(grain+t) = phaseOfGrain(grain+j)
|
phaseOfGrain(grain+t) = phaseOfGrain(grain+j)
|
||||||
phaseOfGrain(grain+j) = m
|
phaseOfGrain(grain+j) = m
|
||||||
|
@ -1539,7 +1541,7 @@ subroutine material_populateGrains
|
||||||
randomOrder = math_range(microstructure_maxNconstituents) ! start out with ordered sequence of constituents
|
randomOrder = math_range(microstructure_maxNconstituents) ! start out with ordered sequence of constituents
|
||||||
call random_number(rndArray) ! as many rnd numbers as (max) constituents
|
call random_number(rndArray) ! as many rnd numbers as (max) constituents
|
||||||
do j = 1_pInt, myNconstituents - 1_pInt ! loop over constituents ...
|
do j = 1_pInt, myNconstituents - 1_pInt ! loop over constituents ...
|
||||||
r = nint(rndArray(j)*(myNconstituents-j)+j+0.5_pReal,pInt) ! ... select one in remaining list
|
r = nint(rndArray(j)*real(myNconstituents-j,pReal)+real(j,pReal)+0.5_pReal,pInt) ! ... select one in remaining list
|
||||||
c = randomOrder(r) ! ... call it "c"
|
c = randomOrder(r) ! ... call it "c"
|
||||||
randomOrder(r) = randomOrder(j) ! ... and exchange with present position in constituent list
|
randomOrder(r) = randomOrder(j) ! ... and exchange with present position in constituent list
|
||||||
grain = sum(NgrainsOfConstituent(1:c-1_pInt)) ! figure out actual starting index in overall/consecutive grain population
|
grain = sum(NgrainsOfConstituent(1:c-1_pInt)) ! figure out actual starting index in overall/consecutive grain population
|
||||||
|
|
|
@ -13,10 +13,10 @@ module math
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
private
|
private
|
||||||
real(pReal), parameter, public :: PI = 3.14159265358979323846264338327950288419716939937510_pReal !< ratio of a circle's circumference to its diameter
|
real(pReal), parameter, public :: PI = 3.141592653589793_pReal !< ratio of a circle's circumference to its diameter
|
||||||
real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree
|
real(pReal), parameter, public :: INDEG = 180.0_pReal/PI !< conversion from radian into degree
|
||||||
real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian
|
real(pReal), parameter, public :: INRAD = PI/180.0_pReal !< conversion from degree into radian
|
||||||
complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)* PI !< Re(0.0), Im(2xPi)
|
complex(pReal), parameter, public :: TWOPIIMG = (0.0_pReal,2.0_pReal)*(PI,0.0_pReal) !< Re(0.0), Im(2xPi)
|
||||||
|
|
||||||
real(pReal), dimension(3,3), parameter, public :: &
|
real(pReal), dimension(3,3), parameter, public :: &
|
||||||
MATH_I3 = reshape([&
|
MATH_I3 = reshape([&
|
||||||
|
@ -723,6 +723,8 @@ end function math_transpose33
|
||||||
! returns all zeroes if not possible, i.e. if det close to zero
|
! returns all zeroes if not possible, i.e. if det close to zero
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
pure function math_inv33(A)
|
pure function math_inv33(A)
|
||||||
|
use prec, only: &
|
||||||
|
dNeq
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
real(pReal),dimension(3,3),intent(in) :: A
|
real(pReal),dimension(3,3),intent(in) :: A
|
||||||
|
@ -735,7 +737,7 @@ pure function math_inv33(A)
|
||||||
|
|
||||||
DetA = A(1,1) * math_inv33(1,1) + A(1,2) * math_inv33(2,1) + A(1,3) * math_inv33(3,1)
|
DetA = A(1,1) * math_inv33(1,1) + A(1,2) * math_inv33(2,1) + A(1,3) * math_inv33(3,1)
|
||||||
|
|
||||||
if (abs(DetA) > tiny(DetA)) then ! use a real threshold here
|
if (dNeq(DetA,0.0_pReal)) then
|
||||||
math_inv33(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2)
|
math_inv33(1,2) = -A(1,2) * A(3,3) + A(1,3) * A(3,2)
|
||||||
math_inv33(2,2) = A(1,1) * A(3,3) - A(1,3) * A(3,1)
|
math_inv33(2,2) = A(1,1) * A(3,3) - A(1,3) * A(3,1)
|
||||||
math_inv33(3,2) = -A(1,1) * A(3,2) + A(1,2) * A(3,1)
|
math_inv33(3,2) = -A(1,1) * A(3,2) + A(1,2) * A(3,1)
|
||||||
|
@ -759,6 +761,8 @@ end function math_inv33
|
||||||
! returns error if not possible, i.e. if det close to zero
|
! returns error if not possible, i.e. if det close to zero
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
pure subroutine math_invert33(A, InvA, DetA, error)
|
pure subroutine math_invert33(A, InvA, DetA, error)
|
||||||
|
use prec, only: &
|
||||||
|
dEq
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
logical, intent(out) :: error
|
logical, intent(out) :: error
|
||||||
|
@ -772,7 +776,7 @@ pure subroutine math_invert33(A, InvA, DetA, error)
|
||||||
|
|
||||||
DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1)
|
DetA = A(1,1) * InvA(1,1) + A(1,2) * InvA(2,1) + A(1,3) * InvA(3,1)
|
||||||
|
|
||||||
if (abs(DetA) <= tiny(DetA)) then
|
if (dEq(DetA,0.0_pReal)) then
|
||||||
InvA = 0.0_pReal
|
InvA = 0.0_pReal
|
||||||
error = .true.
|
error = .true.
|
||||||
else
|
else
|
||||||
|
@ -1277,6 +1281,8 @@ end function math_qNorm
|
||||||
!> @brief quaternion inversion
|
!> @brief quaternion inversion
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
pure function math_qInv(Q)
|
pure function math_qInv(Q)
|
||||||
|
use prec, only: &
|
||||||
|
dNeq
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
real(pReal), dimension(4), intent(in) :: Q
|
real(pReal), dimension(4), intent(in) :: Q
|
||||||
|
@ -1286,8 +1292,7 @@ pure function math_qInv(Q)
|
||||||
math_qInv = 0.0_pReal
|
math_qInv = 0.0_pReal
|
||||||
|
|
||||||
squareNorm = math_qDot(Q,Q)
|
squareNorm = math_qDot(Q,Q)
|
||||||
if (abs(squareNorm) > tiny(squareNorm)) &
|
if (dNeq(squareNorm,0.0_pReal)) math_qInv = math_qConj(Q) / squareNorm
|
||||||
math_qInv = math_qConj(Q) / squareNorm
|
|
||||||
|
|
||||||
end function math_qInv
|
end function math_qInv
|
||||||
|
|
||||||
|
@ -2091,6 +2096,8 @@ end function math_eigenvectorBasisSym33
|
||||||
!> @brief rotational part from polar decomposition of 33 tensor m
|
!> @brief rotational part from polar decomposition of 33 tensor m
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
function math_rotationalPart33(m)
|
function math_rotationalPart33(m)
|
||||||
|
use prec, only: &
|
||||||
|
dEq
|
||||||
use IO, only: &
|
use IO, only: &
|
||||||
IO_warning
|
IO_warning
|
||||||
|
|
||||||
|
@ -2102,12 +2109,12 @@ function math_rotationalPart33(m)
|
||||||
U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m))
|
U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m))
|
||||||
Uinv = math_inv33(U)
|
Uinv = math_inv33(U)
|
||||||
|
|
||||||
if (all(abs(Uinv) <= tiny(Uinv))) then ! math_inv33 returns zero when failed, avoid floating point equality comparison
|
inversionFailed: if (all(dEq(Uinv,0.0_pReal))) then
|
||||||
math_rotationalPart33 = math_I3
|
math_rotationalPart33 = math_I3
|
||||||
call IO_warning(650_pInt)
|
call IO_warning(650_pInt)
|
||||||
else
|
else inversionFailed
|
||||||
math_rotationalPart33 = math_mul33x33(m,Uinv)
|
math_rotationalPart33 = math_mul33x33(m,Uinv)
|
||||||
endif
|
endif inversionFailed
|
||||||
|
|
||||||
end function math_rotationalPart33
|
end function math_rotationalPart33
|
||||||
|
|
||||||
|
|
|
@ -963,7 +963,7 @@ subroutine mesh_build_ipCoordinates
|
||||||
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
|
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
|
||||||
myCoords = myCoords + mesh_cellnode(1:3,mesh_cell(n,i,e))
|
myCoords = myCoords + mesh_cellnode(1:3,mesh_cell(n,i,e))
|
||||||
enddo
|
enddo
|
||||||
mesh_ipCoordinates(1:3,i,e) = myCoords / FE_NcellnodesPerCell(c)
|
mesh_ipCoordinates(1:3,i,e) = myCoords / real(FE_NcellnodesPerCell(c),pReal)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END PARALLEL DO
|
!$OMP END PARALLEL DO
|
||||||
|
@ -990,7 +990,7 @@ pure function mesh_cellCenterCoordinates(ip,el)
|
||||||
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
|
do n = 1_pInt,FE_NcellnodesPerCell(c) ! loop over cell nodes in this cell
|
||||||
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el))
|
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates + mesh_cellnode(1:3,mesh_cell(n,ip,el))
|
||||||
enddo
|
enddo
|
||||||
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / FE_NcellnodesPerCell(c)
|
mesh_cellCenterCoordinates = mesh_cellCenterCoordinates / real(FE_NcellnodesPerCell(c),pReal)
|
||||||
|
|
||||||
end function mesh_cellCenterCoordinates
|
end function mesh_cellCenterCoordinates
|
||||||
|
|
||||||
|
@ -3070,7 +3070,6 @@ use IO, only: &
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
integer(pInt), intent(in) :: fileUnit
|
integer(pInt), intent(in) :: fileUnit
|
||||||
|
|
||||||
#ifndef Spectral
|
#ifndef Spectral
|
||||||
integer(pInt), allocatable, dimension(:) :: chunkPos
|
integer(pInt), allocatable, dimension(:) :: chunkPos
|
||||||
integer(pInt) chunk, Nchunks
|
integer(pInt) chunk, Nchunks
|
||||||
|
@ -3082,10 +3081,9 @@ use IO, only: &
|
||||||
mesh_periodicSurface = .true.
|
mesh_periodicSurface = .true.
|
||||||
#else
|
#else
|
||||||
mesh_periodicSurface = .false.
|
mesh_periodicSurface = .false.
|
||||||
#ifdef Marc4DAMASK
|
#if defined(Marc4DAMASK)
|
||||||
keyword = '$damask'
|
keyword = '$damask'
|
||||||
#endif
|
#elif defined(Abaqus)
|
||||||
#ifdef Abaqus
|
|
||||||
keyword = '**damask'
|
keyword = '**damask'
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
|
@ -3693,6 +3691,7 @@ integer(pInt) function FE_mapElemtype(what)
|
||||||
'c3d20t')
|
'c3d20t')
|
||||||
FE_mapElemtype = 13_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
|
FE_mapElemtype = 13_pInt ! Three-dimensional Arbitrarily Distorted quadratic hexahedral
|
||||||
case default
|
case default
|
||||||
|
FE_mapElemtype = -1_pInt ! error return
|
||||||
call IO_error(error_ID=190_pInt,ext_msg=IO_lc(what))
|
call IO_error(error_ID=190_pInt,ext_msg=IO_lc(what))
|
||||||
end select
|
end select
|
||||||
|
|
||||||
|
@ -3701,6 +3700,7 @@ end function FE_mapElemtype
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief find face-matching element of same type
|
!> @brief find face-matching element of same type
|
||||||
|
!> @details currently not used, check if needed for HDF5 output, otherwise delete
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace)
|
subroutine mesh_faceMatch(elem, face ,matchingElem, matchingFace)
|
||||||
|
|
||||||
|
|
|
@ -981,7 +981,8 @@ end subroutine plastic_disloUCLA_LpAndItsTangent
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
|
subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
|
||||||
use prec, only: &
|
use prec, only: &
|
||||||
tol_math_check
|
tol_math_check, &
|
||||||
|
dEq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
pi
|
pi
|
||||||
use material, only: &
|
use material, only: &
|
||||||
|
@ -1119,7 +1120,7 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
|
||||||
!* Dipole formation
|
!* Dipole formation
|
||||||
EdgeDipMinDistance = &
|
EdgeDipMinDistance = &
|
||||||
plastic_disloUCLA_CEdgeDipMinDistance(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance)
|
plastic_disloUCLA_CEdgeDipMinDistance(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance)
|
||||||
if (abs(tau_slip_pos) <= tiny(0.0_pReal)) then
|
if (dEq(tau_slip_pos,0.0_pReal)) then
|
||||||
DotRhoDipFormation = 0.0_pReal
|
DotRhoDipFormation = 0.0_pReal
|
||||||
else
|
else
|
||||||
EdgeDipDistance = &
|
EdgeDipDistance = &
|
||||||
|
@ -1147,7 +1148,7 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
|
||||||
plastic_disloUCLA_CAtomicVolume(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance)**(3.0_pReal)
|
plastic_disloUCLA_CAtomicVolume(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance)**(3.0_pReal)
|
||||||
VacancyDiffusion = &
|
VacancyDiffusion = &
|
||||||
plastic_disloUCLA_D0(instance)*exp(-plastic_disloUCLA_Qsd(instance)/(kB*Temperature))
|
plastic_disloUCLA_D0(instance)*exp(-plastic_disloUCLA_Qsd(instance)/(kB*Temperature))
|
||||||
if (abs(tau_slip_pos) <= tiny(0.0_pReal)) then
|
if (dEq(tau_slip_pos,0.0_pReal)) then
|
||||||
DotRhoEdgeDipClimb = 0.0_pReal
|
DotRhoEdgeDipClimb = 0.0_pReal
|
||||||
else
|
else
|
||||||
ClimbVelocity = &
|
ClimbVelocity = &
|
||||||
|
@ -1180,7 +1181,8 @@ end subroutine plastic_disloUCLA_dotState
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el)
|
function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el)
|
||||||
use prec, only: &
|
use prec, only: &
|
||||||
tol_math_check
|
tol_math_check, &
|
||||||
|
dEq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
pi
|
pi
|
||||||
use material, only: &
|
use material, only: &
|
||||||
|
@ -1408,7 +1410,7 @@ function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el)
|
||||||
c = c + ns
|
c = c + ns
|
||||||
elseif(plastic_disloUCLA_outputID(o,instance) == stress_exponent_ID) then
|
elseif(plastic_disloUCLA_outputID(o,instance) == stress_exponent_ID) then
|
||||||
do j = 1_pInt, ns
|
do j = 1_pInt, ns
|
||||||
if (abs(gdot_slip_pos(j)+gdot_slip_neg(j))<=tiny(0.0_pReal)) then
|
if (dEq(gdot_slip_pos(j)+gdot_slip_neg(j),0.0_pReal)) then
|
||||||
plastic_disloUCLA_postResults(c+j) = 0.0_pReal
|
plastic_disloUCLA_postResults(c+j) = 0.0_pReal
|
||||||
else
|
else
|
||||||
plastic_disloUCLA_postResults(c+j) = (tau_slip_pos(j)+tau_slip_neg(j))/&
|
plastic_disloUCLA_postResults(c+j) = (tau_slip_pos(j)+tau_slip_neg(j))/&
|
||||||
|
|
|
@ -199,6 +199,9 @@ contains
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_dislotwin_init(fileUnit)
|
subroutine plastic_dislotwin_init(fileUnit)
|
||||||
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
|
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
|
||||||
|
use prec, only: &
|
||||||
|
dEq, &
|
||||||
|
dNeq
|
||||||
use debug, only: &
|
use debug, only: &
|
||||||
debug_level,&
|
debug_level,&
|
||||||
debug_constitutive,&
|
debug_constitutive,&
|
||||||
|
@ -749,8 +752,8 @@ subroutine plastic_dislotwin_init(fileUnit)
|
||||||
if (plastic_dislotwin_Qsd(instance) <= 0.0_pReal) &
|
if (plastic_dislotwin_Qsd(instance) <= 0.0_pReal) &
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')')
|
||||||
if (sum(plastic_dislotwin_Ntwin(:,instance)) > 0_pInt) then
|
if (sum(plastic_dislotwin_Ntwin(:,instance)) > 0_pInt) then
|
||||||
if (abs(plastic_dislotwin_SFE_0K(instance)) <= tiny(0.0_pReal) .and. &
|
if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. &
|
||||||
abs(plastic_dislotwin_dSFE_dT(instance)) <= tiny(0.0_pReal) .and. &
|
dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. &
|
||||||
lattice_structure(phase) == LATTICE_fcc_ID) &
|
lattice_structure(phase) == LATTICE_fcc_ID) &
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
|
||||||
if (plastic_dislotwin_aTolRho(instance) <= 0.0_pReal) &
|
if (plastic_dislotwin_aTolRho(instance) <= 0.0_pReal) &
|
||||||
|
@ -759,8 +762,8 @@ subroutine plastic_dislotwin_init(fileUnit)
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')')
|
||||||
endif
|
endif
|
||||||
if (sum(plastic_dislotwin_Ntrans(:,instance)) > 0_pInt) then
|
if (sum(plastic_dislotwin_Ntrans(:,instance)) > 0_pInt) then
|
||||||
if (abs(plastic_dislotwin_SFE_0K(instance)) <= tiny(0.0_pReal) .and. &
|
if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. &
|
||||||
abs(plastic_dislotwin_dSFE_dT(instance)) <= tiny(0.0_pReal) .and. &
|
dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. &
|
||||||
lattice_structure(phase) == LATTICE_fcc_ID) &
|
lattice_structure(phase) == LATTICE_fcc_ID) &
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
|
||||||
if (plastic_dislotwin_aTolTransFrac(instance) <= 0.0_pReal) &
|
if (plastic_dislotwin_aTolTransFrac(instance) <= 0.0_pReal) &
|
||||||
|
@ -773,8 +776,8 @@ subroutine plastic_dislotwin_init(fileUnit)
|
||||||
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
|
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
|
||||||
plastic_dislotwin_pShearBand(instance) <= 0.0_pReal) &
|
plastic_dislotwin_pShearBand(instance) <= 0.0_pReal) &
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')')
|
||||||
if (abs(plastic_dislotwin_dipoleFormationFactor(instance)) > tiny(0.0_pReal) .and. &
|
if (dNeq(plastic_dislotwin_dipoleFormationFactor(instance), 0.0_pReal) .and. &
|
||||||
plastic_dislotwin_dipoleFormationFactor(instance) /= 1.0_pReal) &
|
dNeq(plastic_dislotwin_dipoleFormationFactor(instance), 1.0_pReal)) &
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOTWIN_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOTWIN_label//')')
|
||||||
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
|
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
|
||||||
plastic_dislotwin_qShearBand(instance) <= 0.0_pReal) &
|
plastic_dislotwin_qShearBand(instance) <= 0.0_pReal) &
|
||||||
|
@ -1628,7 +1631,8 @@ end subroutine plastic_dislotwin_microstructure
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el)
|
subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el)
|
||||||
use prec, only: &
|
use prec, only: &
|
||||||
tol_math_check
|
tol_math_check, &
|
||||||
|
dNeq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
math_Plain3333to99, &
|
math_Plain3333to99, &
|
||||||
math_Mandel6to33, &
|
math_Mandel6to33, &
|
||||||
|
@ -1775,8 +1779,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! Shear banding (shearband) part
|
! Shear banding (shearband) part
|
||||||
if(abs(plastic_dislotwin_sbVelocity(instance)) > tiny(0.0_pReal) .and. &
|
if(dNeq(plastic_dislotwin_sbVelocity(instance), 0.0_pReal) .and. &
|
||||||
abs(plastic_dislotwin_sbResistance(instance)) > tiny(0.0_pReal)) then
|
dNeq(plastic_dislotwin_sbResistance(instance),0.0_pReal)) then
|
||||||
gdot_sb = 0.0_pReal
|
gdot_sb = 0.0_pReal
|
||||||
dgdot_dtausb = 0.0_pReal
|
dgdot_dtausb = 0.0_pReal
|
||||||
call math_eigenValuesVectorsSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error)
|
call math_eigenValuesVectorsSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error)
|
||||||
|
@ -1942,7 +1946,8 @@ end subroutine plastic_dislotwin_LpAndItsTangent
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
|
subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
|
||||||
use prec, only: &
|
use prec, only: &
|
||||||
tol_math_check
|
tol_math_check, &
|
||||||
|
dEq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
pi
|
pi
|
||||||
use material, only: &
|
use material, only: &
|
||||||
|
@ -2043,7 +2048,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
|
||||||
!* Dipole formation
|
!* Dipole formation
|
||||||
EdgeDipMinDistance = &
|
EdgeDipMinDistance = &
|
||||||
plastic_dislotwin_CEdgeDipMinDistance(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)
|
plastic_dislotwin_CEdgeDipMinDistance(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)
|
||||||
if (abs(tau_slip(j)) <= tiny(0.0_pReal)) then
|
if (dEq(tau_slip(j),0.0_pReal)) then
|
||||||
DotRhoDipFormation = 0.0_pReal
|
DotRhoDipFormation = 0.0_pReal
|
||||||
else
|
else
|
||||||
EdgeDipDistance = &
|
EdgeDipDistance = &
|
||||||
|
@ -2071,10 +2076,10 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
|
||||||
plastic_dislotwin_CAtomicVolume(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)**(3.0_pReal)
|
plastic_dislotwin_CAtomicVolume(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)**(3.0_pReal)
|
||||||
VacancyDiffusion = &
|
VacancyDiffusion = &
|
||||||
plastic_dislotwin_D0(instance)*exp(-plastic_dislotwin_Qsd(instance)/(kB*Temperature))
|
plastic_dislotwin_D0(instance)*exp(-plastic_dislotwin_Qsd(instance)/(kB*Temperature))
|
||||||
if (abs(tau_slip(j)) <= tiny(0.0_pReal)) then
|
if (dEq(tau_slip(j),0.0_pReal)) then
|
||||||
DotRhoEdgeDipClimb = 0.0_pReal
|
DotRhoEdgeDipClimb = 0.0_pReal
|
||||||
else
|
else
|
||||||
if (EdgeDipDistance-EdgeDipMinDistance <= tiny(0.0_pReal)) then
|
if (dEq(EdgeDipDistance-EdgeDipMinDistance,0.0_pReal)) then
|
||||||
DotRhoEdgeDipClimb = 0.0_pReal
|
DotRhoEdgeDipClimb = 0.0_pReal
|
||||||
else
|
else
|
||||||
ClimbVelocity = 3.0_pReal*lattice_mu(ph)*VacancyDiffusion*AtomicVolume/ &
|
ClimbVelocity = 3.0_pReal*lattice_mu(ph)*VacancyDiffusion*AtomicVolume/ &
|
||||||
|
@ -2189,7 +2194,8 @@ end subroutine plastic_dislotwin_dotState
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
|
function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
|
||||||
use prec, only: &
|
use prec, only: &
|
||||||
tol_math_check
|
tol_math_check, &
|
||||||
|
dEq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
pi, &
|
pi, &
|
||||||
math_Mandel6to33, &
|
math_Mandel6to33, &
|
||||||
|
@ -2504,11 +2510,8 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
|
||||||
endif
|
endif
|
||||||
|
|
||||||
!* Stress exponent
|
!* Stress exponent
|
||||||
if (abs(gdot_slip(j))<=tiny(0.0_pReal)) then
|
plastic_dislotwin_postResults(c+j) = &
|
||||||
plastic_dislotwin_postResults(c+j) = 0.0_pReal
|
merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq(gdot_slip(j),0.0_pReal))
|
||||||
else
|
|
||||||
plastic_dislotwin_postResults(c+j) = (tau/gdot_slip(j))*dgdot_dtauslip
|
|
||||||
endif
|
|
||||||
enddo ; enddo
|
enddo ; enddo
|
||||||
c = c + ns
|
c = c + ns
|
||||||
case (sb_eigenvalues_ID)
|
case (sb_eigenvalues_ID)
|
||||||
|
|
|
@ -524,6 +524,8 @@ end subroutine plastic_isotropic_LiAndItsTangent
|
||||||
!> @brief calculates the rate of change of microstructure
|
!> @brief calculates the rate of change of microstructure
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
|
subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
|
||||||
|
use prec, only: &
|
||||||
|
dEq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
math_mul6x6
|
math_mul6x6
|
||||||
use material, only: &
|
use material, only: &
|
||||||
|
@ -570,7 +572,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! hardening coefficient
|
! hardening coefficient
|
||||||
if (abs(gamma_dot) > 1e-12_pReal) then
|
if (abs(gamma_dot) > 1e-12_pReal) then
|
||||||
if (abs(param(instance)%tausat_SinhFitA) <= tiny(0.0_pReal)) then
|
if (dEq(param(instance)%tausat_SinhFitA,0.0_pReal)) then
|
||||||
saturation = param(instance)%tausat
|
saturation = param(instance)%tausat
|
||||||
else
|
else
|
||||||
saturation = ( param(instance)%tausat &
|
saturation = ( param(instance)%tausat &
|
||||||
|
|
|
@ -1549,6 +1549,8 @@ end subroutine plastic_nonlocal_aTolState
|
||||||
!> @brief calculates quantities characterizing the microstructure
|
!> @brief calculates quantities characterizing the microstructure
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_nonlocal_microstructure(Fe, Fp, ip, el)
|
subroutine plastic_nonlocal_microstructure(Fe, Fp, ip, el)
|
||||||
|
use prec, only: &
|
||||||
|
dEq
|
||||||
use IO, only: &
|
use IO, only: &
|
||||||
IO_error
|
IO_error
|
||||||
use math, only: &
|
use math, only: &
|
||||||
|
@ -1792,7 +1794,7 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance))
|
||||||
- neighbor_rhoExcess(c,s,neighbors(2))
|
- neighbor_rhoExcess(c,s,neighbors(2))
|
||||||
enddo
|
enddo
|
||||||
invConnections = math_inv33(connections)
|
invConnections = math_inv33(connections)
|
||||||
if (all(abs(invConnections) <= tiny(0.0_pReal))) & ! check for failed in version (math_inv33 returns 0) and avoid floating point equality comparison
|
if (all(dEq(invConnections,0.0_pReal))) &
|
||||||
call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error')
|
call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error')
|
||||||
rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), &
|
rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), &
|
||||||
math_mul33x3(invConnections,rhoExcessDifferences))
|
math_mul33x3(invConnections,rhoExcessDifferences))
|
||||||
|
@ -2200,6 +2202,8 @@ end subroutine plastic_nonlocal_LpAndItsTangent
|
||||||
!> @brief (instantaneous) incremental change of microstructure
|
!> @brief (instantaneous) incremental change of microstructure
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_nonlocal_deltaState(Tstar_v,ip,el)
|
subroutine plastic_nonlocal_deltaState(Tstar_v,ip,el)
|
||||||
|
use prec, only: &
|
||||||
|
dNeq
|
||||||
use debug, only: debug_level, &
|
use debug, only: debug_level, &
|
||||||
debug_constitutive, &
|
debug_constitutive, &
|
||||||
debug_levelBasic, &
|
debug_levelBasic, &
|
||||||
|
@ -2322,8 +2326,8 @@ dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) / (4.0_pReal * pi * abs
|
||||||
|
|
||||||
|
|
||||||
forall (c = 1_pInt:2_pInt)
|
forall (c = 1_pInt:2_pInt)
|
||||||
where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
|
where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
|
||||||
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
|
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
|
||||||
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
|
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
|
||||||
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
|
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
|
||||||
dUpper(1:ns,c))
|
dUpper(1:ns,c))
|
||||||
|
@ -2335,7 +2339,7 @@ deltaDUpper = dUpper - dUpperOld
|
||||||
!*** dissociation by stress increase
|
!*** dissociation by stress increase
|
||||||
deltaRhoDipole2SingleStress = 0.0_pReal
|
deltaRhoDipole2SingleStress = 0.0_pReal
|
||||||
forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal .and. &
|
forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal .and. &
|
||||||
abs(dUpperOld(s,c) - dLower(s,c)) > tiny(0.0_pReal)) &
|
dNeq(dUpperOld(s,c) - dLower(s,c),0.0_pReal)) &
|
||||||
deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) &
|
deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) &
|
||||||
/ (dUpperOld(s,c) - dLower(s,c))
|
/ (dUpperOld(s,c) - dLower(s,c))
|
||||||
|
|
||||||
|
@ -2382,7 +2386,9 @@ end subroutine plastic_nonlocal_deltaState
|
||||||
subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
|
subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
|
||||||
timestep,subfrac, ip,el)
|
timestep,subfrac, ip,el)
|
||||||
|
|
||||||
use prec, only: DAMASK_NaN
|
use prec, only: DAMASK_NaN, &
|
||||||
|
dNeq, &
|
||||||
|
dEq
|
||||||
use numerics, only: numerics_integrationMode, &
|
use numerics, only: numerics_integrationMode, &
|
||||||
numerics_timeSyncing
|
numerics_timeSyncing
|
||||||
use IO, only: IO_error
|
use IO, only: IO_error
|
||||||
|
@ -2616,8 +2622,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) &
|
||||||
dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
|
dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
|
||||||
/ (4.0_pReal * pi * abs(tau))
|
/ (4.0_pReal * pi * abs(tau))
|
||||||
forall (c = 1_pInt:2_pInt)
|
forall (c = 1_pInt:2_pInt)
|
||||||
where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
|
where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
|
||||||
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
|
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
|
||||||
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
|
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
|
||||||
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
|
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
|
||||||
dUpper(1:ns,c))
|
dUpper(1:ns,c))
|
||||||
|
@ -2760,8 +2766,7 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
|
||||||
endif
|
endif
|
||||||
|
|
||||||
if (considerEnteringFlux) then
|
if (considerEnteringFlux) then
|
||||||
if(numerics_timeSyncing .and. (subfrac(1_pInt,neighbor_ip,neighbor_el) /= subfrac(1_pInt,ip,el))) &
|
if(numerics_timeSyncing .and. (dNeq(subfrac(1,neighbor_ip,neighbor_el),subfrac(1,ip,el)))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal
|
||||||
then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal
|
|
||||||
forall (s = 1:ns, t = 1_pInt:4_pInt)
|
forall (s = 1:ns, t = 1_pInt:4_pInt)
|
||||||
|
|
||||||
neighbor_v(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no)
|
neighbor_v(s,t) = plasticState(np)%state0(iV (s,t,neighbor_instance),no)
|
||||||
|
@ -2830,11 +2835,11 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
|
||||||
my_rhoSgl = rhoSgl
|
my_rhoSgl = rhoSgl
|
||||||
my_v = v
|
my_v = v
|
||||||
if(numerics_timeSyncing) then
|
if(numerics_timeSyncing) then
|
||||||
if (abs(subfrac(1_pInt,ip,el))<= tiny(0.0_pReal)) then
|
if (dEq(subfrac(1_pInt,ip,el),0.0_pReal)) then
|
||||||
my_rhoSgl = rhoSgl0
|
my_rhoSgl = rhoSgl0
|
||||||
my_v = v0
|
my_v = v0
|
||||||
elseif (neighbor_n > 0_pInt) then
|
elseif (neighbor_n > 0_pInt) then
|
||||||
if (abs(subfrac(1_pInt,neighbor_ip,neighbor_el))<= tiny(0.0_pReal)) then
|
if (dEq(subfrac(1_pInt,neighbor_ip,neighbor_el),0.0_pReal)) then
|
||||||
my_rhoSgl = rhoSgl0
|
my_rhoSgl = rhoSgl0
|
||||||
my_v = v0
|
my_v = v0
|
||||||
endif
|
endif
|
||||||
|
@ -3078,13 +3083,11 @@ slipDirection(1:3,1:ns) = lattice_sd(1:3, slipSystemLattice(1:ns,instance), ph)
|
||||||
!*** start out fully compatible
|
!*** start out fully compatible
|
||||||
|
|
||||||
my_compatibility = 0.0_pReal
|
my_compatibility = 0.0_pReal
|
||||||
forall(s1 = 1_pInt:ns) &
|
forall(s1 = 1_pInt:ns) my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal
|
||||||
my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal
|
|
||||||
|
|
||||||
|
|
||||||
!*** Loop thrugh neighbors and check whether there is any my_compatibility.
|
!*** Loop thrugh neighbors and check whether there is any my_compatibility.
|
||||||
|
|
||||||
do n = 1_pInt,Nneighbors
|
neighbors: do n = 1_pInt,Nneighbors
|
||||||
neighbor_e = mesh_ipNeighborhood(1,n,i,e)
|
neighbor_e = mesh_ipNeighborhood(1,n,i,e)
|
||||||
neighbor_i = mesh_ipNeighborhood(2,n,i,e)
|
neighbor_i = mesh_ipNeighborhood(2,n,i,e)
|
||||||
|
|
||||||
|
@ -3093,8 +3096,7 @@ do n = 1_pInt,Nneighbors
|
||||||
!* Set surface transmissivity to the value specified in the material.config
|
!* Set surface transmissivity to the value specified in the material.config
|
||||||
|
|
||||||
if (neighbor_e <= 0_pInt .or. neighbor_i <= 0_pInt) then
|
if (neighbor_e <= 0_pInt .or. neighbor_i <= 0_pInt) then
|
||||||
forall(s1 = 1_pInt:ns) &
|
forall(s1 = 1_pInt:ns) my_compatibility(1:2,s1,s1,n) = sqrt(surfaceTransmissivity(instance))
|
||||||
my_compatibility(1:2,s1,s1,n) = sqrt(surfaceTransmissivity(instance))
|
|
||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
@ -3107,10 +3109,8 @@ do n = 1_pInt,Nneighbors
|
||||||
|
|
||||||
neighbor_phase = material_phase(1,neighbor_i,neighbor_e)
|
neighbor_phase = material_phase(1,neighbor_i,neighbor_e)
|
||||||
if (neighbor_phase /= ph) then
|
if (neighbor_phase /= ph) then
|
||||||
if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph)) then
|
if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(ph))&
|
||||||
forall(s1 = 1_pInt:ns) &
|
forall(s1 = 1_pInt:ns) my_compatibility(1:2,s1,s1,n) = 0.0_pReal
|
||||||
my_compatibility(1:2,s1,s1,n) = 0.0_pReal ! = sqrt(0.0)
|
|
||||||
endif
|
|
||||||
cycle
|
cycle
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
@ -3141,33 +3141,33 @@ do n = 1_pInt,Nneighbors
|
||||||
else
|
else
|
||||||
absoluteMisorientation = lattice_qDisorientation(orientation(1:4,1,i,e), &
|
absoluteMisorientation = lattice_qDisorientation(orientation(1:4,1,i,e), &
|
||||||
orientation(1:4,1,neighbor_i,neighbor_e)) ! no symmetry
|
orientation(1:4,1,neighbor_i,neighbor_e)) ! no symmetry
|
||||||
do s1 = 1_pInt,ns ! my slip systems
|
mySlipSystems: do s1 = 1_pInt,ns
|
||||||
do s2 = 1_pInt,ns ! my neighbor's slip systems
|
neighborSlipSystems: do s2 = 1_pInt,ns
|
||||||
my_compatibility(1,s2,s1,n) = math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2))) &
|
my_compatibility(1,s2,s1,n) = math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2))) &
|
||||||
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
|
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
|
||||||
my_compatibility(2,s2,s1,n) = abs(math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2)))) &
|
my_compatibility(2,s2,s1,n) = abs(math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2)))) &
|
||||||
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
|
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
|
||||||
enddo
|
enddo neighborSlipSystems
|
||||||
|
|
||||||
my_compatibilitySum = 0.0_pReal
|
my_compatibilitySum = 0.0_pReal
|
||||||
belowThreshold = .true.
|
belowThreshold = .true.
|
||||||
do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns)))
|
do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns)))
|
||||||
thresholdValue = maxval(my_compatibility(2,1:ns,s1,n), belowThreshold(1:ns)) ! screws always positive
|
thresholdValue = maxval(my_compatibility(2,1:ns,s1,n), belowThreshold(1:ns)) ! screws always positive
|
||||||
nThresholdValues = real(count(my_compatibility(2,1:ns,s1,n) == thresholdValue),pReal)
|
nThresholdValues = real(count(my_compatibility(2,1:ns,s1,n) >= thresholdValue),pReal)
|
||||||
where (my_compatibility(2,1:ns,s1,n) >= thresholdValue) &
|
where (my_compatibility(2,1:ns,s1,n) >= thresholdValue) &
|
||||||
belowThreshold(1:ns) = .false.
|
belowThreshold(1:ns) = .false.
|
||||||
if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
|
if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
|
||||||
where (abs(my_compatibility(1:2,1:ns,s1,n)) == thresholdValue) & ! MD: rather check below threshold?
|
where (abs(my_compatibility(1:2,1:ns,s1,n)) >= thresholdValue) & ! MD: rather check below threshold?
|
||||||
my_compatibility(1:2,1:ns,s1,n) = sign((1.0_pReal - my_compatibilitySum) &
|
my_compatibility(1:2,1:ns,s1,n) = sign((1.0_pReal - my_compatibilitySum) &
|
||||||
/ nThresholdValues, my_compatibility(1:2,1:ns,s1,n))
|
/ nThresholdValues, my_compatibility(1:2,1:ns,s1,n))
|
||||||
my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
|
my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
|
||||||
enddo
|
enddo
|
||||||
where (belowThreshold(1:ns)) my_compatibility(1,1:ns,s1,n) = 0.0_pReal
|
where (belowThreshold(1:ns)) my_compatibility(1,1:ns,s1,n) = 0.0_pReal
|
||||||
where (belowThreshold(1:ns)) my_compatibility(2,1:ns,s1,n) = 0.0_pReal
|
where (belowThreshold(1:ns)) my_compatibility(2,1:ns,s1,n) = 0.0_pReal
|
||||||
enddo ! my slip systems cycle
|
enddo mySlipSystems
|
||||||
endif
|
endif
|
||||||
|
|
||||||
enddo ! neighbor cycle
|
enddo neighbors
|
||||||
|
|
||||||
compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = my_compatibility
|
compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = my_compatibility
|
||||||
|
|
||||||
|
@ -3177,6 +3177,8 @@ end subroutine plastic_nonlocal_updateCompatibility
|
||||||
!* calculates quantities characterizing the microstructure *
|
!* calculates quantities characterizing the microstructure *
|
||||||
!*********************************************************************
|
!*********************************************************************
|
||||||
function plastic_nonlocal_dislocationstress(Fe, ip, el)
|
function plastic_nonlocal_dislocationstress(Fe, ip, el)
|
||||||
|
use prec, only: &
|
||||||
|
dEq
|
||||||
use math, only: math_mul33x33, &
|
use math, only: math_mul33x33, &
|
||||||
math_mul33x3, &
|
math_mul33x3, &
|
||||||
math_inv33, &
|
math_inv33, &
|
||||||
|
@ -3389,7 +3391,7 @@ if (.not. phase_localPlasticity(ph)) then
|
||||||
Rsquare = R * R
|
Rsquare = R * R
|
||||||
Rcube = Rsquare * R
|
Rcube = Rsquare * R
|
||||||
denominator = R * (R + flipSign * lambda)
|
denominator = R * (R + flipSign * lambda)
|
||||||
if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop
|
if (dEq(denominator,0.0_pReal)) exit ipLoop
|
||||||
|
|
||||||
sigma(1,1) = sigma(1,1) - real(side,pReal) &
|
sigma(1,1) = sigma(1,1) - real(side,pReal) &
|
||||||
* flipSign * z / denominator &
|
* flipSign * z / denominator &
|
||||||
|
@ -3434,7 +3436,7 @@ if (.not. phase_localPlasticity(ph)) then
|
||||||
Rsquare = R * R
|
Rsquare = R * R
|
||||||
Rcube = Rsquare * R
|
Rcube = Rsquare * R
|
||||||
denominator = R * (R + flipSign * lambda)
|
denominator = R * (R + flipSign * lambda)
|
||||||
if (abs(denominator)<= tiny(0.0_pReal)) exit ipLoop
|
if (dEq(denominator,0.0_pReal)) exit ipLoop
|
||||||
|
|
||||||
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z &
|
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z &
|
||||||
* (1.0_pReal - lattice_nu(ph)) / denominator &
|
* (1.0_pReal - lattice_nu(ph)) / denominator &
|
||||||
|
@ -3523,6 +3525,8 @@ end function plastic_nonlocal_dislocationstress
|
||||||
!> @brief return array of constitutive results
|
!> @brief return array of constitutive results
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
|
function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
|
||||||
|
use prec, only: &
|
||||||
|
dNeq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
math_mul6x6, &
|
math_mul6x6, &
|
||||||
math_mul33x3, &
|
math_mul33x3, &
|
||||||
|
@ -3639,8 +3643,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) &
|
||||||
dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
|
dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
|
||||||
/ (4.0_pReal * pi * abs(tau))
|
/ (4.0_pReal * pi * abs(tau))
|
||||||
forall (c = 1_pInt:2_pInt)
|
forall (c = 1_pInt:2_pInt)
|
||||||
where(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
|
where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
|
||||||
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)) >= tiny(0.0_pReal)) &
|
abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
|
||||||
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
|
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
|
||||||
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
|
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
|
||||||
dUpper(1:ns,c))
|
dUpper(1:ns,c))
|
||||||
|
|
|
@ -112,6 +112,8 @@ contains
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_phenoplus_init(fileUnit)
|
subroutine plastic_phenoplus_init(fileUnit)
|
||||||
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
|
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
|
||||||
|
use prec, only: &
|
||||||
|
dEq
|
||||||
use debug, only: &
|
use debug, only: &
|
||||||
debug_level, &
|
debug_level, &
|
||||||
debug_constitutive,&
|
debug_constitutive,&
|
||||||
|
@ -481,7 +483,7 @@ subroutine plastic_phenoplus_init(fileUnit)
|
||||||
if (any(plastic_phenoplus_tausat_slip(:,instance) <= 0.0_pReal .and. &
|
if (any(plastic_phenoplus_tausat_slip(:,instance) <= 0.0_pReal .and. &
|
||||||
plastic_phenoplus_Nslip(:,instance) > 0)) &
|
plastic_phenoplus_Nslip(:,instance) > 0)) &
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPLUS_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPLUS_label//')')
|
||||||
if (any(abs(plastic_phenoplus_a_slip(instance)) <= tiny(0.0_pReal) .and. &
|
if (any(dEq(plastic_phenoplus_a_slip(instance),0.0_pReal) .and. &
|
||||||
plastic_phenoplus_Nslip(:,instance) > 0)) &
|
plastic_phenoplus_Nslip(:,instance) > 0)) &
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPLUS_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPLUS_label//')')
|
||||||
if (any(plastic_phenoplus_tau0_twin(:,instance) < 0.0_pReal .and. &
|
if (any(plastic_phenoplus_tau0_twin(:,instance) < 0.0_pReal .and. &
|
||||||
|
@ -873,7 +875,7 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
|
||||||
ENDDO LOOPFINDNEISHEAR
|
ENDDO LOOPFINDNEISHEAR
|
||||||
|
|
||||||
!***calculate the average accumulative shear and use it as cutoff
|
!***calculate the average accumulative shear and use it as cutoff
|
||||||
avg_acshear_ne = SUM(ne_acshear)/ns
|
avg_acshear_ne = sum(ne_acshear)/real(ns,pReal)
|
||||||
|
|
||||||
!***
|
!***
|
||||||
IF (ph==neighbor_ph) THEN
|
IF (ph==neighbor_ph) THEN
|
||||||
|
@ -923,6 +925,8 @@ end subroutine plastic_phenoplus_microstructure
|
||||||
!> @brief calculates plastic velocity gradient and its tangent
|
!> @brief calculates plastic velocity gradient and its tangent
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
|
subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
|
||||||
|
use prec, only: &
|
||||||
|
dNeq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
math_Plain3333to99, &
|
math_Plain3333to99, &
|
||||||
math_Mandel6to33
|
math_Mandel6to33
|
||||||
|
@ -1038,7 +1042,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
|
||||||
(gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
|
(gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
|
||||||
|
|
||||||
! Calculation of the tangent of Lp
|
! Calculation of the tangent of Lp
|
||||||
if (abs(gdot_slip_pos) > tiny(0.0_pReal)) then
|
if (dNeq(gdot_slip_pos,0.0_pReal)) then
|
||||||
dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenoplus_n_slip(instance)/tau_slip_pos
|
dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenoplus_n_slip(instance)/tau_slip_pos
|
||||||
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_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
||||||
|
@ -1046,7 +1050,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
|
||||||
nonSchmid_tensor(m,n,1)
|
nonSchmid_tensor(m,n,1)
|
||||||
endif
|
endif
|
||||||
|
|
||||||
if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then
|
if (dNeq(gdot_slip_neg,0.0_pReal)) then
|
||||||
dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenoplus_n_slip(instance)/tau_slip_neg
|
dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenoplus_n_slip(instance)/tau_slip_neg
|
||||||
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_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
||||||
|
@ -1073,7 +1077,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
|
||||||
Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
|
Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
|
||||||
|
|
||||||
! Calculation of the tangent of Lp
|
! Calculation of the tangent of Lp
|
||||||
if (abs(gdot_twin) > tiny(0.0_pReal)) then
|
if (dNeq(gdot_twin,0.0_pReal)) then
|
||||||
dgdot_dtautwin = gdot_twin*plastic_phenoplus_n_twin(instance)/tau_twin
|
dgdot_dtautwin = gdot_twin*plastic_phenoplus_n_twin(instance)/tau_twin
|
||||||
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_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
||||||
|
|
|
@ -124,6 +124,8 @@ contains
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_phenopowerlaw_init(fileUnit)
|
subroutine plastic_phenopowerlaw_init(fileUnit)
|
||||||
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
|
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
|
||||||
|
use prec, only: &
|
||||||
|
dEq
|
||||||
use debug, only: &
|
use debug, only: &
|
||||||
debug_level, &
|
debug_level, &
|
||||||
debug_constitutive,&
|
debug_constitutive,&
|
||||||
|
@ -487,7 +489,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
|
||||||
if (any(plastic_phenopowerlaw_tausat_slip(:,instance) <= 0.0_pReal .and. &
|
if (any(plastic_phenopowerlaw_tausat_slip(:,instance) <= 0.0_pReal .and. &
|
||||||
plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
|
plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
|
||||||
if (any(abs(plastic_phenopowerlaw_a_slip(instance)) <= tiny(0.0_pReal) .and. &
|
if (any(dEq(plastic_phenopowerlaw_a_slip(instance),0.0_pReal) .and. &
|
||||||
plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
|
plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
|
||||||
call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
|
call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
|
||||||
if (any(plastic_phenopowerlaw_tau0_twin(:,instance) < 0.0_pReal .and. &
|
if (any(plastic_phenopowerlaw_tau0_twin(:,instance) < 0.0_pReal .and. &
|
||||||
|
@ -774,6 +776,8 @@ end subroutine plastic_phenopowerlaw_aTolState
|
||||||
!> @brief calculates plastic velocity gradient and its tangent
|
!> @brief calculates plastic velocity gradient and its tangent
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
|
subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
|
||||||
|
use prec, only: &
|
||||||
|
dNeq
|
||||||
use math, only: &
|
use math, only: &
|
||||||
math_Plain3333to99, &
|
math_Plain3333to99, &
|
||||||
math_Mandel6to33
|
math_Mandel6to33
|
||||||
|
@ -863,7 +867,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
|
||||||
(gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
|
(gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
|
||||||
|
|
||||||
! Calculation of the tangent of Lp
|
! Calculation of the tangent of Lp
|
||||||
if (abs(gdot_slip_pos) > tiny(0.0_pReal)) then
|
if (dNeq(gdot_slip_pos,0.0_pReal)) then
|
||||||
dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenopowerlaw_n_slip(instance)/tau_slip_pos
|
dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenopowerlaw_n_slip(instance)/tau_slip_pos
|
||||||
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_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
||||||
|
@ -871,7 +875,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
|
||||||
nonSchmid_tensor(m,n,1)
|
nonSchmid_tensor(m,n,1)
|
||||||
endif
|
endif
|
||||||
|
|
||||||
if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then
|
if (dNeq(gdot_slip_neg,0.0_pReal)) then
|
||||||
dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenopowerlaw_n_slip(instance)/tau_slip_neg
|
dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenopowerlaw_n_slip(instance)/tau_slip_neg
|
||||||
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_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
||||||
|
@ -898,7 +902,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
|
||||||
Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
|
Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
|
||||||
|
|
||||||
! Calculation of the tangent of Lp
|
! Calculation of the tangent of Lp
|
||||||
if (abs(gdot_twin) > tiny(0.0_pReal)) then
|
if (dNeq(gdot_twin,0.0_pReal)) then
|
||||||
dgdot_dtautwin = gdot_twin*plastic_phenopowerlaw_n_twin(instance)/tau_twin
|
dgdot_dtautwin = gdot_twin*plastic_phenopowerlaw_n_twin(instance)/tau_twin
|
||||||
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_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
|
||||||
|
@ -1009,8 +1013,9 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
|
||||||
enddo nonSchmidSystems
|
enddo nonSchmidSystems
|
||||||
gdot_slip(j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
|
gdot_slip(j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
|
||||||
((abs(tau_slip_pos)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) &
|
((abs(tau_slip_pos)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) &
|
||||||
+(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance))&
|
*sign(1.0_pReal,tau_slip_pos) &
|
||||||
*sign(1.0_pReal,tau_slip_pos)
|
+(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) &
|
||||||
|
*sign(1.0_pReal,tau_slip_neg))
|
||||||
enddo slipSystems1
|
enddo slipSystems1
|
||||||
enddo slipFamilies1
|
enddo slipFamilies1
|
||||||
|
|
||||||
|
|
|
@ -210,8 +210,7 @@ function porosity_phasefield_getFormationEnergy(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
porosity_phasefield_getFormationEnergy = &
|
porosity_phasefield_getFormationEnergy = &
|
||||||
porosity_phasefield_getFormationEnergy/ &
|
porosity_phasefield_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function porosity_phasefield_getFormationEnergy
|
end function porosity_phasefield_getFormationEnergy
|
||||||
|
|
||||||
|
@ -243,8 +242,7 @@ function porosity_phasefield_getSurfaceEnergy(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
porosity_phasefield_getSurfaceEnergy = &
|
porosity_phasefield_getSurfaceEnergy = &
|
||||||
porosity_phasefield_getSurfaceEnergy/ &
|
porosity_phasefield_getSurfaceEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function porosity_phasefield_getSurfaceEnergy
|
end function porosity_phasefield_getSurfaceEnergy
|
||||||
|
|
||||||
|
@ -308,7 +306,7 @@ subroutine porosity_phasefield_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi,
|
||||||
enddo
|
enddo
|
||||||
W_e = W_e + sum(abs(strain*math_mul66x6(C,strain)))
|
W_e = W_e + sum(abs(strain*math_mul66x6(C,strain)))
|
||||||
enddo
|
enddo
|
||||||
W_e = W_e/homogenization_Ngrains(homog)
|
W_e = W_e/real(homogenization_Ngrains(homog),pReal)
|
||||||
|
|
||||||
phiDot = 2.0_pReal*(1.0_pReal - phi)*(1.0_pReal - Cv)*(1.0_pReal - Cv) - &
|
phiDot = 2.0_pReal*(1.0_pReal - phi)*(1.0_pReal - Cv)*(1.0_pReal - Cv) - &
|
||||||
2.0_pReal*phi*(W_e + Cv*porosity_phasefield_getFormationEnergy(ip,el))/ &
|
2.0_pReal*phi*(W_e + Cv*porosity_phasefield_getFormationEnergy(ip,el))/ &
|
||||||
|
@ -350,8 +348,7 @@ function porosity_phasefield_getDiffusion33(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
porosity_phasefield_getDiffusion33 = &
|
porosity_phasefield_getDiffusion33 = &
|
||||||
porosity_phasefield_getDiffusion33/ &
|
porosity_phasefield_getDiffusion33/real(homogenization_Ngrains(homog),pReal)
|
||||||
homogenization_Ngrains(homog)
|
|
||||||
|
|
||||||
end function porosity_phasefield_getDiffusion33
|
end function porosity_phasefield_getDiffusion33
|
||||||
|
|
||||||
|
@ -377,10 +374,12 @@ real(pReal) function porosity_phasefield_getMobility(ip,el)
|
||||||
porosity_phasefield_getMobility = 0.0_pReal
|
porosity_phasefield_getMobility = 0.0_pReal
|
||||||
|
|
||||||
do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
|
do ipc = 1, homogenization_Ngrains(mesh_element(3,el))
|
||||||
porosity_phasefield_getMobility = porosity_phasefield_getMobility + lattice_PorosityMobility(material_phase(ipc,ip,el))
|
porosity_phasefield_getMobility = porosity_phasefield_getMobility &
|
||||||
|
+ lattice_PorosityMobility(material_phase(ipc,ip,el))
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
porosity_phasefield_getMobility = porosity_phasefield_getMobility/homogenization_Ngrains(mesh_element(3,el))
|
porosity_phasefield_getMobility = &
|
||||||
|
porosity_phasefield_getMobility/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
|
|
||||||
end function porosity_phasefield_getMobility
|
end function porosity_phasefield_getMobility
|
||||||
|
|
||||||
|
|
|
@ -21,10 +21,10 @@ module prec
|
||||||
#if (FLOAT==8)
|
#if (FLOAT==8)
|
||||||
integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300)
|
integer, parameter, public :: pReal = 8 !< floating point double precision (was selected_real_kind(15,300), number with 15 significant digits, up to 1e+-300)
|
||||||
#ifdef __INTEL_COMPILER
|
#ifdef __INTEL_COMPILER
|
||||||
real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
|
real(pReal), parameter, public :: DAMASK_NaN = Z'7FF8000000000000' !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html)
|
||||||
#endif
|
#endif
|
||||||
#ifdef __GFORTRAN__
|
#ifdef __GFORTRAN__
|
||||||
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000',pReal) !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html, copy can be found in documentation/Code/Fortran)
|
real(pReal), parameter, public :: DAMASK_NaN = real(Z'7FF8000000000000',pReal) !< quiet NaN for double precision (from http://www.hpc.unimelb.edu.au/doc/f90lrm/dfum_035.html)
|
||||||
#endif
|
#endif
|
||||||
#else
|
#else
|
||||||
NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION
|
NO SUITABLE PRECISION FOR REAL SELECTED, STOPPING COMPILATION
|
||||||
|
@ -115,7 +115,9 @@ module prec
|
||||||
prec_init, &
|
prec_init, &
|
||||||
prec_isNaN, &
|
prec_isNaN, &
|
||||||
dEq, &
|
dEq, &
|
||||||
dNeq
|
cEq, &
|
||||||
|
dNeq, &
|
||||||
|
cNeq
|
||||||
|
|
||||||
contains
|
contains
|
||||||
|
|
||||||
|
@ -180,28 +182,69 @@ end function prec_isNaN
|
||||||
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief equality comparison for double precision
|
!> @brief equality comparison for float with double precision
|
||||||
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
|
! replaces "==" but for certain (relative) tolerance. Counterpart to dNeq
|
||||||
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
|
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
logical elemental pure function dEq(a,b,tol)
|
logical elemental pure function dEq(a,b,tol)
|
||||||
|
|
||||||
|
implicit none
|
||||||
real(pReal), intent(in) :: a,b
|
real(pReal), intent(in) :: a,b
|
||||||
real(pReal), intent(in), optional :: tol
|
real(pReal), intent(in), optional :: tol
|
||||||
real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C
|
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
|
||||||
|
|
||||||
dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
|
dEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
|
||||||
end function dEq
|
end function dEq
|
||||||
|
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
!> @brief inequality comparison for double precision
|
!> @brief inequality comparison for float with double precision
|
||||||
! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
|
! replaces "!=" but for certain (relative) tolerance. Counterpart to dEq
|
||||||
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
|
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
logical elemental pure function dNeq(a,b,tol)
|
logical elemental pure function dNeq(a,b,tol)
|
||||||
|
|
||||||
|
implicit none
|
||||||
real(pReal), intent(in) :: a,b
|
real(pReal), intent(in) :: a,b
|
||||||
real(pReal), intent(in), optional :: tol
|
real(pReal), intent(in), optional :: tol
|
||||||
real(pReal), parameter :: eps = 2.2204460492503131E-16 ! DBL_EPSILON in C
|
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
|
||||||
|
|
||||||
dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
|
dNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
|
||||||
end function dNeq
|
end function dNeq
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
!> @brief equality comparison for complex with double precision
|
||||||
|
! replaces "==" but for certain (relative) tolerance. Counterpart to cNeq
|
||||||
|
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
|
||||||
|
! probably a component wise comparison would be more accurate than the comparsion of the absolute
|
||||||
|
! value
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
logical elemental pure function cEq(a,b,tol)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
complex(pReal), intent(in) :: a,b
|
||||||
|
real(pReal), intent(in), optional :: tol
|
||||||
|
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
|
||||||
|
|
||||||
|
cEq = merge(.True., .False.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
|
||||||
|
end function cEq
|
||||||
|
|
||||||
|
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
!> @brief inequality comparison for complex with double precision
|
||||||
|
! replaces "!=" but for certain (relative) tolerance. Counterpart to cEq
|
||||||
|
! http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
|
||||||
|
! probably a component wise comparison would be more accurate than the comparsion of the absolute
|
||||||
|
! value
|
||||||
|
!--------------------------------------------------------------------------------------------------
|
||||||
|
logical elemental pure function cNeq(a,b,tol)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
complex(pReal), intent(in) :: a,b
|
||||||
|
real(pReal), intent(in), optional :: tol
|
||||||
|
real(pReal), parameter :: eps = 2.220446049250313E-16 ! DBL_EPSILON in C
|
||||||
|
|
||||||
|
cNeq = merge(.False., .True.,abs(a-b) <= merge(tol,eps,present(tol))*maxval(abs([a,b])))
|
||||||
|
end function cNeq
|
||||||
|
|
||||||
end module prec
|
end module prec
|
||||||
|
|
|
@ -102,8 +102,6 @@ module spectral_utilities
|
||||||
real(pReal) :: density
|
real(pReal) :: density
|
||||||
end type tSolutionParams
|
end type tSolutionParams
|
||||||
|
|
||||||
type(tSolutionParams), private :: params
|
|
||||||
|
|
||||||
type, public :: phaseFieldDataBin !< set of parameters defining a phase field
|
type, public :: phaseFieldDataBin !< set of parameters defining a phase field
|
||||||
real(pReal) :: diffusion = 0.0_pReal, & !< thermal conductivity
|
real(pReal) :: diffusion = 0.0_pReal, & !< thermal conductivity
|
||||||
mobility = 0.0_pReal, & !< thermal mobility
|
mobility = 0.0_pReal, & !< thermal mobility
|
||||||
|
@ -265,8 +263,9 @@ subroutine utilities_init()
|
||||||
enddo
|
enddo
|
||||||
elseif (divergence_correction == 2_pInt) then
|
elseif (divergence_correction == 2_pInt) then
|
||||||
do j = 1_pInt, 3_pInt
|
do j = 1_pInt, 3_pInt
|
||||||
if (j /= minloc(geomSize/grid,1) .and. j /= maxloc(geomSize/grid,1)) &
|
if ( j /= int(minloc(geomSize/real(grid,pReal),1),pInt) &
|
||||||
scaledGeomSize = geomSize/geomSize(j)*grid(j)
|
.and. j /= int(maxloc(geomSize/real(grid,pReal),1),pInt)) &
|
||||||
|
scaledGeomSize = geomSize/geomSize(j)*real(grid(j),pReal)
|
||||||
enddo
|
enddo
|
||||||
else
|
else
|
||||||
scaledGeomSize = geomSize
|
scaledGeomSize = geomSize
|
||||||
|
@ -403,7 +402,7 @@ subroutine utilities_updateGamma(C,saveReference)
|
||||||
integer(pInt) :: &
|
integer(pInt) :: &
|
||||||
i, j, k, &
|
i, j, k, &
|
||||||
l, m, n, o
|
l, m, n, o
|
||||||
logical :: ierr
|
logical :: err
|
||||||
|
|
||||||
C_ref = C
|
C_ref = C
|
||||||
if (saveReference) then
|
if (saveReference) then
|
||||||
|
@ -427,7 +426,7 @@ subroutine utilities_updateGamma(C,saveReference)
|
||||||
matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex)
|
matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex)
|
||||||
matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex)
|
matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex)
|
||||||
if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then
|
if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then
|
||||||
call math_invert(6_pInt, matA, matInvA, ierr)
|
call math_invert(6_pInt, matA, matInvA, err)
|
||||||
temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal)
|
temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal)
|
||||||
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||||
gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* &
|
gamma_hat(l,m,n,o,i,j,k-grid3Offset) = temp33_complex(l,n)* &
|
||||||
|
@ -543,7 +542,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
|
||||||
integer(pInt) :: &
|
integer(pInt) :: &
|
||||||
i, j, k, &
|
i, j, k, &
|
||||||
l, m, n, o
|
l, m, n, o
|
||||||
logical :: ierr
|
logical :: err
|
||||||
|
|
||||||
|
|
||||||
if (worldrank == 0_pInt) then
|
if (worldrank == 0_pInt) then
|
||||||
|
@ -563,7 +562,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim)
|
||||||
matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex)
|
matA(1:3,1:3) = real(temp33_complex); matA(4:6,4:6) = real(temp33_complex)
|
||||||
matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex)
|
matA(1:3,4:6) = aimag(temp33_complex); matA(4:6,1:3) = -aimag(temp33_complex)
|
||||||
if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then
|
if (abs(math_det33(matA(1:3,1:3))) > 1e-16) then
|
||||||
call math_invert(6_pInt, matA, matInvA, ierr)
|
call math_invert(6_pInt, matA, matInvA, err)
|
||||||
temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal)
|
temp33_complex = cmplx(matInvA(1:3,1:3),matInvA(1:3,4:6),pReal)
|
||||||
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt) &
|
||||||
gamma_hat(l,m,n,o,1,1,1) = temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k)
|
gamma_hat(l,m,n,o,1,1,1) = temp33_complex(l,n)*conjg(-xi1st(o,i,j,k))*xi1st(m,i,j,k)
|
||||||
|
@ -623,6 +622,8 @@ end subroutine utilities_fourierGreenConvolution
|
||||||
!> @brief calculate root mean square of divergence of field_fourier
|
!> @brief calculate root mean square of divergence of field_fourier
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
real(pReal) function utilities_divergenceRMS()
|
real(pReal) function utilities_divergenceRMS()
|
||||||
|
use IO, only: &
|
||||||
|
IO_error
|
||||||
use numerics, only: &
|
use numerics, only: &
|
||||||
worldrank
|
worldrank
|
||||||
use mesh, only: &
|
use mesh, only: &
|
||||||
|
@ -631,8 +632,8 @@ real(pReal) function utilities_divergenceRMS()
|
||||||
grid3
|
grid3
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
integer(pInt) :: i, j, k
|
integer(pInt) :: i, j, k, ierr
|
||||||
PetscErrorCode :: ierr
|
complex(pReal), dimension(3) :: rescaledGeom
|
||||||
|
|
||||||
external :: &
|
external :: &
|
||||||
MPI_Allreduce
|
MPI_Allreduce
|
||||||
|
@ -641,6 +642,7 @@ real(pReal) function utilities_divergenceRMS()
|
||||||
write(6,'(/,a)') ' ... calculating divergence ................................................'
|
write(6,'(/,a)') ' ... calculating divergence ................................................'
|
||||||
flush(6)
|
flush(6)
|
||||||
endif
|
endif
|
||||||
|
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! calculating RMS divergence criterion in Fourier space
|
! calculating RMS divergence criterion in Fourier space
|
||||||
|
@ -649,22 +651,23 @@ real(pReal) function utilities_divergenceRMS()
|
||||||
do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
|
do i = 2_pInt, grid1Red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
|
||||||
utilities_divergenceRMS = utilities_divergenceRMS &
|
utilities_divergenceRMS = utilities_divergenceRMS &
|
||||||
+ 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
|
+ 2.0_pReal*(sum (real(matmul(tensorField_fourier(1:3,1:3,i,j,k),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
|
||||||
conjg(-xi1st(1:3,i,j,k))*geomSize/scaledGeomSize))**2.0_pReal)& ! --> sum squared L_2 norm of vector
|
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal)& ! --> sum squared L_2 norm of vector
|
||||||
+sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),&
|
+sum(aimag(matmul(tensorField_fourier(1:3,1:3,i,j,k),&
|
||||||
conjg(-xi1st(1:3,i,j,k))*geomSize/scaledGeomSize))**2.0_pReal))
|
conjg(-xi1st(1:3,i,j,k))*rescaledGeom))**2.0_pReal))
|
||||||
enddo
|
enddo
|
||||||
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1)
|
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if grid(1) /= 1)
|
||||||
+ sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
|
+ sum( real(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
|
||||||
conjg(-xi1st(1:3,1,j,k))*geomSize/scaledGeomSize))**2.0_pReal) &
|
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) &
|
||||||
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
|
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,1 ,j,k), &
|
||||||
conjg(-xi1st(1:3,1,j,k))*geomSize/scaledGeomSize))**2.0_pReal) &
|
conjg(-xi1st(1:3,1,j,k))*rescaledGeom))**2.0_pReal) &
|
||||||
+ sum( real(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
|
+ sum( real(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
|
||||||
conjg(-xi1st(1:3,grid1Red,j,k))*geomSize/scaledGeomSize))**2.0_pReal) &
|
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal) &
|
||||||
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
|
+ sum(aimag(matmul(tensorField_fourier(1:3,1:3,grid1Red,j,k), &
|
||||||
conjg(-xi1st(1:3,grid1Red,j,k))*geomSize/scaledGeomSize))**2.0_pReal)
|
conjg(-xi1st(1:3,grid1Red,j,k))*rescaledGeom))**2.0_pReal)
|
||||||
enddo; enddo
|
enddo; enddo
|
||||||
if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
|
if(grid(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
|
||||||
call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
|
call MPI_Allreduce(MPI_IN_PLACE,utilities_divergenceRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
|
||||||
|
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_divergenceRMS')
|
||||||
utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
|
utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
|
||||||
|
|
||||||
|
|
||||||
|
@ -675,6 +678,8 @@ end function utilities_divergenceRMS
|
||||||
!> @brief calculate max of curl of field_fourier
|
!> @brief calculate max of curl of field_fourier
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
real(pReal) function utilities_curlRMS()
|
real(pReal) function utilities_curlRMS()
|
||||||
|
use IO, only: &
|
||||||
|
IO_error
|
||||||
use numerics, only: &
|
use numerics, only: &
|
||||||
worldrank
|
worldrank
|
||||||
use mesh, only: &
|
use mesh, only: &
|
||||||
|
@ -683,9 +688,9 @@ real(pReal) function utilities_curlRMS()
|
||||||
grid3
|
grid3
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
integer(pInt) :: i, j, k, l
|
integer(pInt) :: i, j, k, l, ierr
|
||||||
complex(pReal), dimension(3,3) :: curl_fourier
|
complex(pReal), dimension(3,3) :: curl_fourier
|
||||||
PetscErrorCode :: ierr
|
complex(pReal), dimension(3) :: rescaledGeom
|
||||||
|
|
||||||
external :: &
|
external :: &
|
||||||
MPI_Reduce, &
|
MPI_Reduce, &
|
||||||
|
@ -695,6 +700,7 @@ real(pReal) function utilities_curlRMS()
|
||||||
write(6,'(/,a)') ' ... calculating curl ......................................................'
|
write(6,'(/,a)') ' ... calculating curl ......................................................'
|
||||||
flush(6)
|
flush(6)
|
||||||
endif
|
endif
|
||||||
|
rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal)
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! calculating max curl criterion in Fourier space
|
! calculating max curl criterion in Fourier space
|
||||||
|
@ -703,39 +709,40 @@ real(pReal) function utilities_curlRMS()
|
||||||
do k = 1_pInt, grid3; do j = 1_pInt, grid(2);
|
do k = 1_pInt, grid3; do j = 1_pInt, grid(2);
|
||||||
do i = 2_pInt, grid1Red - 1_pInt
|
do i = 2_pInt, grid1Red - 1_pInt
|
||||||
do l = 1_pInt, 3_pInt
|
do l = 1_pInt, 3_pInt
|
||||||
curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*geomSize(2)/scaledGeomSize(2) &
|
curl_fourier(l,1) = (+tensorField_fourier(l,3,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2) &
|
||||||
-tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*geomSize(3)/scaledGeomSize(3))
|
-tensorField_fourier(l,2,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3))
|
||||||
curl_fourier(l,2) = (+tensorField_fourier(l,1,i,j,k)*xi1st(3,i,j,k)*geomSize(3)/scaledGeomSize(3) &
|
curl_fourier(l,2) = (+tensorField_fourier(l,1,i,j,k)*xi1st(3,i,j,k)*rescaledGeom(3) &
|
||||||
-tensorField_fourier(l,3,i,j,k)*xi1st(1,i,j,k)*geomSize(1)/scaledGeomSize(1))
|
-tensorField_fourier(l,3,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1))
|
||||||
curl_fourier(l,3) = (+tensorField_fourier(l,2,i,j,k)*xi1st(1,i,j,k)*geomSize(1)/scaledGeomSize(1) &
|
curl_fourier(l,3) = (+tensorField_fourier(l,2,i,j,k)*xi1st(1,i,j,k)*rescaledGeom(1) &
|
||||||
-tensorField_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*geomSize(2)/scaledGeomSize(2))
|
-tensorField_fourier(l,1,i,j,k)*xi1st(2,i,j,k)*rescaledGeom(2))
|
||||||
enddo
|
enddo
|
||||||
utilities_curlRMS = utilities_curlRMS + &
|
utilities_curlRMS = utilities_curlRMS + &
|
||||||
2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! Has somewhere a conj. complex counterpart. Therefore count it twice.
|
2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! Has somewhere a conj. complex counterpart. Therefore count it twice.
|
||||||
enddo
|
enddo
|
||||||
do l = 1_pInt, 3_pInt
|
do l = 1_pInt, 3_pInt
|
||||||
curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*geomSize(2)/scaledGeomSize(2) &
|
curl_fourier = (+tensorField_fourier(l,3,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2) &
|
||||||
-tensorField_fourier(l,2,1,j,k)*xi1st(3,1,j,k)*geomSize(3)/scaledGeomSize(3))
|
-tensorField_fourier(l,2,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3))
|
||||||
curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi1st(3,1,j,k)*geomSize(3)/scaledGeomSize(3) &
|
curl_fourier = (+tensorField_fourier(l,1,1,j,k)*xi1st(3,1,j,k)*rescaledGeom(3) &
|
||||||
-tensorField_fourier(l,3,1,j,k)*xi1st(1,1,j,k)*geomSize(1)/scaledGeomSize(1))
|
-tensorField_fourier(l,3,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1))
|
||||||
curl_fourier = (+tensorField_fourier(l,2,1,j,k)*xi1st(1,1,j,k)*geomSize(1)/scaledGeomSize(1) &
|
curl_fourier = (+tensorField_fourier(l,2,1,j,k)*xi1st(1,1,j,k)*rescaledGeom(1) &
|
||||||
-tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*geomSize(2)/scaledGeomSize(2))
|
-tensorField_fourier(l,1,1,j,k)*xi1st(2,1,j,k)*rescaledGeom(2))
|
||||||
enddo
|
enddo
|
||||||
utilities_curlRMS = utilities_curlRMS + &
|
utilities_curlRMS = utilities_curlRMS + &
|
||||||
sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1)
|
sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (DC) does not have a conjugate complex counterpart (if grid(1) /= 1)
|
||||||
do l = 1_pInt, 3_pInt
|
do l = 1_pInt, 3_pInt
|
||||||
curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*geomSize(2)/scaledGeomSize(2) &
|
curl_fourier = (+tensorField_fourier(l,3,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2) &
|
||||||
-tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*geomSize(3)/scaledGeomSize(3))
|
-tensorField_fourier(l,2,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3))
|
||||||
curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*geomSize(3)/scaledGeomSize(3) &
|
curl_fourier = (+tensorField_fourier(l,1,grid1Red,j,k)*xi1st(3,grid1Red,j,k)*rescaledGeom(3) &
|
||||||
-tensorField_fourier(l,3,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*geomSize(1)/scaledGeomSize(1))
|
-tensorField_fourier(l,3,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1))
|
||||||
curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*geomSize(1)/scaledGeomSize(1) &
|
curl_fourier = (+tensorField_fourier(l,2,grid1Red,j,k)*xi1st(1,grid1Red,j,k)*rescaledGeom(1) &
|
||||||
-tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*geomSize(2)/scaledGeomSize(2))
|
-tensorField_fourier(l,1,grid1Red,j,k)*xi1st(2,grid1Red,j,k)*rescaledGeom(2))
|
||||||
enddo
|
enddo
|
||||||
utilities_curlRMS = utilities_curlRMS + &
|
utilities_curlRMS = utilities_curlRMS + &
|
||||||
sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1)
|
sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)! this layer (Nyquist) does not have a conjugate complex counterpart (if grid(1) /= 1)
|
||||||
enddo; enddo
|
enddo; enddo
|
||||||
|
|
||||||
call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
|
call MPI_Allreduce(MPI_IN_PLACE,utilities_curlRMS,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
|
||||||
|
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='utilities_curlRMS')
|
||||||
utilities_curlRMS = sqrt(utilities_curlRMS) * wgt
|
utilities_curlRMS = sqrt(utilities_curlRMS) * wgt
|
||||||
if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
|
if(grid(1) == 1_pInt) utilities_curlRMS = utilities_curlRMS * 0.5_pReal ! counted twice in case of grid(1) == 1
|
||||||
|
|
||||||
|
@ -931,6 +938,10 @@ end subroutine utilities_fourierTensorDivergence
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
||||||
P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC)
|
P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC)
|
||||||
|
use prec, only: &
|
||||||
|
dNeq
|
||||||
|
use IO, only: &
|
||||||
|
IO_error
|
||||||
use debug, only: &
|
use debug, only: &
|
||||||
debug_reset, &
|
debug_reset, &
|
||||||
debug_info
|
debug_info
|
||||||
|
@ -969,10 +980,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
||||||
age
|
age
|
||||||
|
|
||||||
integer(pInt) :: &
|
integer(pInt) :: &
|
||||||
j,k
|
j,k,ierr
|
||||||
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
|
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
|
||||||
real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
|
real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
|
||||||
PetscErrorCode :: ierr
|
|
||||||
|
|
||||||
external :: &
|
external :: &
|
||||||
MPI_Reduce, &
|
MPI_Reduce, &
|
||||||
|
@ -1006,7 +1016,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
||||||
defgradDetMin = min(defgradDetMin,defgradDet)
|
defgradDetMin = min(defgradDetMin,defgradDet)
|
||||||
end do
|
end do
|
||||||
call MPI_reduce(MPI_IN_PLACE,defgradDetMax,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr)
|
call MPI_reduce(MPI_IN_PLACE,defgradDetMax,1,MPI_DOUBLE,MPI_MAX,0,PETSC_COMM_WORLD,ierr)
|
||||||
|
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max')
|
||||||
call MPI_reduce(MPI_IN_PLACE,defgradDetMin,1,MPI_DOUBLE,MPI_MIN,0,PETSC_COMM_WORLD,ierr)
|
call MPI_reduce(MPI_IN_PLACE,defgradDetMin,1,MPI_DOUBLE,MPI_MIN,0,PETSC_COMM_WORLD,ierr)
|
||||||
|
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min')
|
||||||
if (worldrank == 0_pInt) then
|
if (worldrank == 0_pInt) then
|
||||||
write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
|
write(6,'(a,1x,es11.4)') ' max determinant of deformation =', defgradDetMax
|
||||||
write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
|
write(6,'(a,1x,es11.4)') ' min determinant of deformation =', defgradDetMin
|
||||||
|
@ -1032,7 +1044,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, &
|
||||||
end do
|
end do
|
||||||
|
|
||||||
call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
|
call MPI_Allreduce(MPI_IN_PLACE,max_dPdF,81,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
|
||||||
|
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce max')
|
||||||
call MPI_Allreduce(MPI_IN_PLACE,min_dPdF,81,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD,ierr)
|
call MPI_Allreduce(MPI_IN_PLACE,min_dPdF,81,MPI_DOUBLE,MPI_MIN,PETSC_COMM_WORLD,ierr)
|
||||||
|
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='MPI_Allreduce min')
|
||||||
|
|
||||||
C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF)
|
C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF)
|
||||||
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
|
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
|
||||||
|
@ -1187,6 +1201,10 @@ end function utilities_getFreqDerivative
|
||||||
! convolution
|
! convolution
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
subroutine utilities_updateIPcoords(F)
|
subroutine utilities_updateIPcoords(F)
|
||||||
|
use prec, only: &
|
||||||
|
cNeq
|
||||||
|
use IO, only: &
|
||||||
|
IO_error
|
||||||
use math, only: &
|
use math, only: &
|
||||||
math_mul33x3
|
math_mul33x3
|
||||||
use mesh, only: &
|
use mesh, only: &
|
||||||
|
@ -1198,10 +1216,9 @@ subroutine utilities_updateIPcoords(F)
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
|
real(pReal), dimension(3,3,grid(1),grid(2),grid3), intent(in) :: F
|
||||||
integer(pInt) :: i, j, k, m
|
integer(pInt) :: i, j, k, m, ierr
|
||||||
real(pReal), dimension(3) :: step, offset_coords
|
real(pReal), dimension(3) :: step, offset_coords
|
||||||
real(pReal), dimension(3,3) :: Favg
|
real(pReal), dimension(3,3) :: Favg
|
||||||
PetscErrorCode :: ierr
|
|
||||||
external &
|
external &
|
||||||
MPI_Bcast
|
MPI_Bcast
|
||||||
|
|
||||||
|
@ -1213,7 +1230,7 @@ subroutine utilities_updateIPcoords(F)
|
||||||
call utilities_fourierTensorDivergence()
|
call utilities_fourierTensorDivergence()
|
||||||
|
|
||||||
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
|
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid1Red
|
||||||
if (any(abs(xi1st(1:3,i,j,k)) > tiny(0.0_pReal))) &
|
if (any(cNeq(xi1st(1:3,i,j,k),cmplx(0.0_pReal,0.0_pReal)))) &
|
||||||
vectorField_fourier(1:3,i,j,k) = vectorField_fourier(1:3,i,j,k)/ &
|
vectorField_fourier(1:3,i,j,k) = vectorField_fourier(1:3,i,j,k)/ &
|
||||||
sum(conjg(-xi1st(1:3,i,j,k))*xi1st(1:3,i,j,k))
|
sum(conjg(-xi1st(1:3,i,j,k))*xi1st(1:3,i,j,k))
|
||||||
enddo; enddo; enddo
|
enddo; enddo; enddo
|
||||||
|
@ -1223,12 +1240,14 @@ subroutine utilities_updateIPcoords(F)
|
||||||
! average F
|
! average F
|
||||||
if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
|
if (grid3Offset == 0_pInt) Favg = real(tensorField_fourier(1:3,1:3,1,1,1),pReal)*wgt
|
||||||
call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
|
call MPI_Bcast(Favg,9,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
|
||||||
|
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords')
|
||||||
|
|
||||||
!--------------------------------------------------------------------------------------------------
|
!--------------------------------------------------------------------------------------------------
|
||||||
! add average to fluctuation and put (0,0,0) on (0,0,0)
|
! add average to fluctuation and put (0,0,0) on (0,0,0)
|
||||||
step = geomSize/real(grid, pReal)
|
step = geomSize/real(grid, pReal)
|
||||||
if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1)
|
if (grid3Offset == 0_pInt) offset_coords = vectorField_real(1:3,1,1,1)
|
||||||
call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
|
call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
|
||||||
|
if(ierr /=0_pInt) call IO_error(894_pInt, ext_msg='update_IPcoords')
|
||||||
offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords
|
offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords
|
||||||
m = 1_pInt
|
m = 1_pInt
|
||||||
do k = 1_pInt,grid3; do j = 1_pInt,grid(2); do i = 1_pInt,grid(1)
|
do k = 1_pInt,grid3; do j = 1_pInt,grid(2); do i = 1_pInt,grid(1)
|
||||||
|
|
|
@ -236,7 +236,7 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
|
||||||
use material, only: &
|
use material, only: &
|
||||||
homogenization_Ngrains, &
|
homogenization_Ngrains, &
|
||||||
mappingHomogenization, &
|
mappingHomogenization, &
|
||||||
phaseAt, phasememberAt, &
|
phaseAt, &
|
||||||
thermal_typeInstance, &
|
thermal_typeInstance, &
|
||||||
phase_Nsources, &
|
phase_Nsources, &
|
||||||
phase_source, &
|
phase_source, &
|
||||||
|
@ -297,8 +297,8 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
Tdot = Tdot/homogenization_Ngrains(homog)
|
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal)
|
||||||
dTdot_dT = dTdot_dT/homogenization_Ngrains(homog)
|
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal)
|
||||||
|
|
||||||
end subroutine thermal_adiabatic_getSourceAndItsTangent
|
end subroutine thermal_adiabatic_getSourceAndItsTangent
|
||||||
|
|
||||||
|
@ -336,8 +336,7 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
thermal_adiabatic_getSpecificHeat = &
|
thermal_adiabatic_getSpecificHeat = &
|
||||||
thermal_adiabatic_getSpecificHeat/ &
|
thermal_adiabatic_getSpecificHeat/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function thermal_adiabatic_getSpecificHeat
|
end function thermal_adiabatic_getSpecificHeat
|
||||||
|
|
||||||
|
@ -375,8 +374,7 @@ function thermal_adiabatic_getMassDensity(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
thermal_adiabatic_getMassDensity = &
|
thermal_adiabatic_getMassDensity = &
|
||||||
thermal_adiabatic_getMassDensity/ &
|
thermal_adiabatic_getMassDensity/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function thermal_adiabatic_getMassDensity
|
end function thermal_adiabatic_getMassDensity
|
||||||
|
|
||||||
|
|
|
@ -190,7 +190,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
|
||||||
use material, only: &
|
use material, only: &
|
||||||
homogenization_Ngrains, &
|
homogenization_Ngrains, &
|
||||||
mappingHomogenization, &
|
mappingHomogenization, &
|
||||||
phaseAt, phasememberAt, &
|
phaseAt, &
|
||||||
thermal_typeInstance, &
|
thermal_typeInstance, &
|
||||||
phase_Nsources, &
|
phase_Nsources, &
|
||||||
phase_source, &
|
phase_source, &
|
||||||
|
@ -252,8 +252,8 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
Tdot = Tdot/homogenization_Ngrains(homog)
|
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal)
|
||||||
dTdot_dT = dTdot_dT/homogenization_Ngrains(homog)
|
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal)
|
||||||
|
|
||||||
end subroutine thermal_conduction_getSourceAndItsTangent
|
end subroutine thermal_conduction_getSourceAndItsTangent
|
||||||
|
|
||||||
|
@ -291,8 +291,7 @@ function thermal_conduction_getConductivity33(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
thermal_conduction_getConductivity33 = &
|
thermal_conduction_getConductivity33 = &
|
||||||
thermal_conduction_getConductivity33/ &
|
thermal_conduction_getConductivity33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function thermal_conduction_getConductivity33
|
end function thermal_conduction_getConductivity33
|
||||||
|
|
||||||
|
@ -330,8 +329,7 @@ function thermal_conduction_getSpecificHeat(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
thermal_conduction_getSpecificHeat = &
|
thermal_conduction_getSpecificHeat = &
|
||||||
thermal_conduction_getSpecificHeat/ &
|
thermal_conduction_getSpecificHeat/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function thermal_conduction_getSpecificHeat
|
end function thermal_conduction_getSpecificHeat
|
||||||
|
|
||||||
|
@ -369,8 +367,7 @@ function thermal_conduction_getMassDensity(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
thermal_conduction_getMassDensity = &
|
thermal_conduction_getMassDensity = &
|
||||||
thermal_conduction_getMassDensity/ &
|
thermal_conduction_getMassDensity/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function thermal_conduction_getMassDensity
|
end function thermal_conduction_getMassDensity
|
||||||
|
|
||||||
|
|
|
@ -219,7 +219,7 @@ subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv
|
||||||
use material, only: &
|
use material, only: &
|
||||||
homogenization_Ngrains, &
|
homogenization_Ngrains, &
|
||||||
mappingHomogenization, &
|
mappingHomogenization, &
|
||||||
phaseAt, phasememberAt, &
|
phaseAt, &
|
||||||
phase_source, &
|
phase_source, &
|
||||||
phase_Nsources, &
|
phase_Nsources, &
|
||||||
SOURCE_vacancy_phenoplasticity_ID, &
|
SOURCE_vacancy_phenoplasticity_ID, &
|
||||||
|
@ -266,8 +266,8 @@ subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
CvDot = CvDot/homogenization_Ngrains(mappingHomogenization(2,ip,el))
|
CvDot = CvDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
|
||||||
dCvDot_dCv = dCvDot_dCv/homogenization_Ngrains(mappingHomogenization(2,ip,el))
|
dCvDot_dCv = dCvDot_dCv/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
|
||||||
|
|
||||||
end subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent
|
end subroutine vacancyflux_cahnhilliard_getSourceAndItsTangent
|
||||||
|
|
||||||
|
@ -301,8 +301,7 @@ function vacancyflux_cahnhilliard_getMobility33(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
vacancyflux_cahnhilliard_getMobility33 = &
|
vacancyflux_cahnhilliard_getMobility33 = &
|
||||||
vacancyflux_cahnhilliard_getMobility33/ &
|
vacancyflux_cahnhilliard_getMobility33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function vacancyflux_cahnhilliard_getMobility33
|
end function vacancyflux_cahnhilliard_getMobility33
|
||||||
|
|
||||||
|
@ -336,8 +335,7 @@ function vacancyflux_cahnhilliard_getDiffusion33(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
vacancyflux_cahnhilliard_getDiffusion33 = &
|
vacancyflux_cahnhilliard_getDiffusion33 = &
|
||||||
vacancyflux_cahnhilliard_getDiffusion33/ &
|
vacancyflux_cahnhilliard_getDiffusion33/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function vacancyflux_cahnhilliard_getDiffusion33
|
end function vacancyflux_cahnhilliard_getDiffusion33
|
||||||
|
|
||||||
|
@ -371,8 +369,7 @@ real(pReal) function vacancyflux_cahnhilliard_getFormationEnergy(ip,el)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
vacancyflux_cahnhilliard_getFormationEnergy = &
|
vacancyflux_cahnhilliard_getFormationEnergy = &
|
||||||
vacancyflux_cahnhilliard_getFormationEnergy/ &
|
vacancyflux_cahnhilliard_getFormationEnergy/real(homogenization_Ngrains(mesh_element(3,el)),pReal)
|
||||||
homogenization_Ngrains(mesh_element(3,el))
|
|
||||||
|
|
||||||
end function vacancyflux_cahnhilliard_getFormationEnergy
|
end function vacancyflux_cahnhilliard_getFormationEnergy
|
||||||
|
|
||||||
|
@ -408,7 +405,7 @@ real(pReal) function vacancyflux_cahnhilliard_getEntropicCoeff(ip,el)
|
||||||
vacancyflux_cahnhilliard_getEntropicCoeff = &
|
vacancyflux_cahnhilliard_getEntropicCoeff = &
|
||||||
vacancyflux_cahnhilliard_getEntropicCoeff* &
|
vacancyflux_cahnhilliard_getEntropicCoeff* &
|
||||||
temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))/ &
|
temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el))/ &
|
||||||
homogenization_Ngrains(material_homog(ip,el))
|
real(homogenization_Ngrains(material_homog(ip,el)),pReal)
|
||||||
|
|
||||||
end function vacancyflux_cahnhilliard_getEntropicCoeff
|
end function vacancyflux_cahnhilliard_getEntropicCoeff
|
||||||
|
|
||||||
|
@ -467,8 +464,8 @@ subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent(KPot, dKPot_dC
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
KPot = KPot/homogenization_Ngrains(material_homog(ip,el))
|
KPot = KPot/real(homogenization_Ngrains(material_homog(ip,el)),pReal)
|
||||||
dKPot_dCv = dKPot_dCv/homogenization_Ngrains(material_homog(ip,el))
|
dKPot_dCv = dKPot_dCv/real(homogenization_Ngrains(material_homog(ip,el)),pReal)
|
||||||
|
|
||||||
end subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent
|
end subroutine vacancyflux_cahnhilliard_KinematicChemPotAndItsTangent
|
||||||
|
|
||||||
|
|
|
@ -235,7 +235,7 @@ subroutine vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv,
|
||||||
use material, only: &
|
use material, only: &
|
||||||
homogenization_Ngrains, &
|
homogenization_Ngrains, &
|
||||||
mappingHomogenization, &
|
mappingHomogenization, &
|
||||||
phaseAt, phasememberAt, &
|
phaseAt, &
|
||||||
phase_source, &
|
phase_source, &
|
||||||
phase_Nsources, &
|
phase_Nsources, &
|
||||||
SOURCE_vacancy_phenoplasticity_ID, &
|
SOURCE_vacancy_phenoplasticity_ID, &
|
||||||
|
@ -282,8 +282,8 @@ subroutine vacancyflux_isochempot_getSourceAndItsTangent(CvDot, dCvDot_dCv, Cv,
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
CvDot = CvDot/homogenization_Ngrains(mappingHomogenization(2,ip,el))
|
CvDot = CvDot/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
|
||||||
dCvDot_dCv = dCvDot_dCv/homogenization_Ngrains(mappingHomogenization(2,ip,el))
|
dCvDot_dCv = dCvDot_dCv/real(homogenization_Ngrains(mappingHomogenization(2,ip,el)),pReal)
|
||||||
|
|
||||||
end subroutine vacancyflux_isochempot_getSourceAndItsTangent
|
end subroutine vacancyflux_isochempot_getSourceAndItsTangent
|
||||||
|
|
||||||
|
|
Loading…
Reference in New Issue