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

do not store dPdF at the crystallite level

This commit is contained in:
Martin Diehl 2020-09-30 10:53:05 +02:00
parent 7f8613f6ad
commit 57174d0aba
3 changed files with 63 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
examples/ConfigFiles/numerics.yaml
View File

@ -3,27 +3,27 @@
homogenization: homogenization:
mech: mech:
RGC: RGC:
atol: 1.0e+4 # absolute tolerance of RGC residuum (in Pa) atol: 1.0e+4 # absolute tolerance of RGC residuum (in Pa)
rtol: 1.0e-3 # relative ... rtol: 1.0e-3 # relative ...
amax: 1.0e+10 # absolute upper-limit of RGC residuum (in Pa) amax: 1.0e+10 # absolute upper-limit of RGC residuum (in Pa)
rmax: 1.0e+2 # relative ... rmax: 1.0e+2 # relative ...
perturbpenalty: 1.0e-7 # perturbation for computing penalty tangent perturbpenalty: 1.0e-7 # perturbation for computing penalty tangent
relevantmismatch: 1.0e-5 # minimum threshold of mismatch relevantmismatch: 1.0e-5 # minimum threshold of mismatch
viscositypower: 1.0e+0 # power (sensitivity rate) of numerical viscosity in RGC scheme 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) 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 # 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) 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) maxrelaxationrate: 1.0e+0 # threshold of maximum relaxation vector increment (if exceed this then cutback)
maxvoldiscrepancy: 1.0e-5 # maximum allowable relative volume discrepancy maxvoldiscrepancy: 1.0e-5 # maximum allowable relative volume discrepancy
voldiscrepancymod: 1.0e+12 voldiscrepancymod: 1.0e+12
discrepancypower: 5.0 discrepancypower: 5.0
generic: generic:
subStepMin: 1.0e-3 # minimum (relative) size of sub-step allowed during cutback in homogenization 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) 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) stepIncrease: 1.5 # increase of next substep size when previous substep converged in homogenization (value higher than 1)
nMPstate: 10 # materialpoint state loop limit nMPstate: 10 # materialpoint state loop limit
grid: grid:
eps_div_atol: 1.0e-3 # absolute tolerance for fulfillment of stress equilibrium eps_div_atol: 1.0e-3 # absolute tolerance for fulfillment of stress equilibrium
@ -39,30 +39,30 @@ grid:
itmax: 250 # Maximum iteration number itmax: 250 # Maximum iteration number
itmin: 2 # Minimum 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_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 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) maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
maxStaggeredIter: 10 # max number of field level staggered iterations maxStaggeredIter: 10 # max number of field level staggered iterations
memory_efficient: 1 # Precalculate Gamma-operator (81 double per point) 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) 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 derivative: continuous # Approximation used for derivatives in Fourier space
solver: Basic # Type of spectral solver (BasicPETSc/Polarisation/FEM) solver: Basic # Type of spectral solver (BasicPETSc/Polarisation/FEM)
petsc_options: -snes_type ngmres -snes_ngmres_anderson # PetSc solver options 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 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 beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
mesh: mesh:
maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc) maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
maxStaggeredIter: 10 # max number of field level staggered iterations maxStaggeredIter: 10 # max number of field level staggered iterations
structorder: 2 # order of displacement shape functions (when mesh is defined) 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) integrationorder: 2 # order of quadrature rule required (when mesh is defined)
itmax: 250 # Maximum iteration number itmax: 250 # Maximum iteration number
itmin: 2 # Minimum iteration number itmin: 2 # Minimum iteration number
eps_struct_atol: 1.0e-10 # absolute tolerance for mechanical equilibrium eps_struct_atol: 1.0e-10 # absolute tolerance for mechanical equilibrium
eps_struct_rtol: 1.0e-4 # relative tolerance for mechanical equilibrium eps_struct_rtol: 1.0e-4 # relative tolerance for mechanical equilibrium
crystallite: crystallite:
subStepMin: 1.0e-3 # minimum (relative) size of sub-step allowed during cutback in 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) subStepSize: 0.25 # size of substep when cutback introduced in crystallite (value between 0 and 1)
@ -84,5 +84,3 @@ generic:
charLength: 1.0 # characteristic length scale for gradient problems. charLength: 1.0 # characteristic length scale for gradient problems.
random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed. random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed.
residualStiffness: 1.0e-6 # non-zero residual damage. residualStiffness: 1.0e-6 # non-zero residual damage.

41
src/crystallite.f90
View File

@ -69,8 +69,6 @@ module crystallite
real(pReal), dimension(:,:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_partionedF !< def grad to be reached at end of homog inc 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 :: & logical, dimension(:,:,:), allocatable, public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
logical, dimension(:,:,:), allocatable :: & logical, dimension(:,:,:), allocatable :: &
@ -183,8 +181,6 @@ subroutine crystallite_init
crystallite_subFp0,crystallite_subFi0, & crystallite_subFp0,crystallite_subFi0, &
source = crystallite_partionedF) 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_dt(cMax,iMax,eMax),source=0.0_pReal)
allocate(crystallite_subdt,crystallite_subFrac,crystallite_subStep, & allocate(crystallite_subdt,crystallite_subFrac,crystallite_subStep, &
source = crystallite_dt) source = crystallite_dt)
@ -293,7 +289,6 @@ subroutine crystallite_init
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
devNull = crystallite_stress() devNull = crystallite_stress()
call crystallite_stressTangent
#ifdef DEBUG #ifdef DEBUG
if (debugCrystallite%basic) then if (debugCrystallite%basic) then
@ -566,12 +561,14 @@ end subroutine crystallite_restore
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF). !> @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 c, & !< counter in constituent loop
i, & !< counter in integration point loop i, & !< counter in integration point loop
e, & !< counter in element loop e !< counter in element loop
integer :: &
o, & o, &
p p
@ -593,12 +590,6 @@ subroutine crystallite_stressTangent
real(pReal), dimension(9,9):: temp_99 real(pReal), dimension(9,9):: temp_99
logical :: error 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, & call constitutive_SandItsTangents(devNull,dSdFe,dSdFi, &
crystallite_Fe(1:3,1:3,c,i,e), & 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_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)) 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 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 enddo
do o=1,3; do p=1,3 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) & 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), & + matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
dFpinvdF(1:3,1:3,p,o)),temp_33_1) & dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), & + matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), &
transpose(invFp)) & transpose(invFp)) &
+ matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o))) + matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo enddo; enddo
enddo; enddo end function crystallite_stressTangent
enddo elementLooping
!$OMP END PARALLEL DO
end subroutine crystallite_stressTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

26
src/homogenization.f90
View File

@ -374,8 +374,6 @@ subroutine materialpoint_stressAndItsTangent(dt)
enddo cutBackLooping enddo cutBackLooping
call crystallite_stressTangent
if (.not. terminallyIll ) then if (.not. terminallyIll ) then
call crystallite_orientations() ! calculate crystal orientations call crystallite_orientations() ! calculate crystal orientations
!$OMP PARALLEL DO !$OMP PARALLEL DO
@ -437,11 +435,16 @@ function updateState(subdt,subF,ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
integer :: c
logical, dimension(2) :: updateState logical, dimension(2) :: updateState
real(pReal) :: dPdFs(3,3,3,3,homogenization_Ngrains(material_homogenizationAt(el)))
updateState = .true. updateState = .true.
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c=1,homogenization_Ngrains(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo
updateState = & updateState = &
updateState .and. & updateState .and. &
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
@ -449,7 +452,7 @@ function updateState(subdt,subF,ip,el)
crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),& crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),&
subF,& subF,&
subdt, & subdt, &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & dPdFs, &
ip, & ip, &
el) el)
end select chosenHomogenization end select chosenHomogenization
@ -483,26 +486,35 @@ subroutine averageStressAndItsTangent(ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number 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))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) 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 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(& call mech_isostrain_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,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_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))) homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization 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(& call mech_RGC_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,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_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))) homogenization_typeInstance(material_homogenizationAt(el)))
end select chosenHomogenization end select chosenHomogenization