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

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).

This commit is contained in:
Denny Tjahjanto 2009-11-17 13:42:38 +00:00
parent cb88019aa6
commit 6702cb3baa
5 changed files with 75 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
code/IO.f90
View File

@ -1047,6 +1047,10 @@ endfunction
msg = 'Non-positive penalty perturbation in RGC' msg = 'Non-positive penalty perturbation in RGC'
case (277) case (277)
msg = 'Non-positive relevant mismatch in RGC' 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 !* Error message when selected perturbation method is not defined
case (299) case (299)

2
code/homogenization.f90
View File

@ -567,10 +567,12 @@ function homogenization_updateState(&
!* RGC homogenization added <<<updated 31.07.2009>>> !* RGC homogenization added <<<updated 31.07.2009>>>
case (homogenization_RGC_label) case (homogenization_RGC_label)
homogenization_updateState = homogenization_RGC_updateState( homogenization_state(ip,el), & homogenization_updateState = homogenization_RGC_updateState( homogenization_state(ip,el), &
homogenization_subState0(ip,el), &
crystallite_P(:,:,:,ip,el), & crystallite_P(:,:,:,ip,el), &
crystallite_partionedF(:,:,:,ip,el), & crystallite_partionedF(:,:,:,ip,el), &
crystallite_partionedF0(:,:,:,ip,el),& crystallite_partionedF0(:,:,:,ip,el),&
materialpoint_subF(:,:,ip,el),& materialpoint_subF(:,:,ip,el),&
materialpoint_subdt(ip,el), &
crystallite_dPdF(:,:,:,:,:,ip,el), & crystallite_dPdF(:,:,:,:,:,ip,el), &
ip, & ip, &
el) el)

56
code/homogenization_RGC.f90
View File

@ -241,11 +241,13 @@ endsubroutine
!******************************************************************** !********************************************************************
function homogenization_RGC_updateState(& function homogenization_RGC_updateState(&
state, & ! my state state, & ! my state
state0, & ! my state at the beginning of increment
! !
P, & ! array of current grain stresses P, & ! array of current grain stresses
F, & ! array of current grain deformation gradients F, & ! array of current grain deformation gradients
F0, & ! array of initial grain deformation gradients F0, & ! array of initial grain deformation gradients
avgF, & ! average deformation gradient avgF, & ! average deformation gradient
dt, & ! time increment
dPdF, & ! array of current grain stiffnesses dPdF, & ! array of current grain stiffnesses
ip, & ! my integration point ip, & ! my integration point
el & ! my element el & ! my element
@ -254,18 +256,23 @@ function homogenization_RGC_updateState(&
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use math, only: math_invert use math, only: math_invert
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,homogenization_typeInstance,homogenization_Ngrains use material, only: homogenization_maxNgrains,homogenization_typeInstance, &
use numerics, only: absTol_RGC,relTol_RGC,absMax_RGC,relMax_RGC,pPert_RGC 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 use FEsolving, only: theInc,cycleCounter,theTime
implicit none implicit none
!* Definition of variables !* Definition of variables
type(p_vec), intent(inout) :: state 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,3,3,homogenization_maxNgrains), intent(in) :: dPdF
real(pReal), dimension (3,3), intent(in) :: avgF 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 logical, dimension(2) :: homogenization_RGC_updateState
integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID
integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc
@ -278,9 +285,8 @@ function homogenization_RGC_updateState(&
logical error,RGCdebug,RGCdebugJacobi,RGCcheck logical error,RGCdebug,RGCdebugJacobi,RGCcheck
! !
integer(pInt), parameter :: nFace = 6 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 real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
RGCcheck = (el == 1 .and. ip == 1) RGCcheck = (el == 1 .and. ip == 1)
@ -297,7 +303,8 @@ function homogenization_RGC_updateState(&
allocate(resid(3*nIntFaceTot)); resid = 0.0_pReal allocate(resid(3*nIntFaceTot)); resid = 0.0_pReal
allocate(tract(nIntFaceTot,3)); tract = 0.0_pReal allocate(tract(nIntFaceTot,3)); tract = 0.0_pReal
allocate(relax(3*nIntFaceTot)); relax = state%p(1:3*nIntFaceTot) 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 !* Debugging the obtained state
if (RGCdebug) then if (RGCdebug) then
write(6,'(x,a30)')'Obtained state: ' 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_getInterface(intFaceP,2*faceID(1)-1,iGr3P)
call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal call homogenization_RGC_interfaceNormal(normP,intFaceP) ! get the interface normal
do i = 1,3 ! compute the traction balance at the interface do i = 1,3 ! compute the traction balance at the interface
do j = 1,3 tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1))/dt))**ratePower_RGC, &
tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP))*normP(j) & drelax(i+3*(iNum-1)))
+ (P(i,j,iGrN) + R(i,j,iGrN))*normN(j) do j = 1,3
resid(i+3*(iNum-1)) = tract(iNum,i) ! map into 1D residual array tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP))*normP(j) &
enddo + (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 enddo
!* Debugging the residual stress !* Debugging the residual stress
if (RGCdebug) then if (RGCdebug) then
@ -470,7 +479,7 @@ function homogenization_RGC_updateState(&
call flush(6) call flush(6)
endif 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(pmatrix(3*nIntFaceTot,3*nIntFaceTot)); pmatrix = 0.0_pReal
allocate(p_relax(3*nIntFaceTot)); p_relax = 0.0_pReal allocate(p_relax(3*nIntFaceTot)); p_relax = 0.0_pReal
allocate(p_resid(3*nIntFaceTot)); p_resid = 0.0_pReal allocate(p_resid(3*nIntFaceTot)); p_resid = 0.0_pReal
@ -514,8 +523,23 @@ function homogenization_RGC_updateState(&
call flush(6) call flush(6)
endif 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) !* 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 if (RGCdebugJacobi) then
write(6,'(x,a30)')'Jacobian matrix (total)' write(6,'(x,a30)')'Jacobian matrix (total)'
do i = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot
@ -537,6 +561,7 @@ function homogenization_RGC_updateState(&
endif endif
!* Calculate the state update (i.e., new relaxation vectors) !* Calculate the state update (i.e., new relaxation vectors)
drelax = 0.0_pReal
do i = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot
do j = 1,3*nIntFaceTot do j = 1,3*nIntFaceTot
drelax(i) = drelax(i) - jnverse(i,j)*resid(j) drelax(i) = drelax(i) - jnverse(i,j)*resid(j)
@ -544,7 +569,8 @@ function homogenization_RGC_updateState(&
enddo enddo
relax = relax + drelax relax = relax + drelax
state%p(1:3*nIntFaceTot) = relax 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 ! 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' ! 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./) homogenization_RGC_updateState = (/.true.,.false./)

3
code/numerics.config
View File

@ -32,5 +32,8 @@ aMax_RGC 1.0e+10 # absolute upper-limit of RGC residuum (in
rMax_RGC 1.0e+3 # relative ... rMax_RGC 1.0e+3 # relative ...
perturbPenalty_RGC 1.0e-7 # perturbation for computing penalty tangent perturbPenalty_RGC 1.0e-7 # perturbation for computing penalty tangent
relevantMismatch_RGC 1.0e-5 # minimum threshold of mismatch 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 fixed_seed 1234 # put any number larger than zero, integer, if you want to have a pseudo random distribution

30
code/numerics.f90
View File

@ -28,13 +28,16 @@ real(pReal) relevantStrain, & ! strain
rTol_crystalliteStress, & ! relative tolerance in crystallite stress loop rTol_crystalliteStress, & ! relative tolerance in crystallite stress loop
aTol_crystalliteStress, & ! absolute tolerance in crystallite stress loop aTol_crystalliteStress, & ! absolute tolerance in crystallite stress loop
!* RGC parameters: added <<<updated 31.07.2009>>> !* RGC parameters: added <<<updated 17.11.2009>>>
absTol_RGC, & ! absolute tolerance of RGC residuum absTol_RGC, & ! absolute tolerance of RGC residuum
relTol_RGC, & ! relative tolerance of RGC residuum relTol_RGC, & ! relative tolerance of RGC residuum
absMax_RGC, & ! absolute maximum of RGC residuum absMax_RGC, & ! absolute maximum of RGC residuum
relMax_RGC, & ! relative maximum of RGC residuum relMax_RGC, & ! relative maximum of RGC residuum
pPert_RGC, & ! perturbation for computing RGC penalty tangent 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 <<<updated 27.08.2009>>> !* Random seeding parameters: added <<<updated 27.08.2009>>>
integer(pInt) fixedSeed ! fixed seeding for pseudo-random number generator integer(pInt) fixedSeed ! fixed seeding for pseudo-random number generator
@ -98,13 +101,16 @@ subroutine numerics_init()
rTol_crystalliteStress = 1.0e-6_pReal 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) aTol_crystalliteStress = 1.0e-8_pReal ! residuum is in Lp (hence strain on the order of 1e-8 here)
!* RGC parameters: added <<<updated 31.07.2009>>> with moderate setting !* RGC parameters: added <<<updated 17.11.2009>>> with moderate setting
absTol_RGC = 1.0e+4 absTol_RGC = 1.0e+4
relTol_RGC = 1.0e-3 relTol_RGC = 1.0e-3
absMax_RGC = 1.0e+10 absMax_RGC = 1.0e+10
relMax_RGC = 1.0e+2 relMax_RGC = 1.0e+2
pPert_RGC = 1.0e-7 pPert_RGC = 1.0e-7
xSmoo_RGC = 1.0e-5 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 <<<updated 27.08.2009>>> !* Random seeding parameters: added <<<updated 27.08.2009>>>
fixedSeed = 0_pInt fixedSeed = 0_pInt
@ -164,7 +170,7 @@ subroutine numerics_init()
case ('atol_crystallitestress') case ('atol_crystallitestress')
aTol_crystalliteStress = IO_floatValue(line,positions,2) aTol_crystalliteStress = IO_floatValue(line,positions,2)
!* RGC parameters: added <<<updated 31.07.2009>>> !* RGC parameters: added <<<updated 17.11.2009>>>
case ('atol_rgc') case ('atol_rgc')
absTol_RGC = IO_floatValue(line,positions,2) absTol_RGC = IO_floatValue(line,positions,2)
case ('rtol_rgc') case ('rtol_rgc')
@ -177,6 +183,12 @@ subroutine numerics_init()
pPert_RGC = IO_floatValue(line,positions,2) pPert_RGC = IO_floatValue(line,positions,2)
case ('relevantmismatch_rgc') case ('relevantmismatch_rgc')
xSmoo_RGC = IO_floatValue(line,positions,2) 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 <<<updated 27.08.2009>>> !* Random seeding parameters: added <<<updated 27.08.2009>>>
case ('fixed_seed') case ('fixed_seed')
@ -218,13 +230,16 @@ subroutine numerics_init()
write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate write(6,'(a24,x,i8)') 'nMPstate: ',nMPstate
write(6,*) write(6,*)
!* RGC parameters: added <<<updated 31.07.2009>>> !* RGC parameters: added <<<updated 17.11.2009>>>
write(6,'(a24,x,e8.1)') 'aTol_RGC: ',absTol_RGC 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)') 'rTol_RGC: ',relTol_RGC
write(6,'(a24,x,e8.1)') 'aMax_RGC: ',absMax_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)') 'rMax_RGC: ',relMax_RGC
write(6,'(a24,x,e8.1)') 'perturbPenalty_RGC: ',pPert_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)') '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,*) write(6,*)
!* Random seeding parameters: added <<<updated 27.08.2009>>> !* Random seeding parameters: added <<<updated 27.08.2009>>>
@ -254,13 +269,16 @@ subroutine numerics_init()
if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270) if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270)
if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271) if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271)
!* RGC parameters: added <<<updated 31.07.2009>>> !* RGC parameters: added <<<updated 17.11.2009>>>
if (absTol_RGC <= 0.0_pReal) call IO_error(272) if (absTol_RGC <= 0.0_pReal) call IO_error(272)
if (relTol_RGC <= 0.0_pReal) call IO_error(273) if (relTol_RGC <= 0.0_pReal) call IO_error(273)
if (absMax_RGC <= 0.0_pReal) call IO_error(274) if (absMax_RGC <= 0.0_pReal) call IO_error(274)
if (relMax_RGC <= 0.0_pReal) call IO_error(275) if (relMax_RGC <= 0.0_pReal) call IO_error(275)
if (pPert_RGC <= 0.0_pReal) call IO_error(276) !! oops !! if (pPert_RGC <= 0.0_pReal) call IO_error(276) !! oops !!
if (xSmoo_RGC <= 0.0_pReal) call IO_error(277) 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!' if (fixedSeed <= 0_pInt) write(6,'(a)') 'Random is random!'
endsubroutine endsubroutine