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

changed reference stiffness to min max avg from volume avg. using terminally ill information in convergence check

This commit is contained in:
Pratheek Shanthraj 2013-03-06 14:31:13 +00:00
parent 48e57fc3cb
commit 1f4d7c2ca4
4 changed files with 117 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
code/DAMASK_spectral_solverAL.f90
View File

@ -57,7 +57,7 @@ module DAMASK_spectral_solverAL
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
C = 0.0_pReal, & !< current average stiffness
C = 0.0_pReal, C_minmaxAvg = 0.0_pReal, & !< current average stiffness
C_lastInc = 0.0_pReal, & !< previous average stiffness
S = 0.0_pReal, & !< current compliance (filled up with zeros)
C_scale = 0.0_pReal, &
@ -186,25 +186,25 @@ subroutine AL_init(temperature)
read (777,rec=1) C_lastInc
close (777)
call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(temp3333_Real))
read (777,rec=1) temp3333_Real
read (777,rec=1) C_minmaxAvg
close (777)
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,reshape(&
F,[3,3,res(1),res(2),res(3)])),[3,1,mesh_NcpElems])
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real2,&
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,temp3333_Real,temp3333_Real2,&
temp33_Real,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! reference stiffness
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
temp3333_Real = temp3333_Real2
C = temp3333_Real2
C_minmaxAvg = temp3333_Real2
C = temp3333_Real
endif
call Utilities_updateGamma(temp3333_Real,.True.)
C_scale = temp3333_Real
S_scale = math_invSym3333(temp3333_Real)
call Utilities_updateGamma(temp3333_Real2,.True.)
C_scale = temp3333_Real2
S_scale = math_invSym3333(temp3333_Real2)
end subroutine AL_init
@ -219,7 +219,8 @@ type(tSolutionState) function &
itmax
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
math_rotate_backward33, &
math_invSym3333
use mesh, only: &
res, &
geomdim, &
@ -342,7 +343,11 @@ use mesh, only: &
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C)
if (update_gamma) call Utilities_updateGamma(C,restartWrite)
if (update_gamma) then
call Utilities_updateGamma(C_minmaxAvg,restartWrite)
C_scale = C_minmaxAvg
S_scale = math_invSym3333(C_minmaxAvg)
endif
ForwardData = .True.
mask_stress = P_BC%maskFloat
@ -478,7 +483,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(F_lastInc,F - residual_F_tau/polarBeta,params%temperature,params%timeinc, &
residual_F,C,P_av,ForwardData,params%rotation_BC)
residual_F,C,C_minmaxAvg,P_av,ForwardData,params%rotation_BC)
ForwardData = .False.
!--------------------------------------------------------------------------------------------------
@ -520,6 +525,8 @@ subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr)
err_p_tol, &
err_stress_tolrel, &
err_stress_tolabs
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
@ -538,7 +545,8 @@ subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr)
Converged = (it > itmin .and. &
all([ err_f/sqrt(sum((F_aim-math_I3)**2.0_pReal))/err_f_tol, &
err_p/sqrt(sum((F_aim-math_I3)**2.0_pReal))/err_p_tol, &
err_stress/err_stress_tol] < 1.0_pReal))
err_stress/err_stress_tol] < 1.0_pReal)) &
.or. terminallyIll
if (Converged) then
reason = 1

10
code/DAMASK_spectral_solverBasic.f90
View File

@ -35,7 +35,7 @@ module DAMASK_spectral_SolverBasic
F_aim_lastInc = math_I3, & !< deformation gradient aim last increment
F_aimDot = 0.0_pReal !< assumed rate
real(pReal), private,dimension(3,3,3,3) :: &
C = 0.0_pReal, & !< average stiffness
C = 0.0_pReal, C_minmaxAvg = 0.0_pReal, & !< average stiffness
C_lastInc = 0.0_pReal !< average stiffness last increment
contains
@ -128,9 +128,9 @@ subroutine basic_init(temperature)
close (777)
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomdim,F),[3,1,mesh_NcpElems])
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness
call Utilities_constitutiveResponse(F,F,temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
temp3333_Real = C
temp3333_Real = C_minmaxAvg
endif
call Utilities_updateGamma(temp3333_Real,.True.)
@ -270,7 +270,7 @@ type(tSolutionState) function &
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C)
if (update_gamma) call Utilities_updateGamma(C,restartWrite)
if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite)
!--------------------------------------------------------------------------------------------------
! iteration till converged
@ -294,7 +294,7 @@ type(tSolutionState) function &
F_aim_lastIter = F_aim
basic_solution%termIll = .false.
call Utilities_constitutiveResponse(F_lastInc,F,temperature_bc,timeinc,&
P,C,P_av,ForwardData,rotation_BC)
P,C,C_minmaxAvg,P_av,ForwardData,rotation_BC)
basic_solution%termIll = terminallyIll
terminallyIll = .false.
ForwardData = .false.

15
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -53,7 +53,7 @@ module DAMASK_spectral_SolverBasicPETSc
F_aimDot=0.0_pReal
character(len=1024), private :: incInfo
real(pReal), private, dimension(3,3,3,3) :: &
C = 0.0_pReal, C_lastInc= 0.0_pReal, &
C = 0.0_pReal, C_minmaxAvg = 0.0_pReal, C_lastInc= 0.0_pReal, &
S = 0.0_pReal
real(pReal), private :: err_stress, err_div
@ -182,10 +182,10 @@ subroutine basicPETSc_init(temperature)
call Utilities_constitutiveResponse(&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3)
temperature,0.0_pReal,P,C,C_minmaxAvg,temp33_Real,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back into PETSc
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
temp3333_Real = C
temp3333_Real = C_minmaxAvg
endif
call Utilities_updateGamma(temp3333_Real,.True.)
@ -302,7 +302,7 @@ type(tSolutionState) function &
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C)
if (update_gamma) call Utilities_updateGamma(C,restartWrite)
if (update_gamma) call Utilities_updateGamma(C_minmaxAvg,restartWrite)
ForwardData = .True.
@ -393,7 +393,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(F_lastInc,x_scal,params%temperature,params%timeinc, &
f_scal,C,P_av,ForwardData,params%rotation_BC)
f_scal,C,C_minmaxAvg,P_av,ForwardData,params%rotation_BC)
ForwardData = .false.
!--------------------------------------------------------------------------------------------------
@ -433,6 +433,8 @@ subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ier
math_mul33x33, &
math_eigenvalues33, &
math_transpose33
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
@ -453,7 +455,8 @@ subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ier
pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,math_transpose33(P_av)))))
Converged = (it >= itmin .and. &
all([ err_div/pAvgDivL2/err_div_tol, &
err_stress/err_stress_tol ] < 1.0_pReal))
err_stress/err_stress_tol ] < 1.0_pReal)) &
.or. terminallyIll
if (Converged) then
reason = 1

86
code/DAMASK_spectral_utilities.f90
View File

@ -573,6 +573,66 @@ real(pReal) function utilities_divergenceRMS()
endif
end function utilities_divergenceRMS
!--------------------------------------------------------------------------------------------------
!> @brief calculate max of curl of field_fourier
!--------------------------------------------------------------------------------------------------
real(pReal) function utilities_curlRMS()
use math !< must use the whole module for use of FFTW
use mesh, only: &
res, &
res1_red, &
wgt
implicit none
integer(pInt) :: i, j, k, l
complex(pReal), dimension(3,3) :: curl_fourier
real(pReal) :: curl_abs
write(6,'(/,a)') ' ... calculating curl ................................................'
flush(6)
!--------------------------------------------------------------------------------------------------
! calculating max curl criterion in Fourier space
utilities_curlRMS = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2);
do i = 2_pInt, res1_red - 1_pInt
do l = 1_pInt, 3_pInt
curl_fourier(l,1) = (field_fourier(i,j,k,l,3)*xi(2,i,j,k)&
- field_fourier(i,j,k,l,2)*xi(3,i,j,k))*TWOPIIMG
curl_fourier(l,2) = (-field_fourier(i,j,k,l,3)*xi(1,i,j,k)&
+field_fourier(i,j,k,l,1)*xi(3,i,j,k) )*TWOPIIMG
curl_fourier(l,3) = ( field_fourier(i,j,k,l,2)*xi(1,i,j,k)&
-field_fourier(i,j,k,l,1)*xi(2,i,j,k) )*TWOPIIMG
enddo
utilities_curlRMS = utilities_curlRMS + &
2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)
enddo
do l = 1_pInt, 3_pInt
curl_fourier = (field_fourier(1,j,k,l,3)*xi(2,1,j,k)&
- field_fourier(1,j,k,l,2)*xi(3,1,j,k))*TWOPIIMG
curl_fourier = (-field_fourier(1,j,k,l,3)*xi(1,1,j,k)&
+field_fourier(1,j,k,l,1)*xi(3,1,j,k) )*TWOPIIMG
curl_fourier = ( field_fourier(1,j,k,l,2)*xi(1,1,j,k)&
-field_fourier(1,j,k,l,1)*xi(2,1,j,k) )*TWOPIIMG
enddo
utilities_curlRMS = utilities_curlRMS + &
2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)
do l = 1_pInt, 3_pInt
curl_fourier = (field_fourier(res1_red,j,k,l,3)*xi(2,res1_red,j,k)&
- field_fourier(res1_red,j,k,l,2)*xi(3,res1_red,j,k))*TWOPIIMG
curl_fourier = (-field_fourier(res1_red,j,k,l,3)*xi(1,res1_red,j,k)&
+field_fourier(res1_red,j,k,l,1)*xi(3,res1_red,j,k) )*TWOPIIMG
curl_fourier = ( field_fourier(res1_red,j,k,l,2)*xi(1,res1_red,j,k)&
-field_fourier(res1_red,j,k,l,1)*xi(2,res1_red,j,k) )*TWOPIIMG
enddo
utilities_curlRMS = utilities_curlRMS + &
2.0_pReal*sum(real(curl_fourier)**2.0_pReal + aimag(curl_fourier)**2.0_pReal)
enddo; enddo
utilities_curlRMS = sqrt(utilities_curlRMS) *wgt
end function utilities_curlRMS
!--------------------------------------------------------------------------------------------------
@ -677,7 +737,7 @@ end function utilities_maskedCompliance
!> @brief calculates constitutive response
!--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
P,C,P_av,forwardData,rotation_BC)
P,C_volAvg,C_minmaxAvg,P_av,forwardData,rotation_BC)
use debug, only: &
debug_reset, &
debug_info
@ -713,7 +773,7 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
logical, intent(in) :: forwardData !< age results
real(pReal), intent(in), dimension(3,3) :: rotation_BC !< rotation of load frame
real(pReal),intent(out), dimension(3,3,3,3) :: C !< average stiffness
real(pReal),intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal),intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal),intent(out), dimension(3,3,res(1),res(2),res(3)) :: P !< PK stress
@ -723,6 +783,9 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
real(pReal), dimension(3,3,3,3) :: dPdF !< d P / d F
real(pReal), dimension(6) :: sigma !< cauchy stress in mandel notation
real(pReal), dimension(6,6) :: dsde !< d sigma / d Epsilon
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
real(pReal) :: max_dPdF_norm, min_dPdF_norm
integer(pInt) :: k
write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
calcMode = CPFEM_CALCRESULTS
@ -764,9 +827,26 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
call CPFEM_general(calcMode,F_lastInc(1:3,1:3,1,1,1), F(1:3,1:3,1,1,1), & ! first call calculates everything
temperature,timeinc,1_pInt,1_pInt,sigma,dsde,P(1:3,1:3,1,1,1),dPdF)
max_dPdF = 0.0_pReal
max_dPdF_norm = 0.0_pReal
min_dPdF = huge(1.0_pReal)
min_dPdF_norm = huge(1.0_pReal)
do k = 1_pInt, res(3)*res(2)*res(3)
if (max_dPdF_norm < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then
max_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)
max_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)
endif
if (min_dPdF_norm > sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)) then
min_dPdF = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)
min_dPdF_norm = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,k)**2.0_pReal)
endif
enddo
P = reshape(materialpoint_P, [3,3,res(1),res(2),res(3)])
C = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) * wgt
C_minmaxAvg = 0.5_pReal*(max_dPdF + min_dPdF)
call debug_info()
restartWrite = .false. ! reset restartWrite status