From 7f8613f6ad32b56aa3c7c144ff103b18ce595a92 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Sep 2020 10:24:24 +0200 Subject: [PATCH 1/3] always update dPdF (was the default anyways) --- examples/ConfigFiles/numerics.yaml | 1 - src/CPFEM.f90 | 11 +---------- src/grid/spectral_utilities.f90 | 2 +- src/homogenization.f90 | 5 ++--- src/mesh/FEM_utilities.f90 | 2 +- 5 files changed, 5 insertions(+), 16 deletions(-) diff --git a/examples/ConfigFiles/numerics.yaml b/examples/ConfigFiles/numerics.yaml index 4180a7d8a..2ff31930e 100644 --- a/examples/ConfigFiles/numerics.yaml +++ b/examples/ConfigFiles/numerics.yaml @@ -78,7 +78,6 @@ 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: diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 05255e2a1..ec36ffe12 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -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 diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index dc6430c72..b32e8aa5f 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -827,7 +827,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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 0e9ca386b..9ba5ae6de 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -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,7 +374,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt) enddo cutBackLooping - if(updateJaco) call crystallite_stressTangent + call crystallite_stressTangent if (.not. terminallyIll ) then call crystallite_orientations() ! calculate crystal orientations diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 4d9786112..4927d0c1c 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -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 From 57174d0abaf3258b74b909df2dcce0ac4e9477d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Sep 2020 10:53:05 +0200 Subject: [PATCH 2/3] 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 From 88b92bf3b1aa3efb98b09fd17770eedb31c85f58 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 3 Oct 2020 08:34:10 +0200 Subject: [PATCH 3/3] updated test results changed calling sequence leads to slightly different results --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 3b498f5cb..64e62f805 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 3b498f5cb3c50e669588106de1b4cdc4c03ffff1 +Subproject commit 64e62f805b5ad784e3397ee5f735aaeb3cc134c2