cleaned and restructured output on screen, changed L_inf norm to L_2 for calculation of error in Fourier space.

removed unnecessary function from math.f90
added some documents/manuals regarding spectral method
This commit is contained in:
Martin Diehl 2011-08-31 14:37:01 +00:00
parent f527bdcc78
commit 0469d37fc3
9 changed files with 743 additions and 60 deletions

View File

@ -92,7 +92,7 @@ program DAMASK_spectral
real(pReal), dimension(3,3) :: ones, zeroes, temp33_Real, & real(pReal), dimension(3,3) :: ones, zeroes, temp33_Real, &
pstress, pstress_av, defgrad_av,& pstress, pstress_av, defgrad_av,&
defgradAim, defgradAimOld, defgradAimCorr,& defgradAim, defgradAimOld, defgradAimCorr,&
mask_stress, mask_defgrad, deltaF mask_stress, mask_defgrad, fDot
real(pReal), dimension(3,3,3,3) :: dPdF, c0_reference, c_current, s_prev, c_prev ! stiffness and compliance real(pReal), dimension(3,3,3,3) :: dPdF, c0_reference, c_current, s_prev, c_prev ! stiffness and compliance
real(pReal), dimension(6) :: cstress ! cauchy stress real(pReal), dimension(6) :: cstress ! cauchy stress
real(pReal), dimension(6,6) :: dsde ! small strain stiffness real(pReal), dimension(6,6) :: dsde ! small strain stiffness
@ -354,6 +354,7 @@ program DAMASK_spectral
enddo; enddo; enddo enddo; enddo; enddo
c0_reference = c_current * wgt ! linear reference material stiffness c0_reference = c_current * wgt ! linear reference material stiffness
c_prev = c0_reference c_prev = c0_reference
do k = 1, resolution(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) do k = 1, resolution(3) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
k_s(3) = k-1 k_s(3) = k-1
if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3) if(k > resolution(3)/2+1) k_s(3) = k_s(3)-resolution(3)
@ -452,7 +453,7 @@ program DAMASK_spectral
allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal
timeinc = bc_timeIncrement(loadcase)/bc_steps(loadcase) ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used timeinc = bc_timeIncrement(loadcase)/bc_steps(loadcase) ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used
deltaF = bc_deformation(:,:,loadcase) ! only valid for given fDot. will be overwritten later in case L is given fDot = bc_deformation(:,:,loadcase) ! only valid for given fDot. will be overwritten later in case L is given
!************************************************************* !*************************************************************
! 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)
@ -471,24 +472,24 @@ program DAMASK_spectral
endif endif
time = time + timeinc time = time + timeinc
if (velGradApplied(loadcase)) & ! calculate deltaF from given L and current F if (velGradApplied(loadcase)) & ! calculate fDot from given L and current F
deltaF = math_mul33x33(bc_deformation(:,:,loadcase), defgradAim) fDot = math_mul33x33(bc_deformation(:,:,loadcase), defgradAim)
!winding forward of deformation aim !winding forward of deformation aim
temp33_Real = defgradAim temp33_Real = defgradAim
defgradAim = defgradAim & defgradAim = defgradAim &
+ guessmode * mask_stress * (defgradAim - defgradAimOld) & + guessmode * mask_stress * (defgradAim - defgradAimOld) &
+ mask_defgrad * deltaF * timeinc + mask_defgrad * fDot * timeinc
defgradAimOld = temp33_Real defgradAimOld = temp33_Real
! update local deformation gradient ! update local deformation gradient
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)
temp33_Real = defgrad(i,j,k,:,:) temp33_Real = defgrad(i,j,k,:,:)
if (velGradApplied(loadcase)) & ! use velocity gradient to calculate new deformation gradient (if not guessing) if (velGradApplied(loadcase)) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
deltaF = math_mul33x33(bc_deformation(:,:,loadcase),defgradold(i,j,k,:,:)) fDot = math_mul33x33(bc_deformation(:,:,loadcase),defgradold(i,j,k,:,:))
defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:) & ! decide if guessing along former trajectory or apply homogeneous addon defgrad(i,j,k,:,:) = defgrad(i,j,k,:,:) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))& ! guessing... + guessmode * (defgrad(i,j,k,:,:) - defgradold(i,j,k,:,:))& ! guessing...
+ (1.0_pReal-guessmode) * mask_defgrad * deltaF *timeinc ! apply the prescribed value where deformation is given if not guessing + (1.0_pReal-guessmode) * mask_defgrad * fDot *timeinc ! apply the prescribed value where deformation is given if not guessing
defgradold(i,j,k,:,:) = temp33_Real defgradold(i,j,k,:,:) = temp33_Real
enddo; enddo; enddo enddo; enddo; enddo
guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase guessmode = 1.0_pReal ! keep guessing along former trajectory during same loadcase
@ -530,12 +531,17 @@ program DAMASK_spectral
(err_div > err_div_tol .or. & (err_div > err_div_tol .or. &
err_stress > err_stress_tol)) err_stress > err_stress_tol))
iter = iter + 1_pInt iter = iter + 1_pInt
if (iter == itmax) not_converged_counter = not_converged_counter + 1 print '(A)', '************************************************************'
print '(A,/)', '============================================================' print '(3(A,I5.5,tr2)A)', '**** Loadcase = ',loadcase, 'Step = ',step, 'Iteration = ',iter,'****'
print '(3(A,I5.5,tr2)/)', 'Loadcase = ',loadcase, ' Step = ',step, ' Iteration = ',iter print '(A)', '************************************************************'
workfft = 0.0_pReal ! needed because of the padding for FFTW workfft = 0.0_pReal ! needed because of the padding for FFTW
!************************************************************* !*************************************************************
print '(A,/)', '== Update Stress Field (Constitutive Evaluation P(F)) ==' do n = 1,3; do m = 1,3
defgrad_av(m,n) = sum(defgrad(:,:,:,m,n)) * wgt
enddo; enddo
print '(a,/,3(3(f12.7,x)/))', 'Deformation Gradient:',math_transpose3x3(defgrad_av)
print '(A,/)', '== Update Stress Field (Constitutive Evaluation P(F)) ======'
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 ielem = ielem + 1
@ -561,38 +567,33 @@ program DAMASK_spectral
do n = 1,3; do m = 1,3 do n = 1,3; do m = 1,3
pstress_av(m,n) = sum(workfft(1:resolution(1),:,:,m,n)) * wgt pstress_av(m,n) = sum(workfft(1:resolution(1),:,:,m,n)) * wgt
defgrad_av(m,n) = sum(defgrad(: ,:,:,m,n)) * wgt
enddo; enddo enddo; enddo
print '(a,/,3(3(f12.7,x)/))', 'Deformation Gradient:',math_transpose3x3(defgrad_av)
print '(a,/,3(3(f12.7,x)/))', 'Piola-Kirchhoff Stress / MPa: ',math_transpose3x3(pstress_av)/1.e6 print '(a,/,3(3(f12.7,x)/))', 'Piola-Kirchhoff Stress / MPa: ',math_transpose3x3(pstress_av)/1.e6
err_stress_tol = 0.0_pReal err_stress_tol = 0.0_pReal
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 * (pstress_av - bc_stress(1:3,1:3,loadcase)))) ! maximum deviaton (tensor norm not applicable) err_stress = maxval(abs(mask_stress * (pstress_av - bc_stress(1:3,1:3,loadcase)))) ! maximum deviaton (tensor norm not applicable)
do m = 1,3 err_stress_tol = maxval(abs(mask_defgrad * pstress_av)) * err_stress_tolrel ! don't use any tensor norm because the comparison should be coherent
err_stress_tol = max(err_stress_tol, sum(abs(pstress_av(m,1:3)))) ! L_inf norm of average stress err_stress_tol = err_stress_tol * err_stress_tolrel ! weighting by relative criterion
enddo print '(A,/)', '== Correcting Deformation Gradient to Fullfill BCs ========='
err_stress_tol = err_stress_tol * err_stress_tolrel ! weighting by relative criterion print '(2(a,E10.5)/)', 'Error Stress = ',err_stress, ', Tol. = ', err_stress_tol
print '(A,/)', '== Correcting Deformation Gradient to Fullfill BCs ==' defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av - bc_stress(1:3,1:3,loadcase)))) ! residual on given stress components
print '(2(a,E10.5)/)', 'Error Stress = ',err_stress, ' Tol. = ', err_stress_tol
defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av - bc_stress(1:3,1:3,loadcase)))) ! residual on given stress components
defgradAim = defgradAim + defgradAimCorr defgradAim = defgradAim + defgradAimCorr
print '(a,/,3(3(f12.7,x)/))', 'Deformation Aim: ',math_transpose3x3(defgradAim)
print '(a,x,f12.7,/)' , 'Determinant of Deformation Aim: ', math_det3x3(defgradAim)
endif endif
print '(A,/)', '== Calculating Equilibrium Using Spectral Method ==' print '(A,/)', '== Calculating Equilibrium Using Spectral Method ==========='
err_div = 0.0_pReal
p_hat_avg = 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_inf norm of average stress in fourier space
p_hat_avg = max(p_hat_avg,sum(abs(workfft(1,1,1,m,1:3)))) ! ignore imaginary part as it is always zero for real only input)) p_hat_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(workfft(1,1,1,1:3,1:3),& ! L_2 norm of average stress in fourier space,
enddo math_transpose3x3(workfft(1,1,1,1:3,1:3)))))) ! ignore imaginary part as it is always zero for real only input))
err_div = 0.0_pReal
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
err_div = err_div + maxval(abs(math_mul33x3_complex(workfft(i*2-1,j,k,1:3,1:3)+& ! avg of L_inf norm of div(stress) in fourier space (Suquet small strain) err_div = err_div + sqrt(sum((math_mul33x3_complex(workfft(i*2-1,j,k,1:3,1:3)+& ! avg of L_2 norm of div(stress) in fourier space (Suquet small strain)
workfft(i*2, j,k,1:3,1:3)*img,xi(1:3,i,j,k)))) workfft(i*2, j,k,1:3,1:3)*img,xi(1:3,i,j,k)))**2.0))
enddo; enddo; enddo enddo; enddo; enddo
err_div = err_div*wgt/p_hat_avg*(minval(geomdimension)*wgt**(-1/4)) ! weigthting, multiplying by minimum dimension to get rid of dimension dependency err_div = err_div*wgt/p_hat_avg*(minval(geomdimension)*wgt**(-1/4)) ! weigthting, multiplying by minimum dimension to get rid of dimension dependency and phenomenologigal factor wgt**(-1/4) to get rid of resolution dependency
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 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
@ -638,28 +639,26 @@ program DAMASK_spectral
defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim(m,n) - defgrad_av(m,n)) ! anticipated target minus current state defgrad(:,:,:,m,n) = defgrad(:,:,:,m,n) + (defgradAim(m,n) - defgrad_av(m,n)) ! anticipated target minus current state
enddo; enddo enddo; enddo
do m = 1,3; do n = 1,3 print '(2(a,E10.5)/)', 'Error Divergence = ',err_div, ', Tol. = ', err_div_tol
defgrad_av(m,n) = sum(defgrad(:,:,:,m,n))*wgt
enddo; enddo
print '(a,/,3(3(f12.7,x)/))', 'Deformation Gradient:',math_transpose3x3(defgrad_av)
print '(2(a,E10.5)/)', 'Error Divergence = ',err_div, ' Tol. = ', err_div_tol
enddo ! end looping when convergency is achieved enddo ! end looping when convergency is achieved
c_prev = c_current*wgt ! calculate stiffness for next step c_prev = c_current*wgt ! calculate stiffness for next step
if (mod(step,bc_frequency(loadcase)) == 0_pInt) & ! at output frequency if (mod(step,bc_frequency(loadcase)) == 0_pInt) & ! at output frequency
write(538) materialpoint_results(:,1,:) ! write result to file write(538) materialpoint_results(:,1,:) ! write result to file
print '(2(A,I5.5),A,/)', 'Step = ',step, ' of Loadcase = ',loadcase, ' Converged' if(err_div<err_div_tol .and. err_stress<err_stress_tol) then
print '(a,x,f12.7,/)' , 'Determinant of Deformation Aim: ', math_det3x3(defgradAim) print '(2(A,I5.5),A,/)', '== Step = ',step, ' of Loadcase = ',loadcase, ' Converged =============='
print '(a,/,3(3(f12.7,x)/))', 'Deformation Aim: ',math_transpose3x3(defgradAim) else
print '(a,/,3(3(f12.7,x)/))', 'Deformation Gradient:',math_transpose3x3(defgrad_av) print '(2(A,I5.5),A,/)', '== Step = ',step, ' of Loadcase = ',loadcase, ' NOT Converged =========='
print '(a,/,3(3(f10.4,x)/))', 'Piola-Kirchhoff Stress / MPa (prev. Iteration): ',math_transpose3x3(pstress_av)/1.e6 not_converged_counter = not_converged_counter + 1
print '(A,/)', '************************************************************' endif
enddo ! end looping over steps in current loadcase enddo ! end looping over steps in current loadcase
deallocate(c_reduced) deallocate(c_reduced)
deallocate(s_reduced) deallocate(s_reduced)
enddo ! end looping over loadcases enddo ! end looping over loadcases
print '(a,i10,a)', 'A Total of ', not_converged_counter, ' Steps did not Converge!' print '(A,/)', '############################################################'
print '(a,i5.5,a)', 'A Total of ', not_converged_counter, ' Steps did not Converge!'
close(538) close(538)
call dfftw_destroy_plan(plan_fft(1)); call dfftw_destroy_plan(plan_fft(2)) call dfftw_destroy_plan(plan_fft(1)); call dfftw_destroy_plan(plan_fft(2))

View File

@ -1166,21 +1166,6 @@ pure function math_transpose3x3(A)
endfunction math_Plain33to9 endfunction math_Plain33to9
!********************************************************************
! convert 3x3 matrix into vector 9x1
!********************************************************************
pure function math_Plain33to9_logical(m33)
use prec, only: pReal,pInt
implicit none
logical, dimension(3,3), intent(in) :: m33
logical, dimension(9) :: math_Plain33to9_logical
integer(pInt) i
forall (i=1:9) math_Plain33to9_logical(i) = m33(mapPlain(1,i),mapPlain(2,i))
endfunction math_Plain33to9_logical
!******************************************************************** !********************************************************************
! convert Plain 9x1 back to 3x3 matrix ! convert Plain 9x1 back to 3x3 matrix

Binary file not shown.

View File

@ -0,0 +1,226 @@
\documentclass[12pt,numbers,sort&compress]{article}
%% Use the option review to obtain double line spacing
%% \documentclass[authoryear,preprint,review,12pt]{elsarticle}
%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}
%% if you use PostScript figures in your article
%% use the graphics package for simple commands
%% \usepackage{graphics}
%% or use the graphicx package for more complicated commands
%% \usepackage{graphicx}
%% or use the epsfig package if you prefer to use the old commands
%% \usepackage{epsfig}
%% The amssymb package provides various useful mathematical symbols
\usepackage[usenames,dvipsnames,pdftex]{color}
\usepackage{amsmath,amssymb,amsfonts}
\usepackage{siunitx}
%\usepackage{subeqnarray}
\usepackage[hang]{subfigure}
\usepackage{verbatim}
\usepackage{bm}
\usepackage{tikz}
\usetikzlibrary{arrows}
\usepackage{booktabs}
\usepackage{graphicx}
\newcommand{\pathToFigures}{./figures}
\graphicspath{{\pathToFigures/}}
\DeclareGraphicsExtensions{.pdf,.png}
\usepackage[pdftex, % hyper-references for pdftex
bookmarksnumbered=true,% % generate bookmarks with numbers
pagebackref=true,% % generate backref in biblio
colorlinks=true,%
]{hyperref}%
%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}
%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers.
%% \usepackage{lineno}
\newlength{\diagramsize}
\setlength{\diagramsize}{0.4\textwidth}
\newcommand{\question}[1]{\textcolor{Red}{#1}}
\newcommand{\note}[1]{\textcolor{CornflowerBlue}{#1}}
\newcommand{\term}[1]{\textsc{#1}}
\newcommand{\eref}[1]{Eq.~\eqref{#1}}
\newcommand{\Eref}[1]{Eq.~\eqref{#1}}
\newcommand{\erefs}[1]{Eqs.~\eqref{#1}}
\newcommand{\Erefs}[1]{Eqs.~\eqref{#1}}
\newcommand{\fref}[1]{Fig.~\ref{#1}}
\newcommand{\Fref}[1]{Fig.~\ref{#1}}
\newcommand{\frefs}[1]{Figs.~\ref{#1}}
\newcommand{\Frefs}[1]{Figs.~\ref{#1}}
\newcommand{\tref}[1]{Tab.~\ref{#1}}
\newcommand{\Tref}[1]{Tab.~\ref{#1}}
\newcommand{\trefs}[1]{Tabs.~\ref{#1}}
\newcommand{\Trefs}[1]{Tabs.~\ref{#1}}
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\cf}{\textit{cf.}}
\newcommand{\Euler}{\textsc{Euler}}
\newcommand{\Gauss}{\textsc{Gauss}}
\newcommand{\kB}{\ensuremath{k_\text{B}}}
\newcommand{\transpose}[1]{\ensuremath{{#1}^{\mathrm T}}}
\newcommand{\inverse}[1]{\ensuremath{{#1}^{-1}}}
\newcommand{\invtranspose}[1]{\ensuremath{{#1}^{\mathrm{-T}}}}
\newcommand{\sign}[1]{\ensuremath{\operatorname{sgn}\left({#1}\right)}}
\newcommand{\grad}[1][]{\ensuremath{\operatorname{grad}{#1}}}
\newcommand{\Grad}[1][]{\ensuremath{\operatorname{Grad}{#1}}}
\newcommand{\divergence}[1][]{\ensuremath{\operatorname{div}{#1}}}
\newcommand{\Divergence}[1][]{\ensuremath{\operatorname{Div}{#1}}}
\newcommand{\totalder}[2]{\ensuremath{\frac{\inc{#1}}{\inc{#2}}}}
\newcommand{\partialder}[2]{\ensuremath{\frac{\partial{#1}}{\partial{#2}}}}
\newcommand{\inc}[1]{\ensuremath{\text d{#1}}}
\newcommand{\abs}[1]{\ensuremath{\left|{#1}\right|}}
\newcommand{\norm}[1]{\ensuremath{\left|\left|{#1}\right|\right|}}
\newcommand{\avg}[1]{\ensuremath{\bar{#1}}}
\newcommand{\fluct}[1]{\ensuremath{\tilde{#1}}}
\newcommand{\FT}[1]{\ensuremath{\hat{#1}}}
\newcommand{\domain}[1]{\ensuremath{\mathcal{#1}}}
\newcommand{\tnsrfour}[1]{\ensuremath{\mathbb{#1}}}
\newcommand{\tnsr}[1]{\ensuremath{\bm{#1}}}
\newcommand{\vctr}[1]{\ensuremath{\bm{#1}}}
\newcommand{\eyetwo}{\ensuremath{\tnsr I}}
\newcommand{\eyefour}{\ensuremath{\tnsrfour I}}
\newcommand{\stiffness}{\ensuremath{\tnsrfour D}}
\newcommand{\refStiffness}{\ensuremath{\avg{\tnsrfour D}}}
\newcommand{\fPK}{\ensuremath{\tnsr P}}
\newcommand{\sPK}{\ensuremath{\tnsr S}}
\newcommand{\F}[1][]{\ensuremath{\tnsr F^{#1}}}
\newcommand{\Favg}{\ensuremath{\avg{\F}}}
\newcommand{\Ffluct}{\ensuremath{\fluct{\F}}}
\newcommand{\Fp}[1][]{\ensuremath{\tnsr F_\text{p}^{#1}}}
\newcommand{\Fe}[1][]{\ensuremath{\tnsr F_\text{e}^{#1}}}
\newcommand{\Lp}{\ensuremath{\tnsr L_\text{p}}}
\newcommand{\Q}[1]{\ensuremath{\tnsr Q^{(#1)}}}
\newcommand{\x}[2][]{\ensuremath{\vctr x^{(#2)}_\text{#1}}}
\newcommand{\dg}[2][]{\ensuremath{\Delta\vctr g^{(#2)}_\text{#1}}}
\newcommand{\g}[1][]{\ensuremath{\vctr g_\text{#1}}}
\newcommand{\A}[2][]{\ensuremath{A^{(#2)}_\text{#1}}}
\newcommand{\N}[2]{\ensuremath{\varrho^{(#1)}_\text{#2}}}
\newcommand{\Burgers}[1]{\ensuremath{\vctr s^{(#1)}}}
\newcommand{\n}[1]{\ensuremath{\vctr n^{(#1)}}}
\newcommand{\m}[2]{\ensuremath{\vctr m^{(#1)}_{#2}}}
\newcommand{\ld}[1]{\ensuremath{\vctr p^{(#1)}}}
\newcommand{\velocity}[2]{\ensuremath{v^{(#1)}_\text{#2}}}
\newcommand{\avgvelocity}[2]{\ensuremath{{\bar v}^{(#1)}_ \text{#2}}}
\newcommand{\flux}[2]{\ensuremath{\vctr f^{(#1)}_ \text{#2}}}
\newcommand{\averageflux}[2]{\ensuremath{\bar{\vctr f}^{(#1)}_ \text{#2}}}
\newcommand{\interfaceflux}[2]{\ensuremath{\tilde{\vctr f}^{(#1)}_ \text{#2}}}
\newcommand{\transmissivity}[1]{\ensuremath{\chi^{(#1)}}}
\newcommand{\galpha}{\ensuremath{\gamma^{(\alpha)}}}
\newcommand{\dotgalpha}{\ensuremath{\dot{\gamma}^{(\alpha)}}}
\newcommand{\taualpha}{\ensuremath{\tau^{(\alpha)}}}
\newcommand{\taualphamax}{\ensuremath{\hat\tau^{(\alpha)}}}
\newcommand{\density}[2]{\ensuremath{\varrho^{(#1)}_ \text{#2}}}
\newcommand{\densityfunc}[2]{\ensuremath{{\tilde\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\avgdensity}[2]{\ensuremath{{\bar\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\dotdensity}[2]{\ensuremath{\dot{\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\densityexcess}[2]{\ensuremath{\Delta\varrho^{(#1)}_ \text{#2}}}
\newcommand{\cs}[2][]{\ensuremath{\sigma^{(#1)}_ \text{#2}}}
%% Title, authors and addresses
%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for theassociated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for theassociated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for theassociated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}
\title{Calculation of Divergence in Fourier space}
%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{}
%% \address[label1]{}
%% \address[label2]{}
\author{M.~Diehl}
%% \linenumbers
% main text
\begin{document}
\maketitle
% ----------------------------------------------------------------------------------------------------------------------------
\section{Idea}
% ----------------------------------------------------------------------------------------------------------------------------
A field of 3x3 tensors in the form
\begin{equation}
\begin{pmatrix} P_{11} & P_{12} & P_{13}\\ P_{21} & P_{22} & P_{23} \\ P_{31} & P_{32}& P_{33} \end{pmatrix} = \tnsr{P}(\vec{x})
\end{equation}
depending on a discrete 3D vector $\vec x$ (spatial coordinates in real space)
\begin{equation}
\vec x \mathrm{~with~} x_1=1,2,...,N_1\mathrm{,~}x_2=1,2,...,N_2\mathrm{,~and~}x_3=1,2,...,N_3
\end{equation}
is transformed into Fourier space using the FFT (Subroutines provided by FFTW, 9 transforms of rank 3). The result is a 3x3 tensor field depending on the frequencies $\vec \xi$:
\begin{equation}
\hat{\tnsr P}(\vec \xi) = \mathcal{F}\left( \tnsr P(\vec x) \right)
\end{equation}
the divergence of a tensor is defined as:
\begin{equation}
\operatorname{div}\,(\tnsr P) =
\begin{pmatrix}
\frac{\partial P_{11}}{\partial x_1} +\frac{\partial P_{12}}{\partial x_2} +\frac{\partial P_{13}}{\partial x_3}\\[6pt]
\frac{\partial P_{21}}{\partial x_1} +\frac{\partial P_{22}}{\partial x_2} +\frac{\partial P_{23}}{\partial x_3}\\[6pt]
\frac{\partial P_{z1}}{\partial x_1} +\frac{\partial P_{z2}}{\partial x_2} +\frac{\partial P_{33}}{\partial x_3}
\end{pmatrix}
\end{equation}
It is a 3D vector. With
\begin{equation}\label{div}
\mathcal{F}(\frac{d^n f(x)}{dx^n}) = (2\pi i\xi)^n \hat{f}(\xi)
\end{equation}
follows:
\begin{equation}
\mathcal{F}\left(\operatorname{div}\,(\tnsr P) \right) = 2\pi i\vec \xi \times \hat{\tnsr P}(\vec \xi)
\end{equation}
Due to the order in which FFTW stores the array, the components of the 3D frequencies can be calculated by dividing a linear list $0,1,2,...N/2-1,\pm N/2,-(N/2-1),...,-2,-1$ by the sampling length in each dimension.
\section{Questions}
The inverse transform of $\mathcal{F}\left(\operatorname{div}\,(\tnsr P) \right)$ (FFTW, 3 transforms of rank 3) should give the divergence field in real space.
My problem is, that the calculated divergence has imaginary parts after the inverse transform. That seems strange to me, as the divergence can also be computed in real space, e.g. using FDM.
If a complex-to-complex transform is done, $\hat{\tnsr P} (\vec \xi)$ possesses conjugate complex symmetry.
$\operatorname{div}\,(\hat{\tnsr P})(\vec xi) \right)$ should also have conjugate complex symmetry as it would lead to a real-only divergence after the inverse transform. Surprisingly, even for the simple derivative as given in formula eq.(\ref{div}), the complex conjugate symmetry is also not preserved. The multiplication changes the sign of the whole term and not only of the imaginary part.
As the tensor $\tnsr P$ is a real-only tensor, it should be possible to use a real-to-complex transform for the tensor field and a complex-to-real transform for the back transform of the divergence field. Than one of the dimensions should run from $0$ to $N/2+1$ only. But that implies, that the divergence in Fourier space has conjugate compley symmetry.
\end{document}
\endinput

Binary file not shown.

View File

@ -0,0 +1,231 @@
\documentclass[12pt,numbers,sort&compress]{article}
%% Use the option review to obtain double line spacing
%% \documentclass[authoryear,preprint,review,12pt]{elsarticle}
%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}
%% if you use PostScript figures in your article
%% use the graphics package for simple commands
%% \usepackage{graphics}
%% or use the graphicx package for more complicated commands
%% \usepackage{graphicx}
%% or use the epsfig package if you prefer to use the old commands
%% \usepackage{epsfig}
%% The amssymb package provides various useful mathematical symbols
\usepackage[usenames,dvipsnames,pdftex]{color}
\usepackage{amsmath,amssymb,amsfonts}
\usepackage{siunitx}
%\usepackage{subeqnarray}
\usepackage[hang]{subfigure}
\usepackage{verbatim}
\usepackage{bm}
\usepackage{tikz}
\usetikzlibrary{arrows}
\usepackage{booktabs}
\usepackage{graphicx}
\newcommand{\pathToFigures}{./figures}
\graphicspath{{\pathToFigures/}}
\DeclareGraphicsExtensions{.pdf,.png}
\usepackage[pdftex, % hyper-references for pdftex
bookmarksnumbered=true,% % generate bookmarks with numbers
pagebackref=true,% % generate backref in biblio
colorlinks=true,%
]{hyperref}%
%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}
%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers.
%% \usepackage{lineno}
\newlength{\diagramsize}
\setlength{\diagramsize}{0.4\textwidth}
\newcommand{\question}[1]{\textcolor{Red}{#1}}
\newcommand{\note}[1]{\textcolor{CornflowerBlue}{#1}}
\newcommand{\term}[1]{\textsc{#1}}
\newcommand{\eref}[1]{Eq.~\eqref{#1}}
\newcommand{\Eref}[1]{Eq.~\eqref{#1}}
\newcommand{\erefs}[1]{Eqs.~\eqref{#1}}
\newcommand{\Erefs}[1]{Eqs.~\eqref{#1}}
\newcommand{\fref}[1]{Fig.~\ref{#1}}
\newcommand{\Fref}[1]{Fig.~\ref{#1}}
\newcommand{\frefs}[1]{Figs.~\ref{#1}}
\newcommand{\Frefs}[1]{Figs.~\ref{#1}}
\newcommand{\tref}[1]{Tab.~\ref{#1}}
\newcommand{\Tref}[1]{Tab.~\ref{#1}}
\newcommand{\trefs}[1]{Tabs.~\ref{#1}}
\newcommand{\Trefs}[1]{Tabs.~\ref{#1}}
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\cf}{\textit{cf.}}
\newcommand{\Euler}{\textsc{Euler}}
\newcommand{\Gauss}{\textsc{Gauss}}
\newcommand{\kB}{\ensuremath{k_\text{B}}}
\newcommand{\transpose}[1]{\ensuremath{{#1}^{\mathrm T}}}
\newcommand{\inverse}[1]{\ensuremath{{#1}^{-1}}}
\newcommand{\invtranspose}[1]{\ensuremath{{#1}^{\mathrm{-T}}}}
\newcommand{\sign}[1]{\ensuremath{\operatorname{sgn}\left({#1}\right)}}
\newcommand{\grad}[1][]{\ensuremath{\operatorname{grad}{#1}}}
\newcommand{\Grad}[1][]{\ensuremath{\operatorname{Grad}{#1}}}
\newcommand{\divergence}[1][]{\ensuremath{\operatorname{div}{#1}}}
\newcommand{\Divergence}[1][]{\ensuremath{\operatorname{Div}{#1}}}
\newcommand{\totalder}[2]{\ensuremath{\frac{\inc{#1}}{\inc{#2}}}}
\newcommand{\partialder}[2]{\ensuremath{\frac{\partial{#1}}{\partial{#2}}}}
\newcommand{\inc}[1]{\ensuremath{\text d{#1}}}
\newcommand{\abs}[1]{\ensuremath{\left|{#1}\right|}}
\newcommand{\norm}[1]{\ensuremath{\left|\left|{#1}\right|\right|}}
\newcommand{\avg}[1]{\ensuremath{\bar{#1}}}
\newcommand{\fluct}[1]{\ensuremath{\tilde{#1}}}
\newcommand{\FT}[1]{\ensuremath{\hat{#1}}}
\newcommand{\domain}[1]{\ensuremath{\mathcal{#1}}}
\newcommand{\tnsrfour}[1]{\ensuremath{\mathbb{#1}}}
\newcommand{\tnsr}[1]{\ensuremath{\bm{#1}}}
\newcommand{\vctr}[1]{\ensuremath{\bm{#1}}}
\newcommand{\eyetwo}{\ensuremath{\tnsr I}}
\newcommand{\eyefour}{\ensuremath{\tnsrfour I}}
\newcommand{\stiffness}{\ensuremath{\tnsrfour D}}
\newcommand{\refStiffness}{\ensuremath{\avg{\tnsrfour D}}}
\newcommand{\fPK}{\ensuremath{\tnsr P}}
\newcommand{\sPK}{\ensuremath{\tnsr S}}
\newcommand{\F}[1][]{\ensuremath{\tnsr F^{#1}}}
\newcommand{\Favg}{\ensuremath{\avg{\F}}}
\newcommand{\Ffluct}{\ensuremath{\fluct{\F}}}
\newcommand{\Fp}[1][]{\ensuremath{\tnsr F_\text{p}^{#1}}}
\newcommand{\Fe}[1][]{\ensuremath{\tnsr F_\text{e}^{#1}}}
\newcommand{\Lp}{\ensuremath{\tnsr L_\text{p}}}
\newcommand{\Q}[1]{\ensuremath{\tnsr Q^{(#1)}}}
\newcommand{\x}[2][]{\ensuremath{\vctr x^{(#2)}_\text{#1}}}
\newcommand{\dg}[2][]{\ensuremath{\Delta\vctr g^{(#2)}_\text{#1}}}
\newcommand{\g}[1][]{\ensuremath{\vctr g_\text{#1}}}
\newcommand{\A}[2][]{\ensuremath{A^{(#2)}_\text{#1}}}
\newcommand{\N}[2]{\ensuremath{\varrho^{(#1)}_\text{#2}}}
\newcommand{\Burgers}[1]{\ensuremath{\vctr s^{(#1)}}}
\newcommand{\n}[1]{\ensuremath{\vctr n^{(#1)}}}
\newcommand{\m}[2]{\ensuremath{\vctr m^{(#1)}_{#2}}}
\newcommand{\ld}[1]{\ensuremath{\vctr p^{(#1)}}}
\newcommand{\velocity}[2]{\ensuremath{v^{(#1)}_\text{#2}}}
\newcommand{\avgvelocity}[2]{\ensuremath{{\bar v}^{(#1)}_ \text{#2}}}
\newcommand{\flux}[2]{\ensuremath{\vctr f^{(#1)}_ \text{#2}}}
\newcommand{\averageflux}[2]{\ensuremath{\bar{\vctr f}^{(#1)}_ \text{#2}}}
\newcommand{\interfaceflux}[2]{\ensuremath{\tilde{\vctr f}^{(#1)}_ \text{#2}}}
\newcommand{\transmissivity}[1]{\ensuremath{\chi^{(#1)}}}
\newcommand{\galpha}{\ensuremath{\gamma^{(\alpha)}}}
\newcommand{\dotgalpha}{\ensuremath{\dot{\gamma}^{(\alpha)}}}
\newcommand{\taualpha}{\ensuremath{\tau^{(\alpha)}}}
\newcommand{\taualphamax}{\ensuremath{\hat\tau^{(\alpha)}}}
\newcommand{\density}[2]{\ensuremath{\varrho^{(#1)}_ \text{#2}}}
\newcommand{\densityfunc}[2]{\ensuremath{{\tilde\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\avgdensity}[2]{\ensuremath{{\bar\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\dotdensity}[2]{\ensuremath{\dot{\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\densityexcess}[2]{\ensuremath{\Delta\varrho^{(#1)}_ \text{#2}}}
\newcommand{\cs}[2][]{\ensuremath{\sigma^{(#1)}_ \text{#2}}}
%% Title, authors and addresses
%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for theassociated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for theassociated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for theassociated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}
\title{Fourier Transforms}
%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{}
%% \address[label1]{}
%% \address[label2]{}
\author{M.~Diehl}
%% \linenumbers
% main text
\begin{document}
\maketitle
% ----------------------------------------------------------------------------------------------------------------------------
\section{Discrete vs. continuous FT}
% ----------------------------------------------------------------------------------------------------------------------------
continuous Fourier transform
\begin{equation}
\hat{f}(k) = \int \limits_{-\pi}^{\pi} f(x) \cdot e^{-2\pi i k x} \inc x
\end{equation}
discrete Fourier transform
\begin{align}
\hat{f}_k &= \frac{1}{d} \sum\limits_{n=0}^{N-1} f \left( x = \frac{n}{N}d \right) \cdot e^{\left(-2 \pi i \cdot \frac{k}{d} \cdot \frac{n}{N} \cdot d \right)} \cdot \frac{d}{N}\\
&= \frac{1}{N} \sum \limits_{n=0}^{N-1} f \left( x = \frac{n}{N} d \right) \cdot e^{-\frac{2 \pi i}{N} \cdot k \cdot n}
\end{align}
% ----------------------------------------------------------------------------------------------------------------------------
\section{Differentation}
% ----------------------------------------------------------------------------------------------------------------------------
Expression in frequency and angular frequency
\begin{align}
\hat{f}(k) &= \frac{1}{d} \int \limits_0^d f(x) e^{\frac{-2 \pi i}{d} k x}\\
&= \frac{1}{d} \int \limits_0^d f(x) e^{-2 \pi i \xi x} \inc x
\end{align}
\begin{align}
\hat{f}'&= \frac{\partial}{\partial x} \left( \int \limits_{-\infty }^{\infty} \hat{f}(x) \cdot e^{i \xi x} \inc k \right)\\
&= \int \limits_{-\infty}^{\infty} \frac{\partial}{\partial x} \left( \hat{f}(x) \cdot e^{i \xi x} \right) \inc k\\
&= \int \limits_{-\infty}^{\infty} i \xi \cdot \hat{f}(x) \cdot e^{i \xi x} \inc k
\end{align}
\section{Transform}
example with $N=4$ and $x = \mathrm{sin}\left( \frac{n}{N} \cdot 2 \pi \right)$
\begin{align}
X_k &= \sum \limits_{n=0}^{N-1} x_n \cdot e ^{- \frac{2 \pi i }{N} \cdot k \cdot n};~~~k = 0;1;..;N-1\\
&= \sum \limits_{n=0}^{N-1} x_n \cdot \left(\mathrm{cos}\left(- \frac{2 \pi}{N}\cdot k \cdot n \right) + i\cdot \mathrm{sin}\left(- \frac{2 \pi}{N} \cdot k \cdot n\right) \right)\\
X_0 &= \sum \limits_{n=0}^{N-1} x_n e ^\\
&= 0 + 1 + 0 +(-1) = 0 \\
X_1 &= \sum \limits_{n=0}^{N-1} x_n e^{-i \frac{2\pi}{N} \cdot 1 \cdot n} = 0 + e ^ {-i \frac{2\pi}{N} \cdot 1 \cdot 1} + 0 +e ^ {-i \frac{2\pi}{N} \cdot 1 \cdot 3}\\
&= 0 + (- 2i)\\
X_2 &= 0 + 0i\\
X_3 &= 0 + 2i
\end{align}
$X_2$ is Nyquist frequency and has only a real part, $X_3$ is conjugate complex of $X_1$ for real only input.
\section{Inverse Transform}
\begin{align}
x_n &= \frac{1}{N} \sum \limits_{k=0}^{N-1} X_k e^{\frac{2\pi i}{N} \cdot k \cdot n}\\
x_0 &= \frac{1}{4}\left(0 - 2ie^0 + 0 + 2ie^0 \right) = 0\\
x_1 &= \frac{1}{4}\left(0 - 2ie^{\frac{2\pi i}{4}\cdot 1 \cdot 1} + 0 + 2ie^{\frac{2\pi i}{4}\cdot 3 \cdot 1} \right) = 1\\
x_2 &= 0\\
x_1 &= -1
\end{align}
\end{document}
\endinput

View File

@ -0,0 +1,202 @@
\documentclass[12pt,numbers,sort&compress]{article}
%% Use the option review to obtain double line spacing
%% \documentclass[authoryear,preprint,review,12pt]{elsarticle}
%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[final,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%% \documentclass[final,5p,times,twocolumn]{elsarticle}
%% if you use PostScript figures in your article
%% use the graphics package for simple commands
%% \usepackage{graphics}
%% or use the graphicx package for more complicated commands
%% \usepackage{graphicx}
%% or use the epsfig package if you prefer to use the old commands
%% \usepackage{epsfig}
%% The amssymb package provides various useful mathematical symbols
\usepackage[usenames,dvipsnames,pdftex]{color}
\usepackage{amsmath,amssymb,amsfonts}
\usepackage{siunitx}
%\usepackage{subeqnarray}
\usepackage[hang]{subfigure}
\usepackage{verbatim}
\usepackage{bm}
\usepackage{tikz}
\usetikzlibrary{arrows}
\usepackage{booktabs}
\usepackage[latin1]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{graphicx}
\usepackage{natbib}
\newcommand{\pathToFigures}{./figures}
\graphicspath{{\pathToFigures/}}
\DeclareGraphicsExtensions{.pdf,.png}
\usepackage[pdftex, % hyper-references for pdftex
bookmarksnumbered=true,% % generate bookmarks with numbers
%pagebackref=true,% % generate backref in biblio
%colorlinks=true,%
]{hyperref}%
%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}
%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers.
%% \usepackage{lineno}
\newlength{\diagramsize}
\setlength{\diagramsize}{0.4\textwidth}
\newcommand{\question}[1]{\textcolor{Red}{#1}}
\newcommand{\note}[1]{\textcolor{CornflowerBlue}{#1}}
\newcommand{\term}[1]{\textsc{#1}}
\newcommand{\eref}[1]{Eq.~\eqref{#1}}
\newcommand{\Eref}[1]{Eq.~\eqref{#1}}
\newcommand{\erefs}[1]{Eqs.~\eqref{#1}}
\newcommand{\Erefs}[1]{Eqs.~\eqref{#1}}
\newcommand{\fref}[1]{Fig.~\ref{#1}}
\newcommand{\Fref}[1]{Fig.~\ref{#1}}
\newcommand{\frefs}[1]{Figs.~\ref{#1}}
\newcommand{\Frefs}[1]{Figs.~\ref{#1}}
\newcommand{\tref}[1]{Tab.~\ref{#1}}
\newcommand{\Tref}[1]{Tab.~\ref{#1}}
\newcommand{\trefs}[1]{Tabs.~\ref{#1}}
\newcommand{\Trefs}[1]{Tabs.~\ref{#1}}
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\cf}{\textit{cf.}}
\newcommand{\Euler}{\textsc{Euler}}
\newcommand{\Gauss}{\textsc{Gauss}}
\newcommand{\kB}{\ensuremath{k_\text{B}}}
\newcommand{\transpose}[1]{\ensuremath{{#1}^{\mathrm T}}}
\newcommand{\inverse}[1]{\ensuremath{{#1}^{-1}}}
\newcommand{\invtranspose}[1]{\ensuremath{{#1}^{\mathrm{-T}}}}
\newcommand{\sign}[1]{\ensuremath{\operatorname{sgn}\left({#1}\right)}}
\newcommand{\grad}[1][]{\ensuremath{\operatorname{grad}{#1}}}
\newcommand{\Grad}[1][]{\ensuremath{\operatorname{Grad}{#1}}}
\newcommand{\divergence}[1][]{\ensuremath{\operatorname{div}{#1}}}
\newcommand{\Divergence}[1][]{\ensuremath{\operatorname{Div}{#1}}}
\newcommand{\totalder}[2]{\ensuremath{\frac{\inc{#1}}{\inc{#2}}}}
\newcommand{\partialder}[2]{\ensuremath{\frac{\partial{#1}}{\partial{#2}}}}
\newcommand{\inc}[1]{\ensuremath{\text d{#1}}}
\newcommand{\abs}[1]{\ensuremath{\left|{#1}\right|}}
\newcommand{\norm}[1]{\ensuremath{\left|\left|{#1}\right|\right|}}
\newcommand{\avg}[1]{\ensuremath{\bar{#1}}}
\newcommand{\fluct}[1]{\ensuremath{\tilde{#1}}}
\newcommand{\FT}[1]{\ensuremath{\hat{#1}}}
\newcommand{\domain}[1]{\ensuremath{\mathcal{#1}}}
\newcommand{\tnsrfour}[1]{\ensuremath{\mathbb{#1}}}
\newcommand{\tnsr}[1]{\ensuremath{\bm{#1}}}
\newcommand{\vctr}[1]{\ensuremath{\bm{#1}}}
\newcommand{\eyetwo}{\ensuremath{\tnsr I}}
\newcommand{\eyefour}{\ensuremath{\tnsrfour I}}
\newcommand{\stiffness}{\ensuremath{\tnsrfour D}}
\newcommand{\refStiffness}{\ensuremath{\avg{\tnsrfour D}}}
\newcommand{\fPK}{\ensuremath{\tnsr P}}
\newcommand{\sPK}{\ensuremath{\tnsr S}}
\newcommand{\F}[1][]{\ensuremath{\tnsr F^{#1}}}
\newcommand{\Favg}{\ensuremath{\avg{\F}}}
\newcommand{\Ffluct}{\ensuremath{\fluct{\F}}}
\newcommand{\Fp}[1][]{\ensuremath{\tnsr F_\text{p}^{#1}}}
\newcommand{\Fe}[1][]{\ensuremath{\tnsr F_\text{e}^{#1}}}
\newcommand{\Lp}{\ensuremath{\tnsr L_\text{p}}}
\newcommand{\Q}[1]{\ensuremath{\tnsr Q^{(#1)}}}
\newcommand{\x}[2][]{\ensuremath{\vctr x^{(#2)}_\text{#1}}}
\newcommand{\dg}[2][]{\ensuremath{\Delta\vctr g^{(#2)}_\text{#1}}}
\newcommand{\g}[1][]{\ensuremath{\vctr g_\text{#1}}}
\newcommand{\A}[2][]{\ensuremath{A^{(#2)}_\text{#1}}}
\newcommand{\N}[2]{\ensuremath{\varrho^{(#1)}_\text{#2}}}
\newcommand{\Burgers}[1]{\ensuremath{\vctr s^{(#1)}}}
\newcommand{\n}[1]{\ensuremath{\vctr n^{(#1)}}}
\newcommand{\m}[2]{\ensuremath{\vctr m^{(#1)}_{#2}}}
\newcommand{\ld}[1]{\ensuremath{\vctr p^{(#1)}}}
\newcommand{\velocity}[2]{\ensuremath{v^{(#1)}_\text{#2}}}
\newcommand{\avgvelocity}[2]{\ensuremath{{\bar v}^{(#1)}_ \text{#2}}}
\newcommand{\flux}[2]{\ensuremath{\vctr f^{(#1)}_ \text{#2}}}
\newcommand{\averageflux}[2]{\ensuremath{\bar{\vctr f}^{(#1)}_ \text{#2}}}
\newcommand{\interfaceflux}[2]{\ensuremath{\tilde{\vctr f}^{(#1)}_ \text{#2}}}
\newcommand{\transmissivity}[1]{\ensuremath{\chi^{(#1)}}}
\newcommand{\galpha}{\ensuremath{\gamma^{(\alpha)}}}
\newcommand{\dotgalpha}{\ensuremath{\dot{\gamma}^{(\alpha)}}}
\newcommand{\taualpha}{\ensuremath{\tau^{(\alpha)}}}
\newcommand{\taualphamax}{\ensuremath{\hat\tau^{(\alpha)}}}
\newcommand{\density}[2]{\ensuremath{\varrho^{(#1)}_ \text{#2}}}
\newcommand{\densityfunc}[2]{\ensuremath{{\tilde\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\avgdensity}[2]{\ensuremath{{\bar\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\dotdensity}[2]{\ensuremath{\dot{\varrho}^{(#1)}_ \text{#2}}}
\newcommand{\densityexcess}[2]{\ensuremath{\Delta\varrho^{(#1)}_ \text{#2}}}
\newcommand{\cs}[2][]{\ensuremath{\sigma^{(#1)}_ \text{#2}}}
%% Title, authors and addresses
%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for theassociated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for theassociated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for theassociated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%% \title{Title\tnoteref{label1}}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
%% \ead{email address}
%% \ead[url]{home page}
%% \fntext[label2]{}
%% \cortext[cor1]{}
%% \address{Address\fnref{label3}}
%% \fntext[label3]{}
\title{Implement a remeshing scheme into a spectral method based crystal plasticity code}
%% use optional labels to link authors explicitly to addresses:
%% \author[label1,label2]{}
%% \address[label1]{}
%% \address[label2]{}
\author{M.~Diehl}
%% \linenumbers
% main text
\begin{document}
\maketitle
% ----------------------------------------------------------------------------------------------------------------------------
At MPIE, the flexible crystal plasticity framework ``Düsseldorf Advanced MAterial Simulation Kit'' (DAMASK) is developed.
It consists of different constitutive models, homogenization schemes, and tools for post- and preprocessing \cite{Roters_etal2010}.
It has interfaces to different solvers to the mechanical boundary value problem.
To compute the boundary value problem, commercial FEM software like MSC.Marc or Abaqus or a solver based on a so-called spectral method \cite{Moulinec+Suquet1998,Lebensohn2001}.
Spectral methods have advantages concerning accuracy, performance, and memory efficiency compared to the de-facto standard FEM.
However, their use is limited to periodic boundary conditions due to the approximation of the solution by plane waves.
The spectral method implemented at MPIE uses a finite-strain formulation proposed in \cite{Lahellec_etal2001} that is written in terms of deformation gradient \tnsr F and Piola--Kirchhoff stress \tnsr P and can therefore be used to solve the mechanical boundary value problem in the reference configuration.
Calculations have shown that for inhomogeneous material convergence cannot be achieved any longer at strains larger than ca.~15--20 \%.
We presently believe that this is due to the fact that the regular mesh in the reference configuration is locally heavily deformed to an extent where single points cross the path of neighboring points.
To reach higher strains, a remeshing scheme should be implemented as follows.\begin{enumerate}
\item Write out the current state
\item Approximate the deformed configuration by a regular (undeformed, new) mesh
\item Translate the old state values to the new mesh
\end{enumerate}
\bibliographystyle{unsrtnat}
\bibliography{Masterthesis}
\end{document}
\endinput

View File

@ -0,0 +1,40 @@
To use the spectral method, a Fortran compiler must be installed.
For post an prepocessing, a running python enviroment with numpy must be installed.
Currently the Intel Fortran compiler and the GNU Fortran compiler are supported.
The compilation of the code is automated by using make.
There are several options to tune the compilation in detail
The most important switch is "F90" to choose the Fortran compiler.
It can be set to "ifort" (Intel Fortran compiler) or "gfortran" (GNU compiler).
E.g. "make F90=gfortran"
Multiprocessor support can be switched off by setting "OPENMP=OFF" (default "OPENMP=ON").
The compiler switches are explained in detail in the makefile in %DAMASK%/code
A "make clean" removes all previously compiled code.
If the compilation was succesfull, a file called DAMASK_spectral.exe is found in %DAMASK%/code.
If multiprocessor support is enabled (default), the number of threads can be set to N by setting "DAMASK_NUM_THREADS".
E.g. "export DAMASK_NUM_THREADS=4" for using four cores.
The number of processors can be shown by "echo $DAMASK_NUM_THREADS"
To run the application in the directory containing loadcase, geometry and *.config files, type %DAMASK%/code/DAMASK_spectral.exe geomfile.geom loadcasefile.load
To create geomtry file, use the tools for voronoi tesselation provided in %DAMASK%/processing/pre
At first, install all processing tools by running %DAMASK%/processing/setup/setup_processing.py
!Be carefull, change "ifort" to your compiler name in "os.system('ifort -O3 -parallel -o%s.exe %s'%(os.path.splitext(src)[0],src))" in that file
!And change f90flags="-heap-arrays 500000000" to f90flags="" in "%DAMASK%/processing/post/make_postprocessingMath" for non-Intel compiler
If the post- and preprocessing are correctly installed, you have the file "voronoi_randomSeeding.exe" and "voronoi_tessellation.exe" in "%DAMASK%/processing/pre"
To generate random seeds, run "voronoi_randomSeeding.exe" an follow the instructions on the screen.
You will get a file called *.seeds in your working directory
With the *.seeds file, you can generate geometry files at different resolution and dimensionby running "voronoi_randomSeeding.exe" in "%DAMASK%/processing/pre"
You will get three files called *.spectral, *.geom, *.config .
Ingnore the first one.
Rename the file *.config to "material.config".
This file conatines the information about the orientation of each grain in the generated microstructure.
It DOES NOT contain any information about the constitutive model that should be used.
Copy that information from %DAMASK%/code/config/material.config
A valid "material.config" contains at least one <homogenization> scheme, a <microstructure> for each grain, at leat one <texture>, a <crystallite> and at leat one <phase>