From 6702cb3baae2df98d87d02e6f7cc283191d03fcb Mon Sep 17 00:00:00 2001 From: Denny Tjahjanto Date: Tue, 17 Nov 2009 13:42:38 +0000 Subject: [PATCH] Introduction of a numerical viscosity into the RGC scheme to improve the convergent behavior and the numerical stability of the scheme (see: homogenization_RGC.f90). These changes have been incorporated in all related subroutines (homogenization.f90, IO.f90, numerics.config and numerics.f90). --- code/IO.f90 | 4 +++ code/homogenization.f90 | 2 ++ code/homogenization_RGC.f90 | 56 +++++++++++++++++++++++++++---------- code/numerics.config | 3 ++ code/numerics.f90 | 32 ++++++++++++++++----- 5 files changed, 75 insertions(+), 22 deletions(-) diff --git a/code/IO.f90 b/code/IO.f90 index 2e0e6a2cc..5b98111df 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1047,6 +1047,10 @@ endfunction msg = 'Non-positive penalty perturbation in RGC' case (277) msg = 'Non-positive relevant mismatch in RGC' + case (278) + msg = 'Non-positive positive definite viscosity model in RGC' + case (288) + msg = 'Non-positive maximum threshold of relaxation change in RGC' !* Error message when selected perturbation method is not defined case (299) diff --git a/code/homogenization.f90 b/code/homogenization.f90 index 74d27717f..d3f4adf91 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -567,10 +567,12 @@ function homogenization_updateState(& !* RGC homogenization added <<>> case (homogenization_RGC_label) homogenization_updateState = homogenization_RGC_updateState( homogenization_state(ip,el), & + homogenization_subState0(ip,el), & crystallite_P(:,:,:,ip,el), & crystallite_partionedF(:,:,:,ip,el), & crystallite_partionedF0(:,:,:,ip,el),& materialpoint_subF(:,:,ip,el),& + materialpoint_subdt(ip,el), & crystallite_dPdF(:,:,:,:,:,ip,el), & ip, & el) diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index 2d3ea6ea7..2cdfe2652 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -241,11 +241,13 @@ endsubroutine !******************************************************************** function homogenization_RGC_updateState(& state, & ! my state + state0, & ! my state at the beginning of increment ! P, & ! array of current grain stresses F, & ! array of current grain deformation gradients F0, & ! array of initial grain deformation gradients avgF, & ! average deformation gradient + dt, & ! time increment dPdF, & ! array of current grain stiffnesses ip, & ! my integration point el & ! my element @@ -254,18 +256,23 @@ function homogenization_RGC_updateState(& use prec, only: pReal,pInt,p_vec use math, only: math_invert use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,homogenization_typeInstance,homogenization_Ngrains - use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC + use material, only: homogenization_maxNgrains,homogenization_typeInstance, & + homogenization_Ngrains + use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC, & + maxdRelax_RGC,ratePower_RGC,viscModus_RGC use FEsolving, only: theInc,cycleCounter,theTime implicit none !* Definition of variables type(p_vec), intent(inout) :: state - real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P,F,F0 + type(p_vec), intent(in) :: state0 + real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: P,F,F0 real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF real(pReal), dimension (3,3), intent(in) :: avgF - integer(pInt), intent(in) :: ip,el + real(pReal), intent(in) :: dt + integer(pInt), intent(in) :: ip,el +! logical, dimension(2) :: homogenization_RGC_updateState integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc @@ -278,9 +285,8 @@ function homogenization_RGC_updateState(& logical error,RGCdebug,RGCdebugJacobi,RGCcheck ! integer(pInt), parameter :: nFace = 6 - real(pInt), parameter :: max_drelax = 0.5_pReal ! - real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix + real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax RGCcheck = (el == 1 .and. ip == 1) @@ -297,7 +303,8 @@ function homogenization_RGC_updateState(& allocate(resid(3*nIntFaceTot)); resid = 0.0_pReal allocate(tract(nIntFaceTot,3)); tract = 0.0_pReal allocate(relax(3*nIntFaceTot)); relax = state%p(1:3*nIntFaceTot) - allocate(drelax(3*nIntFaceTot)); drelax = 0.0_pReal + allocate(drelax(3*nIntFaceTot)); & + drelax = state%p(1:3*nIntFaceTot) - state0%p(1:3*nIntFaceTot) !* Debugging the obtained state if (RGCdebug) then write(6,'(x,a30)')'Obtained state: ' @@ -337,11 +344,13 @@ function homogenization_RGC_updateState(& call homogenization_RGC_getInterface(intFaceP,2*faceID(1)-1,iGr3P) call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal do i = 1,3 ! compute the traction balance at the interface - do j = 1,3 - tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP))*normP(j) & - + (P(i,j,iGrN) + R(i,j,iGrN))*normN(j) - resid(i+3*(iNum-1)) = tract(iNum,i) ! map into 1D residual array - enddo + tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1))/dt))**ratePower_RGC, & + drelax(i+3*(iNum-1))) + do j = 1,3 + tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP))*normP(j) & + + (P(i,j,iGrN) + R(i,j,iGrN))*normN(j) + resid(i+3*(iNum-1)) = tract(iNum,i) ! map into 1D residual array + enddo enddo !* Debugging the residual stress if (RGCdebug) then @@ -470,7 +479,7 @@ function homogenization_RGC_updateState(& call flush(6) endif -!* Compute the Jacobian of the stress-like penalty (penalty tangent) using perturbation technique +!* Compute the Jacobian of the stress-like penalty (penalty tangent) allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot)); pmatrix = 0.0_pReal allocate(p_relax(3*nIntFaceTot)); p_relax = 0.0_pReal allocate(p_resid(3*nIntFaceTot)); p_resid = 0.0_pReal @@ -513,9 +522,24 @@ function homogenization_RGC_updateState(& write(6,*)' ' call flush(6) endif + +!* Construct the Jacobian matrix of the numerical viscosity tangent + allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot)); rmatrix = 0.0_pReal + forall (i=1,3*nIntFaceTot) & + rmatrix(i,i) = viscModus_RGC*ratePower_RGC/dt*(abs(drelax(i)/dt))**(ratePower_RGC - 1.0_pReal) +!* Debugging the global Jacobian matrix of numerical viscosity tangent + if (RGCdebugJacobi) then + write(6,'(x,a30)')'Jacobian matrix of penalty' + do i = 1,3*nIntFaceTot + write(6,'(x,100(e10.4,x))')(rmatrix(i,j), j = 1,3*nIntFaceTot) + enddo + write(6,*)' ' + call flush(6) + endif + !* The overall Jacobian matrix (due to constitutive and penalty tangents) - allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix if (RGCdebugJacobi) then write(6,'(x,a30)')'Jacobian matrix (total)' do i = 1,3*nIntFaceTot @@ -537,6 +561,7 @@ function homogenization_RGC_updateState(& endif !* Calculate the state update (i.e., new relaxation vectors) + drelax = 0.0_pReal do i = 1,3*nIntFaceTot do j = 1,3*nIntFaceTot drelax(i) = drelax(i) - jnverse(i,j)*resid(j) @@ -544,7 +569,8 @@ function homogenization_RGC_updateState(& enddo relax = relax + drelax state%p(1:3*nIntFaceTot) = relax - if (any(abs(drelax(:)) > max_drelax)) then + if (any(abs(drelax(:)) > maxdRelax_RGC)) then !* forcing cutback when the incremental change + ! of relaxation vector becomes too large ! state%p(1:3*nIntFaceTot) = 0.0_pReal ! write(6,'(x,a,x,i3,x,a,x,i3,x,a)')'RGC_updateState: ip',ip,'| el',el,'relaxation vectors reset to zero' homogenization_RGC_updateState = (/.true.,.false./) diff --git a/code/numerics.config b/code/numerics.config index e5fcbc966..a45b68a1c 100644 --- a/code/numerics.config +++ b/code/numerics.config @@ -32,5 +32,8 @@ aMax_RGC 1.0e+10 # absolute upper-limit of RGC residuum (in rMax_RGC 1.0e+3 # relative ... perturbPenalty_RGC 1.0e-7 # perturbation for computing penalty tangent relevantMismatch_RGC 1.0e-5 # minimum threshold of mismatch +viscosityRate_RGC 1.0e+0 # power (sensitivity rate) of numerical viscosity in RGC scheme +viscosityModulus_RGC 0.0e+0 # stress modulus of RGC numerical viscosity (zero = without numerical viscosity) +maxRelaxation_RGC 1.0e+0 # threshold of maximum relaxation vector increment (if exceed this then cutback) fixed_seed 1234 # put any number larger than zero, integer, if you want to have a pseudo random distribution diff --git a/code/numerics.f90 b/code/numerics.f90 index ef4773553..9689817f2 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -28,13 +28,16 @@ real(pReal) relevantStrain, & ! strain rTol_crystalliteStress, & ! relative tolerance in crystallite stress loop aTol_crystalliteStress, & ! absolute tolerance in crystallite stress loop -!* RGC parameters: added <<>> +!* RGC parameters: added <<>> absTol_RGC, & ! absolute tolerance of RGC residuum relTol_RGC, & ! relative tolerance of RGC residuum absMax_RGC, & ! absolute maximum of RGC residuum relMax_RGC, & ! relative maximum of RGC residuum pPert_RGC, & ! perturbation for computing RGC penalty tangent - xSmoo_RGC ! RGC penalty smoothing parameter (hyperbolic tangent) + xSmoo_RGC, & ! RGC penalty smoothing parameter (hyperbolic tangent) + ratePower_RGC, & ! power (sensitivity rate) of numerical viscosity in RGC scheme + viscModus_RGC, & ! stress modulus of RGC numerical viscosity + maxdRelax_RGC ! threshold of maximum relaxation vector increment (if exceed this then cutback) !* Random seeding parameters: added <<>> integer(pInt) fixedSeed ! fixed seeding for pseudo-random number generator @@ -98,14 +101,17 @@ subroutine numerics_init() rTol_crystalliteStress = 1.0e-6_pReal aTol_crystalliteStress = 1.0e-8_pReal ! residuum is in Lp (hence strain on the order of 1e-8 here) -!* RGC parameters: added <<>> with moderate setting +!* RGC parameters: added <<>> with moderate setting absTol_RGC = 1.0e+4 relTol_RGC = 1.0e-3 absMax_RGC = 1.0e+10 relMax_RGC = 1.0e+2 pPert_RGC = 1.0e-7 xSmoo_RGC = 1.0e-5 - + ratePower_RGC = 1.0e+0 ! Newton viscosity (linear model) + viscModul_RGC = 0.0e+0 ! No viscosity is applied + maxdRelax_RGC = 1.0e+0 + !* Random seeding parameters: added <<>> fixedSeed = 0_pInt @@ -164,7 +170,7 @@ subroutine numerics_init() case ('atol_crystallitestress') aTol_crystalliteStress = IO_floatValue(line,positions,2) -!* RGC parameters: added <<>> +!* RGC parameters: added <<>> case ('atol_rgc') absTol_RGC = IO_floatValue(line,positions,2) case ('rtol_rgc') @@ -177,6 +183,12 @@ subroutine numerics_init() pPert_RGC = IO_floatValue(line,positions,2) case ('relevantmismatch_rgc') xSmoo_RGC = IO_floatValue(line,positions,2) + case ('viscosityrate_rgc') + ratePower_RGC = IO_floatValue(line,positions,2) + case ('viscositymodulus_rgc') + viscModul_RGC = IO_floatValue(line,positions,2) + case ('maxrelaxation_rgc') + maxdRelax_RGC = IO_floatValue(line,positions,2) !* Random seeding parameters: added <<>> case ('fixed_seed') @@ -218,13 +230,16 @@ subroutine numerics_init() write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate write(6,*) -!* RGC parameters: added <<>> +!* RGC parameters: added <<>> write(6,'(a24,x,e8.1)') 'aTol_RGC: ',absTol_RGC write(6,'(a24,x,e8.1)') 'rTol_RGC: ',relTol_RGC write(6,'(a24,x,e8.1)') 'aMax_RGC: ',absMax_RGC write(6,'(a24,x,e8.1)') 'rMax_RGC: ',relMax_RGC write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_RGC write(6,'(a24,x,e8.1)') 'relevantMismatch_RGC: ',xSmoo_RGC + write(6,'(a24,x,e8.1)') 'viscosityrate_RGC: ',ratePower_RGC + write(6,'(a24,x,e8.1)') 'viscositymodulus_RGC: ',viscModul_RGC + write(6,'(a24,x,e8.1)') 'maxrelaxation_RGC: ',maxdRelax_RGC write(6,*) !* Random seeding parameters: added <<>> @@ -254,13 +269,16 @@ subroutine numerics_init() if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270) if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271) -!* RGC parameters: added <<>> +!* RGC parameters: added <<>> if (absTol_RGC <= 0.0_pReal) call IO_error(272) if (relTol_RGC <= 0.0_pReal) call IO_error(273) if (absMax_RGC <= 0.0_pReal) call IO_error(274) if (relMax_RGC <= 0.0_pReal) call IO_error(275) if (pPert_RGC <= 0.0_pReal) call IO_error(276) !! oops !! if (xSmoo_RGC <= 0.0_pReal) call IO_error(277) + if (ratePower_RGC < 0.0_pReal) call IO_error(278) + if (viscModul_RGC < 0.0_pReal) call IO_error(278) + if (maxdRelax_RGC <= 0.0_pReal) call IO_error(288) if (fixedSeed <= 0_pInt) write(6,'(a)') 'Random is random!' endsubroutine