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

corrected calculation of divergence in Fourier space, removed normalization of normdyad (was useless), now using correct compliance for calculation of stress BC.

This commit is contained in:
Martin Diehl 2011-07-07 10:03:55 +00:00
parent 18a5841bc5
commit a561ef1cf5
1 changed files with 76 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

159
code/DAMASK_spectral.f90
View File

@ -56,7 +56,7 @@ program DAMASK_spectral
implicit none implicit none
include 'include/fftw3.f' ! header file for fftw3 (declaring variables). Library files are also needed include 'include/fftw3.f' ! header file for fftw3 (declaring variables). Library files are also needed
! compile FFTW 3.2.2 with ./configure --enable-threads ! compile FFTW 3.2.2 with ./configure --enable-threads
! variables to read from loadcase and geom file ! variables to read from loadcase and geom file
real(pReal), dimension(9) :: valuevector ! stores information temporarily from loadcase file real(pReal), dimension(9) :: valuevector ! stores information temporarily from loadcase file
integer(pInt), parameter :: maxNchunksInput = 26 ! 5 identifiers, 18 values for the matrices and 3 scalars integer(pInt), parameter :: maxNchunksInput = 26 ! 5 identifiers, 18 values for the matrices and 3 scalars
@ -90,7 +90,8 @@ program DAMASK_spectral
pstress, pstress_av, cstress_av, defgrad_av,& pstress, pstress_av, cstress_av, defgrad_av,&
defgradAim, defgradAimOld, defgradAimCorr, defgradAimCorrPrev,& defgradAim, defgradAimOld, defgradAimCorr, defgradAimCorrPrev,&
mask_stress, mask_defgrad mask_stress, mask_defgrad
real(pReal), dimension(3,3,3,3) :: dPdF, c0, s0 real(pReal), dimension(3,3,3,3) :: dPdF, c0, s0, c0_temp
real(pReal), dimension(9,9) :: s099
real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation real(pReal), dimension(6) :: cstress ! cauchy stress in Mandel notation
real(pReal), dimension(6,6) :: dsde, c066, s066 ! Mandel notation of 4th order tensors real(pReal), dimension(6,6) :: dsde, c066, s066 ! Mandel notation of 4th order tensors
real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold
@ -99,14 +100,14 @@ program DAMASK_spectral
! variables storing information for spectral method ! variables storing information for spectral method
complex(pReal) :: img complex(pReal) :: img
complex(pReal), dimension(3,3) :: temp33_Complex complex(pReal), dimension(3,3) :: temp33_Complex
real(pReal), dimension(3,3) :: xinormdyad real(pReal), dimension(3,3) :: xidyad
real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat
real(pReal), dimension(3) :: xi, xi_central real(pReal), dimension(:,:,:,:), allocatable :: xi
integer(pInt), dimension(3) :: k_s integer(pInt), dimension(3) :: k_s
integer*8, dimension(2) :: plan_fft integer*8, dimension(2) :: plan_fft
! loop variables, convergence etc. ! loop variables, convergence etc.
real(pReal) guessmode, err_div, err_stress, err_defgrad, sigma0 real(pReal) guessmode, err_div, err_stress, err_defgrad, pHatAv
integer(pInt) i, j, k, l, m, n, p integer(pInt) i, j, k, l, m, n, p
integer(pInt) loadcase, ielem, iter, calcmode, CPFEM_mode, ierr integer(pInt) loadcase, ielem, iter, calcmode, CPFEM_mode, ierr
logical errmatinv logical errmatinv
@ -281,7 +282,8 @@ program DAMASK_spectral
allocate (defgrad (resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal allocate (defgrad (resolution(1),resolution(2),resolution(3),3,3)); defgrad = 0.0_pReal
allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal allocate (defgradold(resolution(1),resolution(2),resolution(3),3,3)); defgradold = 0.0_pReal
allocate (coordinates(3,resolution(1),resolution(2),resolution(3))); coordinates = 0.0_pReal allocate (coordinates(3,resolution(1),resolution(2),resolution(3))); coordinates = 0.0_pReal
allocate (xi (3,resolution(1)/2+1,resolution(2),resolution(3))); xi = 0.0_pReal
wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal) wgt = 1.0_pReal/real(resolution(1)*resolution(2)*resolution(3), pReal)
defgradAim = math_I3 defgradAim = math_I3
@ -301,49 +303,43 @@ program DAMASK_spectral
c066 = c066 + dsde c066 = c066 + dsde
enddo; enddo; enddo enddo; enddo; enddo
c066 = c066 * wgt c066 = c066 * wgt
c0 = math_mandel66to3333(c066) c0 = math_mandel66to3333(c066) ! linear reference material stiffness
call math_invert(6, c066, s066,i, errmatinv)
if(errmatinv) call IO_error(800) ! Matrix inversion error do k = 1, resolution(3) ! calculation of discrete frequencies, order as in FFTW (wrap around)
s0 = math_mandel66to3333(s066) k_s(3) = k-1
if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3)
do j = 1, resolution(2)
k_s(2) = j-1
if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2)
do i = 1, resolution(1)/2+1
k_s(1) = i-1
xi(3,i,j,k) = 0.0_pReal ! 2D case
if(resolution(3) > 1) xi(3,i,j,k) = real(k_s(3), pReal)*2*pi/geomdimension(3) ! 3D case ToDo: Check if to multiply by 2 pi?
xi(2,i,j,k) = real(k_s(2), pReal)*2*pi/geomdimension(2)
xi(1,i,j,k) = real(k_s(1), pReal)*2*pi/geomdimension(1)
enddo; enddo; enddo
if(memory_efficient) then ! allocate just single fourth order tensor if(memory_efficient) then ! allocate just single fourth order tensor
allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal
else ! precalculation of gamma_hat field else ! precalculation of gamma_hat field
allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 0.0_pReal allocate (gamma_hat(resolution(1)/2+1,resolution(2),resolution(3),3,3,3,3)); gamma_hat = 0.0_pReal
do k = 1, resolution(3) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1
k_s(3) = k-1 if (any(xi(:,i,j,k) /= 0.0_pReal)) then
if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) do l = 1,3; do m = 1,3
do j = 1, resolution(2) xidyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k)
k_s(2) = j-1 enddo; enddo
if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2) temp33_Real = math_inv3x3(math_mul3333xx33(c0, xidyad))
do i = 1, resolution(1)/2+1 else
k_s(1) = i-1 xidyad = 0.0_pReal
xi(3) = 0.0_pReal ! for the 2D case temp33_Real = 0.0_pReal
if(resolution(3) > 1) xi(3) = real(k_s(3), pReal)/geomdimension(3) ! 3D case endif
xi(2) = real(k_s(2), pReal)/geomdimension(2) do l=1,3; do m=1,3; do n=1,3; do p=1,3
xi(1) = real(k_s(1), pReal)/geomdimension(1) gamma_hat(i,j,k, l,m,n,p) = - 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *&
if (any(xi /= 0.0_pReal)) then (xidyad(m,p)+xidyad(p,m))
do l = 1,3; do m = 1,3 enddo; enddo; enddo; enddo
xinormdyad(l,m) = xi(l)*xi(m)/sum(xi**2) ! unit sphere, unit vectors in Fourier space
enddo; enddo
temp33_Real = math_inv3x3(math_mul3333xx33(c0, xinormdyad))
else
xinormdyad = 0.0_pReal
temp33_Real = 0.0_pReal
endif
do l=1,3; do m=1,3; do n=1,3; do p=1,3
gamma_hat(i,j,k, l,m,n,p) = - 0.25*(temp33_Real(l,n)+temp33_Real(n,l)) *&
(xinormdyad(m,p)+xinormdyad(p,m))
enddo; enddo; enddo; enddo
enddo; enddo; enddo enddo; enddo; enddo
endif endif
! calculate xi for the calculation of divergence in Fourier space (central frequency)
xi_central(3) = 0.0_pReal ! 2D case
if(resolution(3) > 1) xi_central(3) = real(resolution(3)/2, pReal)/geomdimension(3) ! 3D case
xi_central(2) = real(resolution(2)/2, pReal)/geomdimension(2)
xi_central(1) = real(resolution(1)/2, pReal)/geomdimension(1)
allocate (workfft(resolution(1)+2,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal allocate (workfft(resolution(1)+2,resolution(2),resolution(3),3,3)); workfft = 0.0_pReal
! Initialization of fftw (see manual on fftw.org for more details) ! Initialization of fftw (see manual on fftw.org for more details)
@ -392,7 +388,7 @@ program DAMASK_spectral
! loop oper steps defined in input file for current loadcase ! loop oper steps defined in input file for current loadcase
do step = 1, bc_steps(loadcase) do step = 1, bc_steps(loadcase)
!************************************************************* !*************************************************************
if (bc_logscale(loadcase) == 1) then ! loglinear scale if (bc_logscale(loadcase) == 1) then ! loglinear scale
if (loadcase == 1) then ! 1st loadcase of loglinear scale if (loadcase == 1) then ! 1st loadcase of loglinear scale
if (step == 1) then ! 1st step of 1st loadcase of loglinear scale if (step == 1) then ! 1st step of 1st loadcase of loglinear scale
timeinc = bc_timeIncrement(1)*(2.0**(1 - bc_steps(1))) ! assume 1st step is equal to 2nd timeinc = bc_timeIncrement(1)*(2.0**(1 - bc_steps(1))) ! assume 1st step is equal to 2nd
@ -444,7 +440,7 @@ program DAMASK_spectral
print*, ' ' print*, ' '
print '(3(A,I5.5,tr2))', ' Loadcase = ',loadcase, ' Step = ',step, ' Iteration = ',iter print '(3(A,I5.5,tr2))', ' Loadcase = ',loadcase, ' Step = ',step, ' Iteration = ',iter
cstress_av = 0.0_pReal cstress_av = 0.0_pReal
workfft = 0.0_pReal !needed because of the padding for FFTW workfft = 0.0_pReal ! needed because of the padding for FFTW
!************************************************************* !*************************************************************
! adjust defgrad to fulfill BCs ! adjust defgrad to fulfill BCs
@ -460,6 +456,7 @@ program DAMASK_spectral
cstress,dsde, pstress, dPdF) cstress,dsde, pstress, dPdF)
enddo; enddo; enddo enddo; enddo; enddo
c0_temp = 0.0_pReal !for calculation of s0
ielem = 0_pInt ielem = 0_pInt
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
@ -469,9 +466,14 @@ program DAMASK_spectral
temperature,timeinc,ielem,1_pInt,& temperature,timeinc,ielem,1_pInt,&
cstress,dsde, pstress, dPdF) cstress,dsde, pstress, dPdF)
CPFEM_mode = 2_pInt CPFEM_mode = 2_pInt
workfft(i,j,k,:,:) = pstress ! prepare Fourier transform of first P--K stress c0_temp = c0_temp + dPdF
workfft(i,j,k,:,:) = pstress ! build up average P-K stress
cstress_av = cstress_av + math_mandel6to33(cstress) ! build up average Cauchy stress cstress_av = cstress_av + math_mandel6to33(cstress) ! build up average Cauchy stress
enddo; enddo; enddo enddo; enddo; enddo
call math_invert(9, math_plain3333to99(c0_temp),s099,i,errmatinv)
if(errmatinv) call IO_error(800,ext_msg = "problem in c0 inversion")
s0 = math_plain99to3333(s099) *real(resolution(1)*resolution(2)*resolution(3), pReal) ! average s0 for calculation of BC
cstress_av = cstress_av * wgt cstress_av = cstress_av * wgt
do n = 1,3; do m = 1,3 do n = 1,3; do m = 1,3
@ -486,7 +488,7 @@ program DAMASK_spectral
defgradAimCorrPrev = defgradAimCorr defgradAimCorrPrev = defgradAimCorr
defgradAimCorr = -mask_stress * math_mul3333xx33(s0, (mask_stress*(pstress_av - bc_stress(:,:,loadcase)))) defgradAimCorr = -mask_stress * math_mul3333xx33(s0, (mask_stress*(pstress_av - bc_stress(:,:,loadcase))))
do m=1,3; do n =1,3 ! calculate damper (correction is far too strong) do m=1,3; do n =1,3 ! calculate damper (correction is far too strong) !ToDo:
if ( sign(1.0_pReal,defgradAimCorr(m,n))/=sign(1.0_pReal,defgradAimCorrPrev(m,n))) then if ( sign(1.0_pReal,defgradAimCorr(m,n))/=sign(1.0_pReal,defgradAimCorrPrev(m,n))) then
damper(m,n) = max(0.01_pReal,damper(m,n)*0.8) damper(m,n) = max(0.01_pReal,damper(m,n)*0.8)
else else
@ -537,47 +539,38 @@ program DAMASK_spectral
enddo; enddo enddo; enddo
print *, 'Calculating equilibrium using spectral method' print *, 'Calculating equilibrium using spectral method'
err_div = 0.0_pReal; sigma0 = 0.0_pReal err_div = 0.0_pReal; pHatAv = 0.0_pReal
call dfftw_execute_dft_r2c(plan_fft(1),workfft,workfft) ! FFT of pstress call dfftw_execute_dft_r2c(plan_fft(1),workfft,workfft) ! FFT of pstress
do m = 1,3 ! L infinity Norm of stress tensor do m = 1,3 ! L infinity Norm of stress tensor
sigma0 = max(sigma0, sum(abs(workfft(1,1,1,m,:) + workfft(2,1,1,m,:)*img))) pHatAv = max(pHatAv, sum(abs(workfft(1,1,1,m,:) + workfft(2,1,1,m,:)*img)))
enddo enddo
err_div = (maxval(abs(math_mul33x3_complex(workfft(resolution(1)+1,resolution(2)/2+1,resolution(3)/2+1,:,:)+& ! L infinity norm of div(stress) do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1
workfft(resolution(1)+2,resolution(2)/2+1,resolution(3)/2+1,:,:)*img,xi_central)))) err_div = max(err_div, maxval(abs(math_mul33x3_complex(workfft(i*2-1,j,k,:,:)+& ! L infinity norm of div(stress)
err_div = err_div/sigma0 ! weighting of error workfft(i*2, j,k,:,:)*img,xi(:,i,j,k)))))
enddo; enddo; enddo
err_div = err_div/pHatAv ! Criterion as supposed in Suquet 2001
if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat
do k = 1, resolution(3) do k = 1, resolution(3); do j = 1, resolution(2) ;do i = 1, resolution(1)/2+1
k_s(3) = k-1 if (any(xi(:,i,j,k) /= 0.0_pReal)) then
if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) do l = 1,3; do m = 1,3
do j = 1, resolution(2) xidyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k)
k_s(2) = j-1 enddo; enddo
if(j > resolution(2)/2+1) k_s(2) = k_s(2)-resolution(2) temp33_Real = math_inv3x3(math_mul3333xx33(c0, xidyad))
do i = 1, resolution(1)/2+1 else
k_s(1) = i-1 xidyad = 0.0_pReal
xi(3) = 0.0_pReal ! for the 2D case temp33_Real = 0.0_pReal
if(resolution(3) > 1) xi(3) = real(k_s(3), pReal)/geomdimension(3) ! 3D case endif
xi(2) = real(k_s(2), pReal)/geomdimension(2) do l=1,3; do m=1,3; do n=1,3; do p=1,3
xi(1) = real(k_s(1), pReal)/geomdimension(1) gamma_hat(1,1,1, l,m,n,p) = - 0.25_pReal*(temp33_Real(l,n)+temp33_Real(n,l))*&
if (any(xi(:) /= 0.0_pReal)) then (xidyad(m,p) +xidyad(p,m))
do l = 1,3; do m = 1,3 enddo; enddo; enddo; enddo
xinormdyad(l,m) = xi(l)*xi(m)/sum(xi**2) do m = 1,3; do n = 1,3
enddo; enddo temp33_Complex(m,n) = sum(gamma_hat(1,1,1,m,n,:,:) *(workfft(i*2-1,j,k,:,:)&
temp33_Real = math_inv3x3(math_mul3333xx33(c0, xinormdyad)) +workfft(i*2 ,j,k,:,:)*img))
else enddo; enddo
xinormdyad = 0.0_pReal workfft(i*2-1,j,k,:,:) = real (temp33_Complex) ! change of av strain
temp33_Real = 0.0_pReal workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex)
endif
do l=1,3; do m=1,3; do n=1,3; do p=1,3
gamma_hat(1,1,1, l,m,n,p) = - 0.25_pReal*(temp33_Real(l,n)+temp33_Real(n,l))*&
(xinormdyad(m,p) +xinormdyad(p,m))
enddo; enddo; enddo; enddo
do m = 1,3; do n = 1,3
temp33_Complex(m,n) = sum(gamma_hat(1,1,1,m,n,:,:) *(workfft(i*2-1,j,k,:,:)&
+workfft(i*2 ,j,k,:,:)*img))
enddo; enddo
workfft(i*2-1,j,k,:,:) = real (temp33_Complex) ! change of strain
workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex)
enddo; enddo; enddo enddo; enddo; enddo
else ! use precalculated gamma-operator else ! use precalculated gamma-operator
do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1 do k = 1, resolution(3); do j = 1, resolution(2); do i = 1, resolution(1)/2+1
@ -585,7 +578,7 @@ program DAMASK_spectral
temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n,:,:) *(workfft(i*2-1,j,k,:,:)& temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n,:,:) *(workfft(i*2-1,j,k,:,:)&
+ workfft(i*2 ,j,k,:,:)*img)) + workfft(i*2 ,j,k,:,:)*img))
enddo; enddo enddo; enddo
workfft(i*2-1,j,k,:,:) = real (temp33_Complex) ! change of strain workfft(i*2-1,j,k,:,:) = real (temp33_Complex) ! change of av strain
workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex) workfft(i*2 ,j,k,:,:) = aimag(temp33_Complex)
enddo; enddo; enddo enddo; enddo; enddo
endif endif