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

improved update of ip coordinates for spectral solver. 

do not need to create a new fftw plan all the time, using data already defined for the convolution
This commit is contained in:
Martin Diehl 2015-03-12 22:28:33 +00:00
parent 263f997cf2
commit 00cba25a44
2 changed files with 98 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -111,13 +111,11 @@ subroutine basicPETSc_init(temperature)
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
utilities_updateIPcoords, &
grid, &
grid1Red, &
wgt, &
geomSize
use mesh, only: &
mesh_ipCoordinates, &
mesh_deformedCoordsFFT
geomSize
use math, only: &
math_invSym3333
@ -194,8 +192,7 @@ subroutine basicPETSc_init(temperature)
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(F_lastInc, &
reshape(F(0:8,0:grid(1)-1_pInt,0:grid(2)-1_pInt,0:grid(3)-1_pInt),[3,3,grid(1),grid(2),grid(3)]), &
temperature, &
@ -474,11 +471,9 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
geomSize, &
Utilities_calculateRate, &
Utilities_forwardField, &
utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use mesh, only: &
mesh_ipCoordinates,&
mesh_deformedCoordsFFT
use IO, only: &
IO_write_JobRealFile
use FEsolving, only: &
@ -523,8 +518,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
write (777,rec=1) C_volAvgLastInc
close(777)
endif
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call utilities_updateIPcoords(F)
if (cutBack) then
F_aim = F_aim_lastInc
F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)])
@ -546,8 +541,7 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
mesh_ipCoordinates = reshape(mesh_deformedCoordsFFT(geomSize,reshape(&
F,[3,3,grid(1),grid(2),grid(3)])),[3,1,product(grid)])
call utilities_updateIPcoords(F)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid(3)]))
F_lastInc2 = F_lastInc

96
code/DAMASK_spectral_utilities.f90
View File

@ -37,8 +37,13 @@ module DAMASK_spectral_utilities
!--------------------------------------------------------------------------------------------------
! debug fftw
complex(pReal),private, dimension(:,:,:), pointer :: scalarField_real, & !< scalar field real representation for debug of FFTW
scalarField_fourier !< scalar field complex representation for debug of FFTW
complex(pReal),private, dimension(:,:,:), pointer :: scalarField_real, & !< scalar field real representation for debug of FFTW
scalarField_fourier !< scalar field complex representation for debug of FFTW
!--------------------------------------------------------------------------------------------------
! geometry reconstruction
!real(pReal), private, dimension(:,:,:,:), pointer :: coords_real
!complex(pReal), private, dimension(:,:,:,:), pointer :: coords_fourier
!--------------------------------------------------------------------------------------------------
! debug divergence
@ -52,7 +57,8 @@ module DAMASK_spectral_utilities
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
planDiv!, & !< FFTW plan in case of debugging divergence calculation
!planCoords !< FFTW plan for geometry reconstruction
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
@ -100,7 +106,8 @@ module DAMASK_spectral_utilities
utilities_constitutiveResponse, &
utilities_calculateRate, &
utilities_forwardField, &
utilities_destroy
utilities_destroy, &
utilities_updateIPcoords
private :: &
utilities_getFilter
@ -162,7 +169,8 @@ subroutine utilities_init()
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_fourierC, & !< field cotaining data for FFTW in fourier space when debugging FFTW (no in place)
div !< 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)
!coords_fftw
write(6,'(/,a)') ' <<<+- DAMASK_spectral_utilities init -+>>>'
write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -223,6 +231,9 @@ subroutine utilities_init()
tensorField = fftw_alloc_complex(int(grid1Red*grid(2)*grid(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, [grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3),3,3])! place a pointer for a real representation on tensorField
call c_f_pointer(tensorField, field_fourier,[grid1Red, grid(2),grid(3),3,3])! place a pointer for a complex representation on tensorField
!coords_fftw = fftw_alloc_complex(int(grid1Red*grid(2)*grid(3)*3_pInt,C_SIZE_T)) ! allocate aligned data using a C function, C_SIZE_T is of type integer(8)
!call c_f_pointer(tensorField, coords_real,[grid(1)+2_pInt-mod(grid(1),2_pInt),grid(2),grid(3),3]) ! place a pointer for a real representation on coords_fftw
!call c_f_pointer(tensorField, coords_fourier,[grid1Red, grid(2),grid(3),3]) ! place a pointer for a real representation on coords_fftw
!--------------------------------------------------------------------------------------------------
! general initialization of FFTW (see manual on fftw.org for more details)
@ -244,6 +255,17 @@ subroutine utilities_init()
1, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)), & ! striding, product of physical length in the 3 dimensions
fftw_planner_flag) ! planner precision
!--------------------------------------------------------------------------------------------------
! allocation and FFTW initialization
! planCoords = fftw_plan_many_dft_c2r(3_pInt,[grid(3),grid(2) ,grid(1)],3_pInt,&
! coords_fourier,[grid(3),grid(2) ,grid1Red],&
! 1_pInt, grid(3)*grid(2)* grid1Red,&
! coords_real,[grid(3),grid(2) ,grid(1)+2_pInt-mod(grid(1),2_pInt)],&
! 1_pInt, grid(3)*grid(2)*(grid(1)+2_pInt-mod(grid(1),2_pInt)),&
! fftw_planner_flag)
!--------------------------------------------------------------------------------------------------
! depending on debug options, allocate more memory and create additional plans
if (debugDivergence) then
@ -992,8 +1014,72 @@ subroutine utilities_destroy()
if (debugFFTW) call fftw_destroy_plan(planDebugBack)
call fftw_destroy_plan(planForth)
call fftw_destroy_plan(planBack)
!call fftw_destroy_plan(planCoords)
end subroutine utilities_destroy
!--------------------------------------------------------------------------------------------------
!> @brief calculate coordinates in current configuration for given defgrad field
! using integration in Fourier space. Similar as in mesh.f90, but using data already defined for
! convolution
!--------------------------------------------------------------------------------------------------
subroutine utilities_updateIPcoords(F)
use math
use mesh, only: &
mesh_ipCoordinates
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid(3)), intent(in) :: F
integer(pInt) :: i, j, k, m
integer(pInt), dimension(3) :: k_s
real(pReal), dimension(3) :: step, offset_coords, integrator
real(pReal), dimension(3,3) :: Favg
field_real = 0.0_pReal
field_real(1:grid(1),1:grid(2),1:grid(3),1:3,1:3) = reshape(F,[grid(1),grid(2),grid(3),3,3],&
order=[4,5,1,2,3])
call utilities_FFTforward()
integrator = geomsize * 0.5_pReal / PI
step = geomsize/real(grid, pReal)
!--------------------------------------------------------------------------------------------------
! average F
Favg = real(field_fourier(1,1,1,1:3,1:3),pReal)/real(product(grid),pReal)
!--------------------------------------------------------------------------------------------------
! integration in Fourier space
field_fourier = cmplx(0.0_pReal,0.0_pReal,pReal)
do k = 1_pInt, grid(3)
k_s(3) = k-1_pInt
if(k > grid(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-grid(3)
do j = 1_pInt, grid(2)
k_s(2) = j-1_pInt
if(j > grid(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-grid(2)
do i = 1_pInt, grid1Red
k_s(1) = i-1_pInt
do m = 1_pInt,3_pInt
field_fourier(i,j,k,m,1) = sum(field_fourier(i,j,k,m,1:3)*&
cmplx(0.0_pReal,real(k_s,pReal)*integrator,pReal))
enddo
if (any(k_s /= 0_pInt)) field_fourier(i,j,k,1:3,1) = &
field_fourier(i,j,k,1:3,1) / cmplx(-sum(k_s*k_s),0.0_pReal,pReal)
enddo; enddo; enddo
call utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! add average to fluctuation and put (0,0,0) on (0,0,0)
offset_coords = math_mul33x3(Favg,step/2.0_pReal) - field_real(1,1,1,1:3,1)
m = 1_pInt
do k = 1_pInt,grid(3); do j = 1_pInt,grid(2); do i = 1_pInt,grid(1)
mesh_ipCoordinates(1:3,1,m) = field_real(i,j,k,1:3,1) &
+ offset_coords &
+ math_mul33x3(Favg,step*real([i,j,k]-1_pInt,pReal))
m = m+1_pInt
enddo; enddo; enddo
end subroutine utilities_updateIPcoords
end module DAMASK_spectral_utilities