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

allowed arbitrary grid (formerly know as resolution), i.e. any integer>1 in each dimension is possible. 

removed square root correction for divergence (doesn't make sense)
fixed FFT debug of spectral solver
This commit is contained in:
Martin Diehl 2013-04-10 10:19:16 +00:00
parent 4338cd13bc
commit ddcc795461
6 changed files with 227 additions and 197 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
code/DAMASK_spectral_solverAL.f90
View File

@ -527,8 +527,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc
err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual

8
code/DAMASK_spectral_solverBasic.f90
View File

@ -228,7 +228,7 @@ type(tSolutionState) function &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! loop variables, convergence etc.
real(pReal) :: err_div, err_stress real(pReal) :: err_div, err_stress
integer(pInt) :: iter, row, column integer(pInt) :: iter
logical :: ForwardData logical :: ForwardData
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -375,9 +375,9 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress)
err_stress/err_stress_tol ] < 1.0_pReal) err_stress/err_stress_tol ] < 1.0_pReal)
write(6,'(/,a,f10.2,a,es11.5,a,es11.4,a)') ' error divergence = ', & write(6,'(/,a,f10.2,a,es11.5,a,es11.4,a)') ' error divergence = ', &
err_div/pAvgDivL2/err_div_tol, ' (',err_div/pAvgDivL2,' / m, tol =',err_div_tol,')' err_div/pAvgDivL2/err_div_tol, ' (',err_div/pAvgDivL2,' / m, tol =',err_div_tol,')'
write(6,'(a,f8.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', & write(6, '(a,f10.2,a,es11.5,a,es11.4,a)') ' error stress BC = ', &
err_stress/err_stress_tol, ' (',err_stress, ' Pa , tol =',err_stress_tol,')' err_stress/err_stress_tol, ' (',err_stress, ' Pa, tol =',err_stress_tol,')'
flush(6) flush(6)
end function basic_Converged end function basic_Converged

225
code/DAMASK_spectral_utilities.f90
View File

@ -52,14 +52,17 @@ module DAMASK_spectral_utilities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! debug divergence ! debug divergence
real(pReal), private, dimension(:,:,:,:), pointer :: divergence_real !< scalar field real representation for debugging divergence calculation real(pReal), private, dimension(:,:,:,:), pointer :: divReal !< scalar field real representation for debugging divergence calculation
complex(pReal),private, dimension(:,:,:,:), pointer :: divergence_fourier !< scalar field real representation for debugging divergence calculation complex(pReal),private, dimension(:,:,:,:), pointer :: divFourier !< scalar field real representation for debugging divergence calculation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! plans for FFTW ! plans for FFTW
type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back !< plans for FFTW in case of debugging the Fourier transform type(C_PTR), private :: &
type(C_PTR), private :: plan_forward, plan_backward !< plans for FFTW planForth, & !< FFTW plan P(x) to P(k)
type(C_PTR), private :: plan_divergence !< plan for FFTW in case of debugging divergence calculation planBack, & !< FFTW plan F(k) to F(x)
planDebugForth, & !< FFTW plan for scalar field (proof that order of usual transform is correct)
planDebugBack, & !< FFTW plan for scalar field inverse (proof that order of usual transform is correct)
planDiv !< plan for FFTW in case of debugging divergence calculation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables controlling debugging ! variables controlling debugging
@ -73,14 +76,14 @@ module DAMASK_spectral_utilities
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type, public :: tSolutionState !< return type of solution from spectral solver variants type, public :: tSolutionState !< return type of solution from spectral solver variants
logical :: converged = .true. logical :: converged = .true.
logical :: regrid = .false. logical :: regrid = .false.
logical :: termIll = .false. logical :: termIll = .false.
integer(pInt) :: iterationsNeeded = 0_pInt integer(pInt) :: iterationsNeeded = 0_pInt
end type tSolutionState end type tSolutionState
type, public :: tBoundaryCondition !< set of parameters defining a boundary condition type, public :: tBoundaryCondition !< set of parameters defining a boundary condition
real(pReal), dimension(3,3) :: values = 0.0_pReal real(pReal), dimension(3,3) :: values = 0.0_pReal
real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal
logical, dimension(3,3) :: maskLogical = .false. logical, dimension(3,3) :: maskLogical = .false.
@ -159,7 +162,7 @@ subroutine utilities_init()
tensorField, & !< field cotaining data for FFTW in real and fourier space (in place) tensorField, & !< field cotaining data for FFTW in real and fourier space (in place)
scalarField_realC, & !< field cotaining data for FFTW in real space when debugging FFTW (no in place) scalarField_realC, & !< field cotaining data for FFTW in real space when debugging FFTW (no in place)
scalarField_fourierC, & !< field cotaining data for FFTW in fourier space when debugging FFTW (no in place) scalarField_fourierC, & !< field cotaining data for FFTW in fourier space when debugging FFTW (no in place)
divergence !< field cotaining data for FFTW in real and fourier space when debugging divergence (in place) div !< field cotaining data for FFTW in real and fourier space when debugging divergence (in place)
write(6,'(/,a)') ' <<<+- DAMASK_spectral_utilities init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_utilities init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -188,8 +191,8 @@ subroutine utilities_init()
! allocation ! allocation
allocate (xi(3,res1_red,res(2),res(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension allocate (xi(3,res1_red,res(2),res(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension
tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8) tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8)
call c_f_pointer(tensorField, field_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField call c_f_pointer(tensorField, field_real, [res(1)+2_pInt-mod(res(1),2_pInt),res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField
call c_f_pointer(tensorField, field_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField call c_f_pointer(tensorField, field_fourier,[res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details) ! general initialization of FFTW (see manual on fftw.org for more details)
@ -203,29 +206,31 @@ subroutine utilities_init()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! creating plans for the convolution ! creating plans for the convolution
plan_forward = fftw_plan_many_dft_r2c(3, [res(3),res(2) ,res(1)], 9,& ! dimensions, logical length in each dimension in reversed order, no. of transforms planForth = fftw_plan_many_dft_r2c(3,[res(3),res(2) ,res(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms
field_real, [res(3),res(2) ,res(1)+2_pInt],& ! input data, physical length in each dimension in reversed order field_real,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & ! input data, physical length in each dimension in reversed order
1, res(3)*res(2)*(res(1)+2_pInt),& ! striding, product of physical length in the 3 dimensions 1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
field_fourier, [res(3),res(2) ,res1_red],& ! output data, physical length in each dimension in reversed order field_fourier,[res(3),res(2) ,res1_red], & ! output data, physical length in each dimension in reversed order
1, res(3)*res(2)* res1_red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision 1, res(3)*res(2)* res1_red, fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision
plan_backward = fftw_plan_many_dft_c2r(3, [res(3),res(2) ,res(1)], 9,& ! dimensions, logical length in each dimension in reversed order, no. of transforms planBack = fftw_plan_many_dft_c2r(3,[res(3),res(2) ,res(1)], 9, & ! dimensions, logical length in each dimension in reversed order, no. of transforms
field_fourier, [res(3),res(2) ,res1_red],& ! input data, physical length in each dimension in reversed order field_fourier,[res(3),res(2) ,res1_red], & ! input data, physical length in each dimension in reversed order
1, res(3)*res(2)* res1_red,& ! striding, product of physical length in the 3 dimensions 1, res(3)*res(2)* res1_red, & ! striding, product of physical length in the 3 dimensions
field_real, [res(3),res(2) ,res(1)+2_pInt],& ! output data, physical length in each dimension in reversed order field_real,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], & ! output data, physical length in each dimension in reversed order
1, res(3)*res(2)*(res(1)+2_pInt), fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision 1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
fftw_planner_flag) ! planner precision
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! depending on debug options, allocate more memory and create additional plans ! depending on debug options, allocate more memory and create additional plans
if (debugDivergence) then if (debugDivergence) then
divergence = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T)) div = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T))
call c_f_pointer(divergence, divergence_real, [ res(1)+2_pInt,res(2),res(3),3]) call c_f_pointer(div,divReal, [res(1)+2_pInt-mod(res(1),2_pInt),res(2),res(3),3])
call c_f_pointer(divergence, divergence_fourier, [ res1_red, res(2),res(3),3]) call c_f_pointer(div,divFourier,[res1_red, res(2),res(3),3])
plan_divergence = fftw_plan_many_dft_c2r(3,[ res(3),res(2) ,res(1)],3,& planDiv = fftw_plan_many_dft_c2r(3,[res(3),res(2) ,res(1)],3,&
divergence_fourier,[ res(3),res(2) ,res1_red],& divFourier,[res(3),res(2) ,res1_red],&
1, res(3)*res(2)* res1_red,& 1, res(3)*res(2)* res1_red,&
divergence_real,[ res(3),res(2) ,res(1)+2_pInt],& divReal,[res(3),res(2) ,res(1)+2_pInt-mod(res(1),2_pInt)], &
1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) 1, res(3)*res(2)*(res(1)+2_pInt-mod(res(1),2_pInt)), &
fftw_planner_flag)
endif endif
if (debugFFTW) then if (debugFFTW) then
@ -233,9 +238,9 @@ subroutine utilities_init()
scalarField_fourierC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! allocate data for fourier representation (no in place transform) scalarField_fourierC = fftw_alloc_complex(int(res(1)*res(2)*res(3),C_SIZE_T)) ! allocate data for fourier representation (no in place transform)
call c_f_pointer(scalarField_realC, scalarField_real, [res(1),res(2),res(3)]) ! place a pointer for a real representation call c_f_pointer(scalarField_realC, scalarField_real, [res(1),res(2),res(3)]) ! place a pointer for a real representation
call c_f_pointer(scalarField_fourierC, scalarField_fourier, [res(1),res(2),res(3)]) ! place a pointer for a fourier representation call c_f_pointer(scalarField_fourierC, scalarField_fourier, [res(1),res(2),res(3)]) ! place a pointer for a fourier representation
plan_scalarField_forth = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) planDebugForth = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style)
scalarField_real,scalarField_fourier,-1,fftw_planner_flag) ! input, output, forward FFT(-1), planner precision scalarField_real,scalarField_fourier,-1,fftw_planner_flag) ! input, output, forward FFT(-1), planner precision
plan_scalarField_back = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style) planDebugBack = fftw_plan_dft_3d(res(3),res(2),res(1),& ! reversed order (C style)
scalarField_fourier,scalarField_real,+1,fftw_planner_flag) ! input, output, backward (1), planner precision scalarField_fourier,scalarField_real,+1,fftw_planner_flag) ! input, output, backward (1), planner precision
endif endif
@ -326,7 +331,7 @@ end subroutine utilities_updateGamma
!> In case of debugging the FFT, also one component of the tensor (specified by row and column) !> In case of debugging the FFT, also one component of the tensor (specified by row and column)
!> is independetly transformed complex to complex and compared to the whole tensor transform !> is independetly transformed complex to complex and compared to the whole tensor transform
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_FFTforward(row,column) subroutine utilities_FFTforward() !< @ToDo make row and column between randomly between 1 and 3
use math use math
use mesh, only : & use mesh, only : &
scaledDim, & scaledDim, &
@ -334,42 +339,55 @@ subroutine utilities_FFTforward(row,column)
res1_red res1_red
implicit none implicit none
integer(pInt), intent(in), optional :: row, column !< if debug FFTW, compare 3D array field of row and column integer(pInt) :: row, column ! if debug FFTW, compare 3D array field of row and column
integer(pInt), dimension(2:3,2) :: Nyquist ! highest frequencies to be removed (1 if even, 2 if odd)
real(pReal), dimension(2) :: myRand ! random numbers
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! copy one component of the stress field to to a single FT and check for mismatch ! copy one component of the stress field to to a single FT and check for mismatch
if (debugFFTW .and. present(row) .and. present(column)) & if (debugFFTW) then
scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component call random_number(myRand) ! two numbers: 0 <= x < 1
cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) row = nint(myRand(1)*2_pReal + 1_pReal,pInt)
column = nint(myRand(2)*2_pReal + 1_pReal,pInt)
scalarField_real = cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) ! store the selected component
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing the FFT ! doing the FFT
call fftw_execute_dft_r2c(plan_forward,field_real,field_fourier) call fftw_execute_dft_r2c(planForth,field_real,field_fourier)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 FT results ! comparing 1 and 3x3 FT results
if (debugFFTW .and. present(row) .and. present(column)) then if (debugFFTW) then
call fftw_execute_dft(plan_scalarField_forth,scalarField_real,scalarField_fourier) call fftw_execute_dft(planDebugForth,scalarField_real,scalarField_fourier)
write(6,'(/,a,i1,1x,i1,a)') ' .. checking FT results of compontent ', row, column, ' ..' where(abs(scalarField_fourier(1:res1_red,1:res(2),1:res(3))) > tiny(1.0_pReal)) ! avoid division by zero
flush(6) scalarField_fourier(1:res1_red,1:res(2),1:res(3)) = &
write(6,'(/,a,2(es11.4,1x))') ' max FT relative error = ',& ! print real and imaginary part seperately (scalarField_fourier(1:res1_red,1:res(2),1:res(3))-&
maxval( real((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/&
field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& scalarField_fourier(1:res1_red,1:res(2),1:res(3))
scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), & else where
maxval(aimag((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& scalarField_real = cmplx(0.0,0.0,pReal)
field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& end where
scalarField_fourier(1:res1_red,1:res(2),1:res(3)))) write(6,'(/,a,i1,1x,i1,a)') ' .. checking FT results of compontent ', row, column, ' ..'
flush(6) write(6,'(/,a,2(es11.4,1x))') ' max FT relative error = ',& ! print real and imaginary part seperately
endif maxval(real (scalarField_fourier(1:res1_red,1:res(2),1:res(3)))),&
maxval(aimag(scalarField_fourier(1:res1_red,1:res(2),1:res(3))))
flush(6)
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! removing highest frequencies ! removing highest frequencies
field_fourier ( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& Nyquist(2,1:2) = [res(2)/2_pInt + 1_pInt, res(2)/2_pInt + 1_pInt + mod(res(2),2_pInt)]
Nyquist(3,1:2) = [res(3)/2_pInt + 1_pInt, res(3)/2_pInt + 1_pInt + mod(res(3),2_pInt)]
if(res(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
field_fourier (res1_red, 1:res(2), 1:res(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal) = cmplx(0.0_pReal,0.0_pReal,pReal)
field_fourier (1:res1_red, res(2)/2_pInt+1_pInt,1:res(3) ,1:3,1:3)& if(res(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
field_fourier (1:res1_red,Nyquist(2,1):Nyquist(2,2),1:res(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal) = cmplx(0.0_pReal,0.0_pReal,pReal)
if(res(3)>1_pInt) & ! do not delete the whole slice in case of 2D calculation if(res(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
field_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& field_fourier (1:res1_red,1:res(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal) = cmplx(0.0_pReal,0.0_pReal,pReal)
end subroutine utilities_FFTforward end subroutine utilities_FFTforward
@ -381,7 +399,7 @@ end subroutine utilities_FFTforward
!> is independetly transformed complex to complex and compared to the whole tensor transform !> is independetly transformed complex to complex and compared to the whole tensor transform
!> results is weighted by number of points stored in wgt !> results is weighted by number of points stored in wgt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_FFTbackward(row,column) subroutine utilities_FFTbackward()
use math !< must use the whole module for use of FFTW use math !< must use the whole module for use of FFTW
use mesh, only: & use mesh, only: &
wgt, & wgt, &
@ -389,14 +407,19 @@ subroutine utilities_FFTbackward(row,column)
res1_red res1_red
implicit none implicit none
integer(pInt), intent(in), optional :: row, column !< if debug FFTW, compare 3D array field of row and column integer(pInt) :: row, column !< if debug FFTW, compare 3D array field of row and column
integer(pInt) :: i, j, k, m, n integer(pInt) :: i, j, k, m, n
real(pReal), dimension(2) :: myRand
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! unpack FFT data for conj complex symmetric part. This data is not transformed when using c2r ! unpack FFT data for conj complex symmetric part. This data is not transformed when using c2r
if (debugFFTW .and. present(row) .and. present(column)) then if (debugFFTW) then
scalarField_fourier = field_fourier(1:res1_red,1:res(2),1:res(3),row,column) call random_number(myRand) ! two numbers: 0 <= x < 1
do i = 0_pInt, res(1)/2_pInt-2_pInt row = nint(myRand(1)*2_pReal + 1_pReal,pInt)
column = nint(myRand(2)*2_pReal + 1_pReal,pInt)
scalarField_fourier(1:res1_red,1:res(2),1:res(3)) &
= field_fourier(1:res1_red,1:res(2),1:res(3),row,column)
do i = 0_pInt, res(1)/2_pInt-2_pInt + mod(res(1),2_pInt)
m = 1_pInt m = 1_pInt
do k = 1_pInt, res(3) do k = 1_pInt, res(3)
n = 1_pInt n = 1_pInt
@ -412,22 +435,25 @@ subroutine utilities_FFTbackward(row,column)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing the iFFT ! doing the iFFT
call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient call fftw_execute_dft_c2r(planBack,field_fourier,field_real) ! back transform of fluct deformation gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! comparing 1 and 3x3 inverse FT results ! comparing 1 and 3x3 inverse FT results
if (debugFFTW .and. present(row) .and. present(column)) then if (debugFFTW) then
write(6,'(/,a,i1,1x,i1,a)') ' ... checking iFT results of compontent ', row, column, ' ..' call fftw_execute_dft(planDebugBack,scalarField_fourier,scalarField_real)
flush(6) where(abs(real(scalarField_real,pReal)) > tiny(1.0_pReal)) ! avoid division by zero
call fftw_execute_dft(plan_scalarField_back,scalarField_fourier,scalarField_real) scalarField_real = (scalarField_real &
write(6,'(/,a,es11.4)') ' max iFT relative error = ',& - cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column), 0.0, pReal))/ &
maxval((real(scalarField_real(1:res(1),1:res(2),1:res(3)))-& scalarField_real
field_real(1:res(1),1:res(2),1:res(3),row,column))/& else where
real(scalarField_real(1:res(1),1:res(2),1:res(3)))) scalarField_real = cmplx(0.0,0.0,pReal)
flush(6) end where
endif write(6,'(/,a,i1,1x,i1,a)') ' ... checking iFT results of compontent ', row, column, ' ..'
write(6,'(/,a,es11.4)') ' max iFT relative error = ', maxval(real(scalarField_real,pReal))
flush(6)
endif
field_real = field_real * wgt ! normalize the result by number of elements field_real = field_real * wgt ! normalize the result by number of elements
end subroutine utilities_FFTbackward end subroutine utilities_FFTbackward
@ -507,46 +533,47 @@ real(pReal) function utilities_divergenceRMS()
write(6,'(/,a)') ' ... calculating divergence ................................................' write(6,'(/,a)') ' ... calculating divergence ................................................'
flush(6) flush(6)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space ! calculating RMS divergence criterion in Fourier space
utilities_divergenceRMS = 0.0_pReal utilities_divergenceRMS = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2) do k = 1_pInt, res(3); do j = 1_pInt, res(2)
do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
utilities_divergenceRMS = utilities_divergenceRMS & utilities_divergenceRMS = utilities_divergenceRMS &
+ 2.0_pReal*(sum (real(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again + 2.0_pReal*(sum (real(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again
xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector
+sum(aimag(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& +sum(aimag(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),&
xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal))
enddo enddo
utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart utilities_divergenceRMS = utilities_divergenceRMS & ! these two layers (DC and Nyquist) do not have a conjugate complex counterpart (if res(1) /= 1)
+ sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& + sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& + sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum( real(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& + sum( real(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)& xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& + sum(aimag(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),&
xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)
enddo; enddo enddo; enddo
if(res(1) == 1_pInt) utilities_divergenceRMS = utilities_divergenceRMS * 0.5_pReal ! counted twice in case of res(1) == 1
utilities_divergenceRMS = sqrt(utilities_divergenceRMS) *wgt ! RMS in real space calculated with Parsevals theorem from Fourier space utilities_divergenceRMS = sqrt(utilities_divergenceRMS) * wgt ! RMS in real space calculated with Parsevals theorem from Fourier space
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate additional divergence criteria and report ! calculate additional divergence criteria and report
if (debugDivergence) then ! calculate divergence again if (debugDivergence) then ! calculate divergence again
err_div_max = 0.0_pReal err_div_max = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red
temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier
xi(1:3,i,j,k))*TWOPIIMG xi(1:3,i,j,k))*TWOPIIMG
err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal)) err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal))
divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared divFourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared
enddo; enddo; enddo enddo; enddo; enddo
call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real) ! already weighted call fftw_execute_dft_c2r(planDiv,divFourier,divReal) ! already weighted
err_real_div_RMS = sqrt(wgt*sum(divergence_real**2.0_pReal)) ! RMS in real space err_real_div_RMS = sqrt(wgt*sum(divReal**2.0_pReal)) ! RMS in real space
err_real_div_max = sqrt(maxval(sum(divergence_real**2.0_pReal,dim=4))) ! max in real space err_real_div_max = sqrt(maxval(sum(divReal**2.0_pReal,dim=4))) ! max in real space
err_div_max = sqrt( err_div_max) ! max in Fourier space err_div_max = sqrt( err_div_max) ! max in Fourier space
write(6,'(/,1x,a,es11.4)') 'error divergence FT RMS = ',utilities_divergenceRMS write(6,'(/,1x,a,es11.4)') 'error divergence FT RMS = ',utilities_divergenceRMS
write(6,'(1x,a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS write(6,'(1x,a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS
@ -763,10 +790,10 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,temperature,timeinc,&
integer(pInt) :: & integer(pInt) :: &
calcMode, & !< CPFEM mode for calculation calcMode, & !< CPFEM mode for calculation
collectMode !< CPFEM mode for collection collectMode, & !< CPFEM mode for collection
j,k
real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF real(pReal), dimension(3,3,3,3) :: max_dPdF, min_dPdF
real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet real(pReal) :: max_dPdF_norm, min_dPdF_norm, defgradDetMin, defgradDetMax, defgradDet
integer(pInt) :: i,j,k
write(6,'(/,a)') ' ... evaluating constitutive response ......................................' write(6,'(/,a)') ' ... evaluating constitutive response ......................................'
calcMode = CPFEM_CALCRESULTS calcMode = CPFEM_CALCRESULTS
@ -939,11 +966,11 @@ subroutine utilities_destroy()
implicit none implicit none
if (debugDivergence) call fftw_destroy_plan(plan_divergence) if (debugDivergence) call fftw_destroy_plan(planDiv)
if (debugFFTW) call fftw_destroy_plan(plan_scalarField_forth) if (debugFFTW) call fftw_destroy_plan(planDebugForth)
if (debugFFTW) call fftw_destroy_plan(plan_scalarField_back) if (debugFFTW) call fftw_destroy_plan(planDebugBack)
call fftw_destroy_plan(plan_forward) call fftw_destroy_plan(planForth)
call fftw_destroy_plan(plan_backward) call fftw_destroy_plan(planBack)
end subroutine utilities_destroy end subroutine utilities_destroy

2
code/config/numerics.config
View File

@ -66,7 +66,7 @@ itmin 2 # Minimum iteration number
maxCutBack 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc) maxCutBack 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
memory_efficient 1 # Precalculate Gamma-operator (81 double per point) memory_efficient 1 # Precalculate Gamma-operator (81 double per point)
update_gamma 0 # Update Gamma-operator with current dPdF (not possible if memory_efficient=1) update_gamma 0 # Update Gamma-operator with current dPdF (not possible if memory_efficient=1)
divergence_correction 2 # Use dimension-independent divergence criterion divergence_correction 2 # Use size-independent divergence criterion
myspectralsolver basic # Type of spectral solver (basic: basic, basicPETSc: basic with PETSc, AL: augmented Lagrange) myspectralsolver basic # Type of spectral solver (basic: basic, basicPETSc: basic with PETSc, AL: augmented Lagrange)
myfilter none # Type of filtering method to mitigate Gibb's phenomenon (none, cosine, ...) myfilter none # Type of filtering method to mitigate Gibb's phenomenon (none, cosine, ...)
petsc_options -snes_type ngmres -snes_ngmres_anderson # PetSc solver options petsc_options -snes_type ngmres -snes_ngmres_anderson # PetSc solver options

173
code/mesh.f90
View File

@ -444,9 +444,9 @@ subroutine mesh_init(ip,el)
integer(pInt), intent(in) :: el, ip integer(pInt), intent(in) :: el, ip
integer(pInt) :: j integer(pInt) :: j
write(6,'(/,a)') ' <<<+- mesh init -+>>>' write(6,'(/,a)') ' <<<+- mesh init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
write(6,'(a16,a)') ' Current time : ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem) if (allocated(mesh_mapFEtoCPelem)) deallocate(mesh_mapFEtoCPelem)
@ -486,17 +486,12 @@ subroutine mesh_init(ip,el)
do j = 1_pInt, 3_pInt do j = 1_pInt, 3_pInt
if (j /= minloc(geomdim/res,1) .and. j /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(j)*res(j) if (j /= minloc(geomdim/res,1) .and. j /= maxloc(geomdim/res,1)) scaledDim = geomdim/geomdim(j)*res(j)
enddo enddo
elseif (divergence_correction == 3_pInt) then
do j = 1_pInt, 3_pInt
if (j/=minloc(geomdim/sqrt(real(res,pReal)),1) .and. j/=maxloc(geomdim/sqrt(real(res,pReal)),1))&
scaledDim=geomdim/geomdim(j)*sqrt(real(res(j),pReal))
enddo
else else
scaledDim = geomdim scaledDim = geomdim
endif endif
write(6,'(a,3(i12 ))') ' resolution a b c:', res write(6,'(a,3(i12 ))') ' grid a b c: ', res
write(6,'(a,3(f12.5))') ' dimension x y z:', geomdim write(6,'(a,3(f12.5))') ' size x y z: ', geomdim
write(6,'(a,i5,/)') ' homogenization: ', homog write(6,'(a,i5,/)') ' homogenization: ', homog
call mesh_spectral_count_nodesAndElements call mesh_spectral_count_nodesAndElements
call mesh_spectral_count_cpElements call mesh_spectral_count_cpElements
call mesh_spectral_map_elements call mesh_spectral_map_elements
@ -566,8 +561,8 @@ end subroutine mesh_init
!! valid questions (what) are 'elem', 'node' !! valid questions (what) are 'elem', 'node'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
integer(pInt) function mesh_FEasCP(what,myID) integer(pInt) function mesh_FEasCP(what,myID)
use IO, only: &
use IO, only: IO_lc IO_lc
implicit none implicit none
character(len=*), intent(in) :: what character(len=*), intent(in) :: what
@ -654,11 +649,11 @@ end subroutine mesh_build_subNodeCoords
!> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume' !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipVolumes subroutine mesh_build_ipVolumes
use math, only: &
math_volTetrahedron, &
math_areaTriangle
use math, only: math_volTetrahedron, &
math_areaTriangle
implicit none implicit none
integer(pInt) :: e,f,t,i,j,n, & integer(pInt) :: e,f,t,i,j,n, &
NmySubNodes, & !< number of subnodes already found for this (2d) ip cell NmySubNodes, & !< number of subnodes already found for this (2d) ip cell
Ntriangles !< each interface is made up of this many triangles Ntriangles !< each interface is made up of this many triangles
@ -740,8 +735,8 @@ end subroutine mesh_build_ipVolumes
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_build_ipCoordinates subroutine mesh_build_ipCoordinates
use prec, only: &
use prec, only: tol_gravityNodePos tol_gravityNodePos
implicit none implicit none
integer(pInt) :: e,f,t,i,j,k,n integer(pInt) :: e,f,t,i,j,k,n
@ -786,11 +781,10 @@ end subroutine mesh_build_ipCoordinates
!> @brief Calculates cell center coordinates. !> @brief Calculates cell center coordinates.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function mesh_cellCenterCoordinates(i,e) pure function mesh_cellCenterCoordinates(i,e)
use prec, only: &
use prec, only: tol_gravityNodePos tol_gravityNodePos
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: e, & ! element number integer(pInt), intent(in) :: e, & ! element number
i ! integration point number i ! integration point number
@ -900,9 +894,7 @@ function mesh_spectral_getGrid(fileUnit)
if (.not. gotGrid) & if (.not. gotGrid) &
call IO_error(error_ID = 845_pInt, ext_msg='grid') call IO_error(error_ID = 845_pInt, ext_msg='grid')
if(mesh_spectral_getGrid(1) < 2_pInt .or. & ! must be at least 2 if(any(mesh_spectral_getGrid < 1_pInt)) &
mesh_spectral_getGrid(2) < 2_pInt .or. & ! -"-
mesh_spectral_getGrid(3) < 1_pInt) & ! must be at least 1
call IO_error(error_ID = 843_pInt, ext_msg='mesh_spectral_getGrid') call IO_error(error_ID = 843_pInt, ext_msg='mesh_spectral_getGrid')
end function mesh_spectral_getGrid end function mesh_spectral_getGrid
@ -1116,7 +1108,7 @@ subroutine mesh_spectral_count_cpSizes
implicit none implicit none
integer(pInt) :: t integer(pInt) :: t
t = FE_geomtype(FE_mapElemtype('C3D8R')) ! fake 3D hexahedral 8 node 1 IP element t = FE_geomtype(FE_mapElemtype('C3D8R')) ! fake 3D hexahedral 8 node 1 IP element
mesh_maxNnodes = FE_Nnodes(t) mesh_maxNnodes = FE_Nnodes(t)
mesh_maxNips = FE_Nips(t) mesh_maxNips = FE_Nips(t)
@ -1131,8 +1123,8 @@ end subroutine mesh_spectral_count_cpSizes
!! Allocates global arrays 'mesh_node0' and 'mesh_node' !! Allocates global arrays 'mesh_node0' and 'mesh_node'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine mesh_spectral_build_nodes() subroutine mesh_spectral_build_nodes()
use numerics, &
use numerics, only: numerics_unitlength only: numerics_unitlength
implicit none implicit none
integer(pInt) :: n integer(pInt) :: n
@ -1152,7 +1144,7 @@ subroutine mesh_spectral_build_nodes()
/ real(res(3),pReal) / real(res(3),pReal)
end forall end forall
mesh_node = mesh_node0 !why? mesh_node = mesh_node0 ! why?
end subroutine mesh_spectral_build_nodes end subroutine mesh_spectral_build_nodes
@ -1212,8 +1204,8 @@ subroutine mesh_spectral_build_elements(myUnit)
read(myUnit,'(a65536)') line read(myUnit,'(a65536)') line
enddo enddo
allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)) ; mesh_element = 0_pInt allocate (mesh_element (4_pInt+mesh_maxNnodes,mesh_NcpElems)); mesh_element = 0_pInt
allocate (microstructures (1_pInt+maxIntCount)) ; microstructures = 2_pInt allocate (microstructures (1_pInt+maxIntCount)); microstructures = 2_pInt
e = 0_pInt e = 0_pInt
do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!) do while (e < mesh_NcpElems .and. microstructures(1) > 0_pInt) ! fill expected number of elements, stop at end of data (or blank line!)
@ -1239,7 +1231,7 @@ subroutine mesh_spectral_build_elements(myUnit)
enddo enddo
deallocate(microstructures) deallocate(microstructures)
if (e /= mesh_NcpElems) call IO_error(880_pInt,e) if (e /= mesh_NcpElems) call IO_error(880_pInt,e) !@ToDo does that make sense?
end subroutine mesh_spectral_build_elements end subroutine mesh_spectral_build_elements
@ -1857,8 +1849,8 @@ function mesh_deformedCoordsLinear(gDim,F,FavgIn) result(coords)
! report ! report
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Restore geometry using linear integration' write(6,'(a)') ' Restore geometry using linear integration'
write(6,'(a,3(e12.5))') ' Dimension: ', gDim write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(i5))') ' Resolution:', iRes write(6,'(a,3(f12.5))') ' size x y z: ', gDim
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1961,20 +1953,20 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords)
real(pReal), intent(in), dimension(3) :: gDim real(pReal), intent(in), dimension(3) :: gDim
real(pReal), intent(in), dimension(3,3), optional :: FavgIn real(pReal), intent(in), dimension(3,3), optional :: FavgIn
real(pReal), intent(in), dimension(3), optional :: scalingIn real(pReal), intent(in), dimension(3), optional :: scalingIn
! function
! allocatable arrays for fftw c routines ! allocatable arrays for fftw c routines
type(C_PTR) :: fftw_forth, fftw_back type(C_PTR) :: planForth, planBack
type(C_PTR) :: coords_fftw, defgrad_fftw type(C_PTR) :: coords_fftw, defgrad_fftw
real(pReal), dimension(:,:,:,:,:), pointer :: F_real real(pReal), dimension(:,:,:,:,:), pointer :: F_real
complex(pReal), dimension(:,:,:,:,:), pointer :: F_fourier complex(pReal), dimension(:,:,:,:,:), pointer :: F_fourier
real(pReal), dimension(:,:,:,:), pointer :: coords_real real(pReal), dimension(:,:,:,:), pointer :: coords_real
complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier complex(pReal), dimension(:,:,:,:), pointer :: coords_fourier
! other variables ! other variables
integer(pInt) :: i, j, k, m, res1_red integer(pInt) :: i, j, k, m, res1Red
integer(pInt), dimension(3) :: k_s, iRes integer(pInt), dimension(3) :: k_s, iRes
real(pReal), dimension(3) :: scaling, step, offset_coords, integrator real(pReal), dimension(3) :: scaling, step, offset_coords, integrator
real(pReal), dimension(3,3) :: Favg real(pReal), dimension(3,3) :: Favg
integer(pInt), dimension(2:3,2) :: Nyquist ! highest frequencies to be removed (1 if even, 2 if odd)
if (present(scalingIn)) then if (present(scalingIn)) then
where (scalingIn < 0.0_pReal) ! the f2py way to tell it is not present where (scalingIn < 0.0_pReal) ! the f2py way to tell it is not present
@ -1988,53 +1980,56 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords)
iRes = [size(F,3),size(F,4),size(F,5)] iRes = [size(F,3),size(F,4),size(F,5)]
integrator = gDim / 2.0_pReal / PI ! see notes where it is used integrator = gDim / 2.0_pReal / PI ! see notes where it is used
res1_red = iRes(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) res1Red = iRes(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c)
step = gDim/real(iRes, pReal) step = gDim/real(iRes, pReal)
Nyquist(2,1:2) = [res(2)/2_pInt + 1_pInt, res(2)/2_pInt + 1_pInt + mod(res(2),2_pInt)]
Nyquist(3,1:2) = [res(3)/2_pInt + 1_pInt, res(3)/2_pInt + 1_pInt + mod(res(3),2_pInt)]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Restore geometry using FFT-based integration' write(6,'(a)') ' Restore geometry using FFT-based integration'
write(6,'(a,3(e12.5))') ' Dimension: ', gDim write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(i5))') ' Resolution:', iRes write(6,'(a,3(f12.5))') ' size x y z: ', gDim
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! sanity checks ! sanity check
if ((mod(iRes(3),2_pInt)/=0_pInt .and. iRes(3) /= 1_pInt) .or. &
mod(iRes(2),2_pInt)/=0_pInt .or. &
mod(iRes(1),2_pInt)/=0_pInt) &
call IO_error(0_pInt,ext_msg='Resolution in mesh_deformedCoordsFFT')
if (pReal /= C_DOUBLE .or. pInt /= C_INT) & if (pReal /= C_DOUBLE .or. pInt /= C_INT) &
call IO_error(0_pInt,ext_msg='Fortran to C in mesh_deformedCoordsFFT') call IO_error(0_pInt,ext_msg='Fortran to C in mesh_deformedCoordsFFT')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocation and FFTW initialization ! allocation and FFTW initialization
call fftw_set_timelimit(fftw_timelimit) defgrad_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*9_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8)
defgrad_fftw = fftw_alloc_complex(int(res1_red *iRes(2)*iRes(3)*9_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8) coords_fftw = fftw_alloc_complex(int(res1Red *iRes(2)*iRes(3)*3_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8)
call c_f_pointer(defgrad_fftw, F_real, [iRes(1)+2_pInt,iRes(2),iRes(3),3_pInt,3_pInt]) call c_f_pointer(defgrad_fftw, F_real, &
call c_f_pointer(defgrad_fftw, F_fourier,[res1_red ,iRes(2),iRes(3),3_pInt,3_pInt]) [iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt,3_pInt])
coords_fftw = fftw_alloc_complex(int(res1_red *iRes(2)*iRes(3)*3_pInt,C_SIZE_T)) ! C_SIZE_T is of type integer(8) call c_f_pointer(defgrad_fftw, F_fourier, &
call c_f_pointer(coords_fftw, coords_real, [iRes(1)+2_pInt,iRes(2),iRes(3),3_pInt]) [res1Red, iRes(2),iRes(3),3_pInt,3_pInt])
call c_f_pointer(coords_fftw, coords_fourier, [res1_red ,iRes(2),iRes(3),3_pInt]) call c_f_pointer(coords_fftw, coords_real, &
fftw_forth = fftw_plan_many_dft_r2c(3_pInt,[iRes(3),iRes(2) ,iRes(1)],9_pInt,& ! dimensions , length in each dimension in reversed order [iRes(1)+2_pInt-mod(iRes(1),2_pInt),iRes(2),iRes(3),3_pInt])
F_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt],& ! input data , physical length in each dimension in reversed order call c_f_pointer(coords_fftw, coords_fourier, &
1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions [res1Red, iRes(2),iRes(3),3_pInt])
F_fourier,[iRes(3),iRes(2) ,res1_red],&
1_pInt, iRes(3)*iRes(2)* res1_red,fftw_planner_flag)
fftw_back = fftw_plan_many_dft_c2r(3_pInt,[iRes(3),iRes(2) ,iRes(1)],3_pInt,& call fftw_set_timelimit(fftw_timelimit)
coords_fourier,[iRes(3),iRes(2) ,res1_red],& planForth = fftw_plan_many_dft_r2c(3_pInt,[iRes(3),iRes(2) ,iRes(1)],9_pInt,& ! dimensions , length in each dimension in reversed order
1_pInt, iRes(3)*iRes(2)* res1_red,& F_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],& ! input data , physical length in each dimension in reversed order
coords_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt],& 1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),& ! striding , product of physical lenght in the 3 dimensions
1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt),fftw_planner_flag) F_fourier,[iRes(3),iRes(2) ,res1Red],&
forall(k = 1_pInt:iRes(3), j = 1_pInt:iRes(2), i = 1_pInt:iRes(1)) & 1_pInt, iRes(3)*iRes(2)* res1Red,fftw_planner_flag)
F_real(i,j,k,1:3,1:3) = F(1:3,1:3,i,j,k) ! F_real is overwritten during plan creation and is larger (padding)
planBack = fftw_plan_many_dft_c2r(3_pInt,[iRes(3),iRes(2) ,iRes(1)],3_pInt,&
coords_fourier,[iRes(3),iRes(2) ,res1Red],&
1_pInt, iRes(3)*iRes(2)* res1Red,&
coords_real,[iRes(3),iRes(2) ,iRes(1)+2_pInt-mod(iRes(1),2_pInt)],&
1_pInt, iRes(3)*iRes(2)*(iRes(1)+2_pInt-mod(iRes(1),2_pInt)),&
fftw_planner_flag)
F_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3,1:3) = &
reshape(F,[res(1),res(2),res(3),3,3], order = [4,5,1,2,3]) ! F_real is overwritten during plan creatio, is larger (padding) and has different order
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! FFT ! FFT
call fftw_execute_dft_r2c(fftw_forth, F_real, F_fourier) call fftw_execute_dft_r2c(planForth, F_real, F_fourier)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! if no average F is given, compute it in Fourier space ! if no average F is given, compute it in Fourier space
@ -2050,13 +2045,16 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! remove highest frequency in each direction, in third direction only if not 2D ! remove highest frequency in each direction, in third direction only if not 2D
F_fourier( iRes(1)/2_pInt+1_pInt,1:iRes(2) ,1:iRes(3) ,1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal) if(iRes(1)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
F_fourier( 1:res1_red ,iRes(2)/2_pInt+1_pInt,1:iRes(3) ,1:3,1:3) & F_fourier (res1Red, 1:iRes(2), 1:iRes(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal) = cmplx(0.0_pReal,0.0_pReal,pReal)
if(iRes(3)>1_pInt) & if(iRes(2)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
F_fourier(1:res1_red ,1:iRes(2) ,iRes(3)/2_pInt+1_pInt,1:3,1:3) & F_fourier (1:res1Red,Nyquist(2,1):Nyquist(2,2),1:iRes(3), 1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal) = cmplx(0.0_pReal,0.0_pReal,pReal)
if(iRes(3)/=1_pInt) & ! do not delete the whole slice in case of 2D calculation
F_fourier (1:res1Red,1:iRes(2), Nyquist(3,1):Nyquist(3,2),1:3,1:3) &
= cmplx(0.0_pReal,0.0_pReal,pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! integration in Fourier space ! integration in Fourier space
@ -2067,7 +2065,7 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords)
do j = 1_pInt, iRes(2) do j = 1_pInt, iRes(2)
k_s(2) = j-1_pInt k_s(2) = j-1_pInt
if(j > iRes(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-iRes(2) if(j > iRes(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-iRes(2)
do i = 1_pInt, res1_red do i = 1_pInt, res1Red
k_s(1) = i-1_pInt k_s(1) = i-1_pInt
do m = 1_pInt,3_pInt do m = 1_pInt,3_pInt
coords_fourier(i,j,k,m) = sum(F_fourier(i,j,k,m,1:3)*& coords_fourier(i,j,k,m) = sum(F_fourier(i,j,k,m,1:3)*&
@ -2079,12 +2077,11 @@ function mesh_deformedCoordsFFT(gDim,F,FavgIn,scalingIn) result(coords)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! iFFT and freeing memory ! iFFT and freeing memory
call fftw_execute_dft_c2r(fftw_back,coords_fourier,coords_real) call fftw_execute_dft_c2r(planBack,coords_fourier,coords_real)
coords_real = coords_real/real(product(iRes),pReal) coords = reshape(coords_real(1:iRes(1),1:iRes(2),1:iRes(3),1:3), [3,iRes(1),iRes(2),iRes(3)], &
forall(k = 1_pInt:iRes(3), j = 1_pInt:iRes(2), i = 1_pInt:iRes(1)) & order = [2,3,4,1])/real(product(iRes),pReal) ! weight and change order
coords(1:3,i,j,k) = coords_real(i,j,k,1:3) call fftw_destroy_plan(planForth)
call fftw_destroy_plan(fftw_forth) call fftw_destroy_plan(planBack)
call fftw_destroy_plan(fftw_back)
call fftw_free(defgrad_fftw) call fftw_free(defgrad_fftw)
call fftw_free(coords_fftw) call fftw_free(coords_fftw)
@ -2137,8 +2134,8 @@ function mesh_volumeMismatch(gDim,F,nodes) result(vMismatch)
! report and check ! report and check
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Calculating volume mismatch' write(6,'(a)') ' Calculating volume mismatch'
write(6,'(a,3(e12.5))') ' Dimension: ', gDim write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(i5))') ' Resolution:', iRes write(6,'(a,3(f12.5))') ' size x y z: ', gDim
endif endif
if (any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]))& if (any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]))&
@ -2209,8 +2206,8 @@ function mesh_shapeMismatch(gDim,F,nodes,centres) result(sMismatch)
! report and check ! report and check
if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_mesh),debug_levelBasic) /= 0_pInt) then
write(6,'(a)') ' Calculating shape mismatch' write(6,'(a)') ' Calculating shape mismatch'
write(6,'(a,3(e12.5))') ' Dimension: ', gDim write(6,'(a,3(i12 ))') ' grid a b c: ', iRes
write(6,'(a,3(i5))') ' Resolution:', iRes write(6,'(a,3(f12.5))') ' size x y z: ', gDim
endif endif
if(any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]) .or.& if(any([iRes/=size(nodes,2)-1_pInt,iRes/=size(nodes,3)-1_pInt,iRes/=size(nodes,4)-1_pInt]) .or.&
@ -3553,8 +3550,8 @@ subroutine mesh_build_ipAreas
enddo enddo
enddo enddo
end select end select
enddo enddo
end subroutine mesh_build_ipAreas end subroutine mesh_build_ipAreas
@ -4610,6 +4607,12 @@ subroutine mesh_build_FEdata
! indicates which subnodes make up the interfaces enclosing the IP volume. ! indicates which subnodes make up the interfaces enclosing the IP volume.
! The sorting convention is such that the outward pointing normal ! The sorting convention is such that the outward pointing normal
! follows from a right-handed traversal of the face node list. ! follows from a right-handed traversal of the face node list.
! For two-dimensional elements, which only have lines as "interface"
! one nevertheless has to specify each interface by a closed path,
! e.g., 1,2, 2,1, assuming the line connects nodes 1 and 2.
! This will result in zero ipVolume and interfaceArea, but is not
! detrimental at the moment since non-local constitutive laws are
! currently not foreseen in 2D cases.
me = 0_pInt me = 0_pInt
me = me + 1_pInt me = me + 1_pInt

6
code/numerics.f90
View File

@ -107,8 +107,8 @@ module numerics
itmin = 2_pInt, & !< minimum number of iterations itmin = 2_pInt, & !< minimum number of iterations
maxCutBack = 3_pInt, & !< max number of cut backs maxCutBack = 3_pInt, & !< max number of cut backs
regridMode = 0_pInt, & !< 0: no regrid; 1: regrid if DAMASK doesn't converge; 2: regrid if DAMASK or BVP Solver doesn't converge regridMode = 0_pInt, & !< 0: no regrid; 1: regrid if DAMASK doesn't converge; 2: regrid if DAMASK or BVP Solver doesn't converge
divergence_correction = 2_pInt !< correct divergence calculation in fourier space 0: no correction, 1: dimension scaled to 1, 2: dimension scaled to Npoints, 3: dimension scaled to sqrt(Npoints) divergence_correction = 2_pInt !< correct divergence calculation in fourier space 0: no correction, 1: size scaled to 1, 2: size scaled to Npoints
logical, protected , public :: & logical, protected, public :: &
memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate
update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness
#endif #endif
@ -484,7 +484,7 @@ subroutine numerics_init
if (itmax <= 1_pInt) call IO_error(301_pInt,ext_msg='itmax') if (itmax <= 1_pInt) call IO_error(301_pInt,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1_pInt) call IO_error(301_pInt,ext_msg='itmin') if (itmin > itmax .or. itmin < 1_pInt) call IO_error(301_pInt,ext_msg='itmin')
if (divergence_correction < 0_pInt .or. & if (divergence_correction < 0_pInt .or. &
divergence_correction > 3_pInt) call IO_error(301_pInt,ext_msg='divergence_correction') divergence_correction > 2_pInt) call IO_error(301_pInt,ext_msg='divergence_correction')
if (maxCutBack < 0_pInt) call IO_error(301_pInt,ext_msg='maxCutBack') if (maxCutBack < 0_pInt) call IO_error(301_pInt,ext_msg='maxCutBack')
if (update_gamma .and. & if (update_gamma .and. &
.not. memory_efficient) call IO_error(error_ID = 847_pInt) .not. memory_efficient) call IO_error(error_ID = 847_pInt)