From 0f70a19266227b6ec6017551659cb95507cf40b5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 13 Feb 2020 11:17:31 +0100 Subject: [PATCH] Fp matters, not Fp^-1 mathematically absolutely equivalent, but numerically not. Sometikes makes a huge difference in convergence behavior, even though abs(det(Fp)-1) is in the order of 1e-15 --- src/crystallite.f90 | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 66474478f..781448b8c 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -239,9 +239,6 @@ subroutine crystallite_init myNcomponents = homogenization_Ngrains(material_homogenizationAt(e)) do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, myNcomponents crystallite_Fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! plastic def gradient reflects init orientation - crystallite_Fp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) & - / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) - crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e) crystallite_F0(1:3,1:3,c,i,e) = math_I3 crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phaseAt(c,e)) @@ -995,9 +992,9 @@ logical function integrateStress(ipc,ip,el,timeFraction) enddo LiLoop invFp_new = matmul(invFp_current,B) - invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error + Fp_new = Fp_new / math_det33(Fp_new)**(1.0_pReal/3.0_pReal) ! regularize Fe_new = matmul(matmul(F,invFp_new),invFi_new) integrateStress = .true.