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

Merge branch 'no-crystallite-dPdF' into 'development' 

No crystallite dPdF

See merge request damask/DAMASK!246
This commit is contained in:
Franz Roters 2020-10-05 15:20:44 +02:00
parent c9ead8bf59 88b92bf3b1
commit 37dfd32f5e
7 changed files with 68 additions and 82 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 3b498f5cb3c50e669588106de1b4cdc4c03ffff1
Subproject commit 64e62f805b5ad784e3397ee5f735aaeb3cc134c2

63
examples/ConfigFiles/numerics.yaml
View File

@ -3,27 +3,27 @@
homogenization:
mech:
RGC:
atol: 1.0e+4 # absolute tolerance of RGC residuum (in Pa)
rtol: 1.0e-3 # relative ...
amax: 1.0e+10 # absolute upper-limit of RGC residuum (in Pa)
rmax: 1.0e+2 # relative ...
perturbpenalty: 1.0e-7 # perturbation for computing penalty tangent
relevantmismatch: 1.0e-5 # minimum threshold of mismatch
viscositypower: 1.0e+0 # power (sensitivity rate) of numerical viscosity in RGC scheme
viscositymodulus: 0.0e+0 # stress modulus of RGC numerical viscosity (zero = without numerical viscosity)
# suggestion: larger than the aTol_RGC but still far below the expected flow stress of material
refrelaxationrate: 1.0e-3 # reference rate of relaxation (about the same magnitude as straining rate, possibly a bit higher)
maxrelaxationrate: 1.0e+0 # threshold of maximum relaxation vector increment (if exceed this then cutback)
maxvoldiscrepancy: 1.0e-5 # maximum allowable relative volume discrepancy
voldiscrepancymod: 1.0e+12
discrepancypower: 5.0
RGC:
atol: 1.0e+4 # absolute tolerance of RGC residuum (in Pa)
rtol: 1.0e-3 # relative ...
amax: 1.0e+10 # absolute upper-limit of RGC residuum (in Pa)
rmax: 1.0e+2 # relative ...
perturbpenalty: 1.0e-7 # perturbation for computing penalty tangent
relevantmismatch: 1.0e-5 # minimum threshold of mismatch
viscositypower: 1.0e+0 # power (sensitivity rate) of numerical viscosity in RGC scheme
viscositymodulus: 0.0e+0 # stress modulus of RGC numerical viscosity (zero = without numerical viscosity)
# suggestion: larger than the aTol_RGC but still far below the expected flow stress of material
refrelaxationrate: 1.0e-3 # reference rate of relaxation (about the same magnitude as straining rate, possibly a bit higher)
maxrelaxationrate: 1.0e+0 # threshold of maximum relaxation vector increment (if exceed this then cutback)
maxvoldiscrepancy: 1.0e-5 # maximum allowable relative volume discrepancy
voldiscrepancymod: 1.0e+12
discrepancypower: 5.0
generic:
subStepMin: 1.0e-3 # minimum (relative) size of sub-step allowed during cutback in homogenization
subStepSize: 0.25 # size of substep when cutback introduced in homogenization (value between 0 and 1)
stepIncrease: 1.5 # increase of next substep size when previous substep converged in homogenization (value higher than 1)
nMPstate: 10 # materialpoint state loop limit
subStepMin: 1.0e-3 # minimum (relative) size of sub-step allowed during cutback in homogenization
subStepSize: 0.25 # size of substep when cutback introduced in homogenization (value between 0 and 1)
stepIncrease: 1.5 # increase of next substep size when previous substep converged in homogenization (value higher than 1)
nMPstate: 10 # materialpoint state loop limit
grid:
eps_div_atol: 1.0e-3 # absolute tolerance for fulfillment of stress equilibrium
@ -39,30 +39,30 @@ grid:
itmax: 250 # Maximum iteration number
itmin: 2 # Minimum iteration number
fftw_timelimit: -1.0 # timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit
fftw_plan_mode: FFTW_PATIENT # reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag
maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
maxStaggeredIter: 10 # max number of field level staggered iterations
fftw_plan_mode: FFTW_PATIENT # reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag
maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
maxStaggeredIter: 10 # max number of field level staggered iterations
memory_efficient: 1 # Precalculate Gamma-operator (81 double per point)
update_gamma: false # Update Gamma-operator with current dPdF (not possible if memory_efficient=1)
divergence_correction: 2 # Use size-independent divergence criterion
divergence_correction: 2 # Use size-independent divergence criterion
derivative: continuous # Approximation used for derivatives in Fourier space
solver: Basic # Type of spectral solver (BasicPETSc/Polarisation/FEM)
petsc_options: -snes_type ngmres -snes_ngmres_anderson # PetSc solver options
alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
mesh:
maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
maxStaggeredIter: 10 # max number of field level staggered iterations
maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
maxStaggeredIter: 10 # max number of field level staggered iterations
structorder: 2 # order of displacement shape functions (when mesh is defined)
bbarstabilisation: false
bbarstabilisation: false
integrationorder: 2 # order of quadrature rule required (when mesh is defined)
itmax: 250 # Maximum iteration number
itmin: 2 # Minimum iteration number
eps_struct_atol: 1.0e-10 # absolute tolerance for mechanical equilibrium
eps_struct_rtol: 1.0e-4 # relative tolerance for mechanical equilibrium
crystallite:
subStepMin: 1.0e-3 # minimum (relative) size of sub-step allowed during cutback in crystallite
subStepSize: 0.25 # size of substep when cutback introduced in crystallite (value between 0 and 1)
@ -78,12 +78,9 @@ crystallite:
iJacoLpresiduum: 1 # frequency of Jacobian update of residuum in Lp
commercialFEM:
ijacostiffness: 1 # frequency of stiffness update
unitlength: 1 # physical length of one computational length unit
generic:
charLength: 1.0 # characteristic length scale for gradient problems.
random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed.
residualStiffness: 1.0e-6 # non-zero residual damage.

11
src/CPFEM.f90
View File

@ -103,7 +103,6 @@ end subroutine CPFEM_initAll
subroutine CPFEM_init
class(tNode), pointer :: &
num_commercialFEM, &
debug_CPFEM
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
@ -112,12 +111,6 @@ subroutine CPFEM_init
allocate(CPFEM_dcsdE( 6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
allocate(CPFEM_dcsdE_knownGood(6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
!------------------------------------------------------------------------------
! read numerical parameters and do sanity check
num_commercialFEM => config_numerics%get('commercialFEM',defaultVal=emptyDict)
num%iJacoStiffness = num_commercialFEM%get_asInt('ijacostiffness',defaultVal=1)
if (num%iJacoStiffness < 1) call IO_error(301,ext_msg='iJacoStiffness')
!------------------------------------------------------------------------------
! read debug options
@ -161,7 +154,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog, mySource
logical updateJaco ! flag indicating if Jacobian has to be updated
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll
@ -204,12 +196,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
else validCalculation
updateJaco = mod(cycleCounter,num%iJacoStiffness) == 0
FEsolving_execElem = elCP
FEsolving_execIP = ip
if (debugCPFEM%extensive) &
print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
call materialpoint_stressAndItsTangent(updateJaco, dt)
call materialpoint_stressAndItsTangent(dt)
terminalIllness: if (terminallyIll) then

41
src/crystallite.f90
View File

@ -69,8 +69,6 @@ module crystallite
real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_partionedF !< def grad to be reached at end of homog inc
real(pReal), dimension(:,:,:,:,:,:,:), allocatable, public, protected :: &
crystallite_dPdF !< current individual dPdF per grain (end of converged time step)
logical, dimension(:,:,:), allocatable, public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
logical, dimension(:,:,:), allocatable :: &
@ -183,8 +181,6 @@ subroutine crystallite_init
crystallite_subFp0,crystallite_subFi0, &
source = crystallite_partionedF)
allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_dt(cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_subdt,crystallite_subFrac,crystallite_subStep, &
source = crystallite_dt)
@ -293,7 +289,6 @@ subroutine crystallite_init
!$OMP END PARALLEL DO
devNull = crystallite_stress()
call crystallite_stressTangent
#ifdef DEBUG
if (debugCrystallite%basic) then
@ -566,12 +561,14 @@ end subroutine crystallite_restore
!--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF).
!--------------------------------------------------------------------------------------------------
subroutine crystallite_stressTangent
function crystallite_stressTangent(c,i,e) result(dPdF)
integer :: &
real(pReal), dimension(3,3,3,3) :: dPdF
integer, intent(in) :: &
c, & !< counter in constituent loop
i, & !< counter in integration point loop
e, & !< counter in element loop
e !< counter in element loop
integer :: &
o, &
p
@ -593,12 +590,6 @@ subroutine crystallite_stressTangent
real(pReal), dimension(9,9):: temp_99
logical :: error
!$OMP PARALLEL DO PRIVATE(dSdF,dSdFe,dSdFi,dLpdS,dLpdFi,dFpinvdF,dLidS,dLidFi,dFidS,o,p, &
!$OMP invSubFp0,invSubFi0,invFp,invFi, &
!$OMP rhs_3333,lhs_3333,temp_99,temp_33_1,temp_33_2,temp_33_3,temp_33_4,temp_3333,error)
elementLooping: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e))
call constitutive_SandItsTangents(devNull,dSdFe,dSdFi, &
crystallite_Fe(1:3,1:3,c,i,e), &
@ -679,24 +670,20 @@ subroutine crystallite_stressTangent
temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e),invFp)
temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,c,i,e))
crystallite_dPdF(1:3,1:3,1:3,1:3,c,i,e) = 0.0_pReal
dPdF = 0.0_pReal
do p=1,3
crystallite_dPdF(p,1:3,p,1:3,c,i,e) = transpose(temp_33_2)
dPdF(p,1:3,p,1:3) = transpose(temp_33_2)
enddo
do o=1,3; do p=1,3
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) &
+ matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), &
transpose(invFp)) &
+ matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
+ matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), &
transpose(invFp)) &
+ matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo
enddo; enddo
enddo elementLooping
!$OMP END PARALLEL DO
end subroutine crystallite_stressTangent
end function crystallite_stressTangent
!--------------------------------------------------------------------------------------------------

2
src/grid/spectral_utilities.f90
View File

@ -826,7 +826,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
call materialpoint_stressAndItsTangent(timeinc) ! calculate P field
P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P

29
src/homogenization.f90
View File

@ -204,10 +204,9 @@ end subroutine homogenization_init
!--------------------------------------------------------------------------------------------------
!> @brief parallelized calculation of stress and corresponding tangent at material points
!--------------------------------------------------------------------------------------------------
subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
subroutine materialpoint_stressAndItsTangent(dt)
real(pReal), intent(in) :: dt !< time increment
logical, intent(in) :: updateJaco !< initiating Jacobian update
integer :: &
NiterationHomog, &
NiterationMPstate, &
@ -375,8 +374,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo cutBackLooping
if(updateJaco) call crystallite_stressTangent
if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO
@ -438,11 +435,16 @@ function updateState(subdt,subF,ip,el)
integer, intent(in) :: &
ip, & !< integration point
el !< element number
integer :: c
logical, dimension(2) :: updateState
real(pReal) :: dPdFs(3,3,3,3,homogenization_Ngrains(material_homogenizationAt(el)))
updateState = .true.
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c=1,homogenization_Ngrains(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo
updateState = &
updateState .and. &
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
@ -450,7 +452,7 @@ function updateState(subdt,subF,ip,el)
crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),&
subF,&
subdt, &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
dPdFs, &
ip, &
el)
end select chosenHomogenization
@ -484,26 +486,35 @@ subroutine averageStressAndItsTangent(ip,el)
integer, intent(in) :: &
ip, & !< integration point
el !< element number
integer :: c
real(pReal) :: dPdFs(3,3,3,3,homogenization_Ngrains(material_homogenizationAt(el)))
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el)
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_dPdF(1:3,1:3,1:3,1:3,1,ip,el)
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do c = 1, homogenization_Ngrains(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo
call mech_isostrain_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c = 1, homogenization_Ngrains(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo
call mech_RGC_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))
end select chosenHomogenization

2
src/mesh/FEM_utilities.f90
View File

@ -160,7 +160,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
print'(/,a)', ' ... evaluating constitutive response ......................................'
call materialpoint_stressAndItsTangent(.true.,timeinc) ! calculate P field
call materialpoint_stressAndItsTangent(timeinc) ! calculate P field
cutBack = .false. ! reset cutBack status