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

improved AL solver, now using guesses for P(x) to improve performance. Changes (and whole solver) still experimental

This commit is contained in:
Martin Diehl 2012-03-31 12:41:46 +00:00
parent 91a70b0fb3
commit 990f547091
1 changed files with 124 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

204
code/DAMASK_spectral_AL.f90
View File

@ -111,11 +111,11 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), dimension(3,3) :: P_av = 0.0_pReal, P_star_av = 0.0_pReal, P, & real(pReal), dimension(3,3) :: P_star_av = 0.0_pReal, &
F_aim = math_I3, F_aim_lastInc = math_I3, lambda_av, & F_aim = math_I3, F_aim_lastInc = math_I3, lambda_av, &
mask_stress, mask_defgrad, deltaF, F_star_av, & mask_stress, mask_defgrad, deltaF, F_star_av, &
F_aim_lab ! quantities rotated to other coordinate system F_aim_lab ! quantities rotated to other coordinate system
real(pReal), dimension(3,3,3,3) :: dPdF, C_inc0, C=0.0_pReal, S_lastInc, C_lastInc ! stiffness and compliance real(pReal), dimension(3,3,3,3) :: C_inc0, C=0.0_pReal, S_lastInc, C_lastInc ! stiffness and compliance
real(pReal), dimension(6) :: sigma ! cauchy stress real(pReal), dimension(6) :: sigma ! cauchy stress
real(pReal), dimension(6,6) :: dsde ! small strain stiffness real(pReal), dimension(6,6) :: dsde ! small strain stiffness
real(pReal), dimension(9,9) :: s_prev99, c_prev99 ! compliance and stiffness in matrix notation real(pReal), dimension(9,9) :: s_prev99, c_prev99 ! compliance and stiffness in matrix notation
@ -127,7 +127,8 @@ program DAMASK_spectral_AL
type(C_PTR) :: tensorField ! fields in real an fourier space type(C_PTR) :: tensorField ! fields in real an fourier space
real(pReal), dimension(:,:,:,:,:), pointer :: lambda_real, F_real ! fields in real space (pointer) real(pReal), dimension(:,:,:,:,:), pointer :: lambda_real, F_real ! fields in real space (pointer)
complex(pReal), dimension(:,:,:,:,:), pointer :: lambda_fourier, F_fourier ! fields in fourier space (pointer) complex(pReal), dimension(:,:,:,:,:), pointer :: lambda_fourier, F_fourier ! fields in fourier space (pointer)
real(pReal), dimension(:,:,:,:,:), allocatable :: F_lastInc, F_star, lambda real(pReal), dimension(:,:,:,:,:), allocatable :: F_lastInc, F_star, lambda, P, F_star_lastIter
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: dPdF
real(pReal), dimension(:,:,:,:), allocatable :: coordinates real(pReal), dimension(:,:,:,:), allocatable :: coordinates
real(pReal), dimension(:,:,:), allocatable :: temperature real(pReal), dimension(:,:,:), allocatable :: temperature
@ -142,14 +143,15 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! loop variables, convergence etc.
real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc = 1.0_pReal, timeinc_old = 0.0_pReal ! elapsed time, begin of interval, time interval real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc = 1.0_pReal, timeinc_old = 0.0_pReal ! elapsed time, begin of interval, time interval
real(pReal) :: guessmode, err_stress, err_stress_tol, err_f, err_p, err_crit, err_f_point, pstress_av_L2, err_div_rms, err_div real(pReal) :: guessmode, err_stress, err_stress_tol, err_f, err_p, err_crit, &
err_f_point, err_p_point, pstress_av_L2, err_div_rms, err_div
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal
complex(pReal), dimension(3,3) :: temp33_Complex complex(pReal), dimension(3,3) :: temp33_Complex
real(pReal), dimension(3,3) :: temp33_Real real(pReal), dimension(3,3) :: temp33_Real
integer(pInt) :: i, j, k, l, u, v, w, errorID = 0_pInt, ierr integer(pInt) :: i, j, k, l, u, v, w, errorID = 0_pInt, ierr
integer(pInt) :: N_Loadcases, loadcase, inc, iter, ielem, CPFEM_mode, & integer(pInt) :: N_Loadcases, loadcase, inc, iter, ielem, CPFEM_mode, guesses, guessmax=10_pInt,&
totalIncsCounter = 0_pInt,notConvergedCounter = 0_pInt, convergedCounter = 0_pInt totalIncsCounter = 0_pInt,notConvergedCounter = 0_pInt, convergedCounter = 0_pInt
logical :: errmatinv logical :: errmatinv, callCPFEM
character(len=6) :: loadcase_string character(len=6) :: loadcase_string
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -196,8 +198,6 @@ program DAMASK_spectral_AL
100 N_Loadcases = N_n 100 N_Loadcases = N_n
if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check
call IO_error(error_ID=837_pInt,ext_msg = trim(path)) ! error message for incomplete loadcase call IO_error(error_ID=837_pInt,ext_msg = trim(path)) ! error message for incomplete loadcase
!if (N_Loadcases>9000_pInt) stop !discuss with Philip, stop code trouble. suggesting warning
allocate (bc(N_Loadcases)) allocate (bc(N_Loadcases))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -416,7 +416,10 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate more memory ! allocate more memory
allocate (P ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (dPdF ( res(1), res(2),res(3),3,3,3,3), source = 0.0_pReal)
allocate (F_star ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F_star ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (F_star_lastIter ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (lambda ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (lambda ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal)
@ -496,15 +499,15 @@ program DAMASK_spectral_AL
F_real(i,j,k,1:3,1:3) = math_I3; F_lastInc(i,j,k,1:3,1:3) = math_I3 F_real(i,j,k,1:3,1:3) = math_I3; F_lastInc(i,j,k,1:3,1:3) = math_I3
coordinates(i,j,k,1:3) = geomdim/real(res * [i,j,k], pReal) - geomdim/real(2_pInt*res,pReal) coordinates(i,j,k,1:3) = geomdim/real(res * [i,j,k], pReal) - geomdim/real(2_pInt*res,pReal)
call CPFEM_general(3_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),& call CPFEM_general(3_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),&
0.0_pReal,ielem,1_pInt,sigma,dsde,temp33_Real ,dPdF) 0.0_pReal,ielem,1_pInt,sigma,dsde,temp33_Real ,dPdF(i,j,k,1:3,1:3,1:3,1:3))
enddo; enddo; enddo enddo; enddo; enddo
ielem = 0_pInt ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(2_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),& call CPFEM_general(2_pInt,coordinates(i,j,k,1:3),math_I3,math_I3,temperature(i,j,k),&
0.0_pReal,ielem,1_pInt,sigma,dsde,temp33_Real ,dPdF) 0.0_pReal,ielem,1_pInt,sigma,dsde,temp33_Real ,dPdF(i,j,k,1:3,1:3,1:3,1:3))
C = C + dPdF C = C + dPdF(i,j,k,1:3,1:3,1:3,1:3)
enddo; enddo; enddo enddo; enddo; enddo
C_inc0 = C * wgt ! linear reference material stiffness C_inc0 = C * wgt ! linear reference material stiffness
@ -614,14 +617,6 @@ program DAMASK_spectral_AL
F_aim_lastInc = temp33_Real F_aim_lastInc = temp33_Real
F_star_av = F_aim F_star_av = F_aim
!--------------------------------------------------------------------------------------------------
! Initialize / Update lambda to useful value
temp33_real = math_mul3333xx33(C*wgt, F_aim-F_aim_lastInc)
P_av = P_av + temp33_real
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + temp33_real
enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update local deformation gradient ! update local deformation gradient
deltaF = math_rotate_backward33(deltaF,bc(loadcase)%rotation) deltaF = math_rotate_backward33(deltaF,bc(loadcase)%rotation)
@ -634,6 +629,16 @@ program DAMASK_spectral_AL
F_lastInc(i,j,k,1:3,1:3) = temp33_Real F_lastInc(i,j,k,1:3,1:3) = temp33_Real
enddo; enddo; enddo enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! Initialize / Update lambda to useful value
P_star_av = P_star_av + math_mul3333xx33(C*wgt, F_aim-F_aim_lastInc)
lambda_av = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
lambda(i,j,k,1:3,1:3) = P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3), &
F_real(i,j,k,1:3,1:3)-F_lastInc(i,j,k,1:3,1:3))
lambda_av = lambda_av + lambda(i,j,k,1:3,1:3)
enddo; enddo; enddo
lambda_av=lambda_av*wgt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!Initialize pointwise data for AL scheme: ToDo: good choice? !Initialize pointwise data for AL scheme: ToDo: good choice?
@ -678,13 +683,15 @@ program DAMASK_spectral_AL
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
CPFEM_mode = 1_pInt ! winding forward CPFEM_mode = 1_pInt ! winding forward
iter = 0_pInt iter = 0_pInt
err_crit = huge(err_div_tol) ! go into loop err_crit = huge(1.0_pReal) ! go into loop
callCPFEM=.true.
guessmax = 13
guesses = 0
!################################################################################################## !##################################################################################################
! convergence loop (looping over iterations) ! convergence loop (looping over iterations)
!################################################################################################## !##################################################################################################
do while((iter < itmax .and. (err_div > err_div_tol .or. err_stress > err_stress_tol .or. err_crit > 5.0e-4))& do while((iter < itmax .and. (err_crit > 1.0_pReal)) .or. iter < itmin)
.or. iter < itmin)
iter = iter + 1_pInt iter = iter + 1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -697,11 +704,9 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
if(size_reduced > 0_pInt) then ! calculate stress BC if applied if(size_reduced > 0_pInt) then ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (P_av - bc(loadcase)%P))) ! maximum deviaton (tensor norm not applicable) err_stress = maxval(abs(mask_stress * (lambda_av - bc(loadcase)%P))) ! maximum deviaton (tensor norm not applicable)
write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'stress deviation =',& F_aim = F_aim + math_mul3333xx33(S_lastInc,bc(loadcase)%P- lambda_av)
math_transpose33(mask_stress * (P_av - bc(loadcase)%P))/1.0e6_pReal err_stress_tol = maxval(abs(lambda_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
F_aim = F_aim + math_mul3333xx33(S_lastInc,bc(loadcase)%P- P_av)
err_stress_tol = maxval(abs(P_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
else else
err_stress_tol = + huge(1.0_pReal) err_stress_tol = + huge(1.0_pReal)
endif endif
@ -746,6 +751,10 @@ program DAMASK_spectral_AL
enddo; enddo enddo; enddo
err_div_RMS = sqrt(err_div_RMS)*wgt err_div_RMS = sqrt(err_div_RMS)*wgt
! if (err_div < err_div_RMS/pstress_av_L2 .and. guessmax<0) then
! print*, 'increasing div, stopping calc'
! iter = huge(1_pInt)
! endif
err_div = err_div_RMS/pstress_av_L2 err_div = err_div_RMS/pstress_av_L2
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! using gamma operator to update F ! using gamma operator to update F
@ -783,44 +792,69 @@ program DAMASK_spectral_AL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! !
print '(a)', '... update stress field P(F*) and update F* and λ..........................' if(callCPFEM) then
print '(a)', '... calling CPFEM to update P(F*) and F*.........................'
ielem = 0_pInt ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(3_pInt,& ! collect cycle call CPFEM_general(3_pInt,& ! collect cycle
coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),& coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),&
F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,& F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
sigma,dsde, P, dPdF) sigma,dsde, P(i,j,k,1:3,1:3), dPdF(i,j,k,1:3,1:3,1:3,1:3))
enddo; enddo; enddo enddo; enddo; enddo
ielem = 0_pInt ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,&
coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),&
F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
sigma,dsde, P(i,j,k,1:3,1:3), dPdF(i,j,k,1:3,1:3,1:3,1:3))
CPFEM_mode = 2_pInt ! winding forward
temp33_Real = lambda(i,j,k,1:3,1:3) - P(i,j,k,1:3,1:3) &
+ math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3))
F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) + math_mul3333xx33(math_invSym3333(&
C_inc0 + dPdF(i,j,k,1:3,1:3,1:3,1:3)), temp33_Real)
enddo; enddo; enddo
F_star_lastIter = F_star
else
guesses = guesses +1_pInt
print '(a)', '... .linear solution.P(F*) and update F*........................'
print*, '... linear approximation ', guesses, ' of ', guessmax, ' for P(F*) and F*'
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
temp33_Real = lambda(i,j,k,1:3,1:3) - (P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3),&
F_star(i,j,k,1:3,1:3) -F_star_lastIter(i,j,k,1:3,1:3)))&
+ math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3))
F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) + math_mul3333xx33(math_invSym3333(&
C_inc0 + dPdF(i,j,k,1:3,1:3,1:3,1:3)), temp33_Real)
enddo; enddo; enddo
endif
print '(a)', '... update λ..........................'
err_f = 0.0_pReal err_f = 0.0_pReal
err_f_point = 0.0_pReal err_f_point = 0.0_pReal
err_p = 0.0_pReal
err_p_point = 0.0_pReal
F_star_av = 0.0_pReal F_star_av = 0.0_pReal
P_star_av = 0.0_pReal P_star_av = 0.0_pReal
lambda_av = 0.0_pReal lambda_av = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,&
coordinates(i,j,k,1:3),F_lastInc(i,j,k,1:3,1:3), &
F_star(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
sigma,dsde, P,dPdF)
CPFEM_mode = 2_pInt ! winding forward
temp33_Real = lambda(i,j,k,1:3,1:3) - P &
+ math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3)- F_star(i,j,k,1:3,1:3))
F_star(i,j,k,1:3,1:3) = F_star(i,j,k,1:3,1:3) + &
math_mul3333xx33(math_invSym3333(C_inc0 + dPdF), temp33_Real)
lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3) & lambda(i,j,k,1:3,1:3) = lambda(i,j,k,1:3,1:3) + math_mul3333xx33(C_inc0,F_real(i,j,k,1:3,1:3) &
- F_star(i,j,k,1:3,1:3)) - F_star(i,j,k,1:3,1:3))
F_star_av = F_star_av + F_star(i,j,k,1:3,1:3) F_star_av = F_star_av + F_star(i,j,k,1:3,1:3)
lambda_av = lambda_av + lambda(i,j,k,1:3,1:3) lambda_av = lambda_av + lambda(i,j,k,1:3,1:3)
P_star_av = P_star_av + P P_star_av = P_star_av + P(i,j,k,1:3,1:3)
temp33_real = F_star(i,j,k,1:3,1:3) - F_real(i,j,k,1:3,1:3) temp33_real = F_star(i,j,k,1:3,1:3) - F_real(i,j,k,1:3,1:3)
err_f_point = max(err_f_point, maxval(temp33_real)) err_f_point = max(err_f_point, maxval(abs(temp33_real)))
err_f = max(err_f, sqrt(math_mul33xx33(temp33_real,temp33_real))) err_f = max(err_f, sqrt(math_mul33xx33(temp33_real,temp33_real)))
temp33_real = lambda(i,j,k,1:3,1:3) - P temp33_real = lambda(i,j,k,1:3,1:3) - (P(i,j,k,1:3,1:3) + math_mul3333xx33(dPdF(i,j,k,1:3,1:3,1:3,1:3),&
F_star(i,j,k,1:3,1:3) -F_star_lastIter(i,j,k,1:3,1:3)))
err_p_point = max(err_p_point, maxval(abs(temp33_real)))
err_p = max(err_p, sqrt(math_mul33xx33(temp33_real,temp33_real))) err_p = max(err_p, sqrt(math_mul33xx33(temp33_real,temp33_real)))
enddo; enddo; enddo enddo; enddo; enddo
@ -833,52 +867,34 @@ program DAMASK_spectral_AL
lambda_av = lambda_av *wgt lambda_av = lambda_av *wgt
write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'λ / GPa =',& write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'λ / GPa =',&
math_transpose33(lambda_av) /1.e6_pReal math_transpose33(lambda_av) /1.e6_pReal
! print '(a)', '... update stress field P(F) .....................................'
! ielem = 0_pInt
! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
! ielem = ielem + 1_pInt
! call CPFEM_general(3_pInt,& ! collect cycle
! coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),&
! F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
! sigma,dsde,P,dPdF)
! enddo; enddo; enddo
! ielem = 0_pInt
! err_p = 0.0_pReal
! P_av =0.0_pReal
! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
! ielem = ielem + 1_pInt
! call CPFEM_general(2_pInt,&
! coordinates(i,j,k,1:3),F_lastInc(i,j,k,1:3,1:3), &
! F_real(i,j,k,1:3,1:3),temperature(i,j,k),timeinc,ielem,1_pInt,&
! sigma,dsde,P,dPdF)
! P_av = P_av + P
! temp33_real = lambda(i,j,k,1:3,1:3) - P
! err_p = max(err_p, sqrt(math_mul33xx33(temp33_real,temp33_real)))
! enddo; enddo; enddo
! P_av = math_rotate_forward33(P_av * wgt,bc(loadcase)%rotation)
! write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F) / GPa =',&
! math_transpose33(P_av)/1.e6_pReal
! write (*,'(a,/,3(3(es14.7,1x)/))',advance='no') 'P(F*) - P(F) =',&
! math_transpose33(P_star_av - P_av)
P_av = lambda_av
err_f = err_f/sqrt(math_mul33xx33(F_star_av,F_star_av)) err_f = err_f/sqrt(math_mul33xx33(F_star_av,F_star_av))
err_p = err_p/sqrt(math_mul33xx33(P_av,P_av)) err_p = err_p/sqrt(math_mul33xx33(P_star_av,P_star_av))
write(6,'(a,es14.7,1x)') 'error F', err_f
write(6,'(a,es14.7,1x)') 'max abs err F', err_f_point
write(6,'(a,es14.7,1x)') 'error P', err_p
write(6,'(a,es11.4)') 'error divergence FT RMS = ',err_div_RMS
write(6,'(a,es11.4)') 'error div = ',err_div
write(6,'(a,es11.4)') 'error stress = ',err_stress/err_stress_tol
err_crit = max(err_p, err_f) write(6,'(a,es14.7,es14.7)') 'error F', err_f/1e-4, err_f
write(6,'(a,es14.7,es14.7)') 'error P', err_p/1e-3, err_p
write(6,'(a,es14.7,es14.7)') 'error stress = ',err_stress/err_stress_tol, err_stress
write(6,'(a,es14.7,es14.7)') 'error divergence = ',err_div/err_div_tol, err_div
write(6,*) ' '
write(6,'(a,es14.7)') 'error divergence FT RMS = ',err_div_RMS
write(6,'(a,es14.7)') 'max abs err F', err_f_point
write(6,'(a,es14.7)') 'max abs err P', err_p_point
err_crit = max(err_p/1e-3, err_f/1e-4,err_div/err_div_tol,err_stress/err_stress_tol)
print*, 'critical error', err_crit
if (.not. callCPFEM) then
if(err_crit < 1.0_pReal .or. guesses >= guessmax) callCPFEM = .true.
err_crit =huge(1.0_pReal)
else
if(guessmax > 1 .and. iter>2) callCPFEM=.false.
guesses = 0_pInt
guessmax = guessmax -1
endif
enddo ! end looping when convergency is achieved enddo ! end looping when convergency is achieved
print '(a)', '' print '(a)', ''
print '(a)', '==================================================================' print '(a)', '=================================================================='
if(err_crit > err_div_tol .or. err_stress > err_stress_tol) then if(err_crit > 1.0_pReal) then
print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' NOT converged' print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' NOT converged'
notConvergedCounter = notConvergedCounter + 1_pInt notConvergedCounter = notConvergedCounter + 1_pInt
else else
@ -914,7 +930,7 @@ program DAMASK_spectral_AL
notConvergedCounter + convergedCounter, ' increments did not converge!' notConvergedCounter + convergedCounter, ' increments did not converge!'
close(538) close(538)
call fftw_destroy_plan(plan_lambda); call fftw_destroy_plan(plan_correction) call fftw_destroy_plan(plan_lambda); call fftw_destroy_plan(plan_correction)
call quit(0_pInt) call quit(1_pInt)
end program DAMASK_spectral_AL end program DAMASK_spectral_AL
!******************************************************************** !********************************************************************
@ -927,11 +943,21 @@ subroutine quit(stop_id)
implicit none implicit none
integer(pInt), intent(in) :: stop_id integer(pInt), intent(in) :: stop_id
integer, dimension(8) :: dateAndTime ! type default integer
if (stop_id == 0_pInt) stop 0 ! normal termination call date_and_time(values = dateAndTime)
if (stop_id <= 9000_pInt) then ! trigger regridding write(6,'(/,a)') 'DAMASK_spectral_AL terminated on:'
write(6,'(i4)') stop_id write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',&
stop 1 dateAndTime(2),'/',&
dateAndTime(1)
write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',&
dateAndTime(6),':',&
dateAndTime(7)
write(6,'(/,a)') 'Exit code:'
if (stop_id == 1_pInt) stop 1 ! normal termination
if (stop_id <= 0_pInt) then ! trigger regridding
write(6,'(a,i6)') 'restart a', stop_id*(-1_pInt)
stop 2
endif endif
stop 'abnormal termination of DAMASK_spectral' stop 0 ! error
end subroutine end subroutine