From 57174d0abaf3258b74b909df2dcce0ac4e9477d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Sep 2020 10:53:05 +0200 Subject: [PATCH] do not store dPdF at the crystallite level --- examples/ConfigFiles/numerics.yaml | 62 +++++++++++++++--------------- src/crystallite.f90 | 41 +++++++------------- src/homogenization.f90 | 26 +++++++++---- 3 files changed, 63 insertions(+), 66 deletions(-) diff --git a/examples/ConfigFiles/numerics.yaml b/examples/ConfigFiles/numerics.yaml index 2ff31930e..e25c86f5a 100644 --- a/examples/ConfigFiles/numerics.yaml +++ b/examples/ConfigFiles/numerics.yaml @@ -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) @@ -84,5 +84,3 @@ 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. - - diff --git a/src/crystallite.f90 b/src/crystallite.f90 index ad1c2d4a7..0b19afac1 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -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 !-------------------------------------------------------------------------------------------------- diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 9ba5ae6de..de003fbc3 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -374,8 +374,6 @@ subroutine materialpoint_stressAndItsTangent(dt) enddo cutBackLooping - call crystallite_stressTangent - if (.not. terminallyIll ) then call crystallite_orientations() ! calculate crystal orientations !$OMP PARALLEL DO @@ -437,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), & @@ -449,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 @@ -483,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