From 18b89239291d62ff2239107108854d0570772632 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 09:25:06 +0100 Subject: [PATCH] centralize FFTs --- src/grid/grid_damage_spectral.f90 | 19 ++++++---------- src/grid/grid_mech_spectral_basic.f90 | 4 ---- src/grid/grid_mech_spectral_polarisation.f90 | 4 ---- src/grid/grid_thermal_spectral.f90 | 23 +++++++++----------- src/grid/spectral_utilities.f90 | 14 +++++++----- 5 files changed, 25 insertions(+), 39 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index fc1defac3..f09992f55 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -253,11 +253,13 @@ end function grid_damage_spectral_solution subroutine grid_damage_spectral_forward(cutBack) logical, intent(in) :: cutBack + integer :: i, j, k, ce DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: phi_PETSc PetscErrorCode :: err_PETSc + if (cutBack) then phi_current = phi_lastInc phi_stagInc = phi_lastInc @@ -284,7 +286,7 @@ end subroutine grid_damage_spectral_forward !-------------------------------------------------------------------------------------------------- -!> @brief forms the spectral damage residual vector +!> @brief Construct the residual vector. !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,r,dummy,err_PETSc) @@ -297,7 +299,8 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & r PetscObject :: dummy - PetscErrorCode :: err_PETSc + PetscErrorCode, intent(out) :: err_PETSc + integer :: i, j, k, ce @@ -305,17 +308,13 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) !-------------------------------------------------------------------------------------------------- ! evaluate polarization field scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_current - call utilities_FFTscalarForward - call utilities_fourierScalarGradient !< calculate gradient of damage field - call utilities_FFTvectorBackward + call utilities_fourierScalarGradient() ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 vectorField_real(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField_real(1:3,i,j,k)) end do; end do; end do - call utilities_FFTvectorForward - call utilities_fourierVectorDivergence !< calculate damage divergence in fourier field - call utilities_FFTscalarBackward + call utilities_fourierVectorDivergence() ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 @@ -324,11 +323,7 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) + mu_ref*phi_current(i,j,k) end do; end do; end do -!-------------------------------------------------------------------------------------------------- -! convolution of damage field with green operator - call utilities_FFTscalarForward call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) - call utilities_FFTscalarBackward !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 4dca95e80..a065633de 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -544,12 +544,8 @@ subroutine formResidual(in, F, & F_aim = F_aim - deltaF_aim err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask))) -!-------------------------------------------------------------------------------------------------- -! updated deformation gradient using fix point algorithm of basic scheme tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r ! store fPK field for subsequent FFT forward transform - call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier - call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 9c9a83766..fea34d7d0 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -618,11 +618,7 @@ subroutine formResidual(in, FandF_tau, & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) end do; end do; end do -!-------------------------------------------------------------------------------------------------- -! doing convolution in Fourier space - call utilities_FFTtensorForward call utilities_fourierGammaConvolution(params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) - call utilities_FFTtensorBackward !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index b34717027..c6b190254 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -242,11 +242,13 @@ end function grid_thermal_spectral_solution subroutine grid_thermal_spectral_forward(cutBack) logical, intent(in) :: cutBack + integer :: i, j, k, ce DM :: dm_local PetscScalar, dimension(:,:,:), pointer :: T_PETSc PetscErrorCode :: err_PETSc + if (cutBack) then T_current = T_lastInc T_stagInc = T_lastInc @@ -307,7 +309,7 @@ end subroutine grid_thermal_spectral_restartWrite !-------------------------------------------------------------------------------------------------- -!> @brief forms the spectral thermal residual vector +!> @brief Construct the residual vector. !-------------------------------------------------------------------------------------------------- subroutine formResidual(in,x_scal,r,dummy,err_PETSc) @@ -320,24 +322,22 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) X_RANGE,Y_RANGE,Z_RANGE), intent(out) :: & r PetscObject :: dummy - PetscErrorCode :: err_PETSc + PetscErrorCode, intent(out) :: err_PETSc + integer :: i, j, k, ce + T_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field scalarField_real(1:cells(1),1:cells(2),1:cells3) = T_current - call utilities_FFTscalarForward - call utilities_fourierScalarGradient !< calculate gradient of temperature field - call utilities_FFTvectorBackward + call utilities_fourierScalarGradient() ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 vectorField_real(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField_real(1:3,i,j,k)) end do; end do; end do - call utilities_FFTvectorForward - call utilities_fourierVectorDivergence !< calculate temperature divergence in fourier field - call utilities_FFTscalarBackward + call utilities_fourierVectorDivergence() ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 @@ -346,15 +346,12 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) + mu_ref*T_current(i,j,k) end do; end do; end do -!-------------------------------------------------------------------------------------------------- -! convolution of temperature field with green operator - call utilities_FFTscalarForward call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) - call utilities_FFTscalarBackward !-------------------------------------------------------------------------------------------------- ! constructing residual - r = T_current - scalarField_real(1:cells(1),1:cells(2),1:cells3) + r = T_current & + - scalarField_real(1:cells(1),1:cells(2),1:cells3) err_PETSc = 0 end subroutine formResidual diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index db2efbfec..858fc20db 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -116,12 +116,6 @@ module spectral_utilities public :: & spectral_utilities_init, & utilities_updateGamma, & - utilities_FFTtensorForward, & - utilities_FFTtensorBackward, & - utilities_FFTvectorForward, & - utilities_FFTvectorBackward, & - utilities_FFTscalarForward, & - utilities_FFTscalarBackward, & utilities_fourierGammaConvolution, & utilities_fourierGreenConvolution, & utilities_divergenceRMS, & @@ -526,6 +520,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim) print'(/,1x,a)', '... doing gamma convolution ...............................................' flush(IO_STDOUT) + call utilities_FFTtensorForward() !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation (mechanical equilibrium) memoryEfficient: if (num%memory_efficient) then @@ -587,6 +582,7 @@ subroutine utilities_fourierGammaConvolution(fieldAim) end if memoryEfficient if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) + call utilities_FFTtensorBackward() end subroutine utilities_fourierGammaConvolution @@ -603,6 +599,7 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t) !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation + call utilities_FFTscalarForward() !$OMP PARALLEL DO PRIVATE(GreenOp_hat) do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red GreenOp_hat = cmplx(1.0_pReal,0.0_pReal,pReal) & @@ -611,6 +608,7 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t) scalarField_fourier(i,k,j) = scalarField_fourier(i,k,j)*GreenOp_hat end do; end do; end do !$OMP END PARALLEL DO + call utilities_FFTscalarBackward() end subroutine utilities_fourierGreenConvolution @@ -810,9 +808,11 @@ subroutine utilities_fourierScalarGradient() integer :: i, j, k + call utilities_FFTscalarForward() do j = 1, cells2; do k = 1, cells(3); do i = 1,cells1Red vectorField_fourier(1:3,i,k,j) = scalarField_fourier(i,k,j)*xi1st(1:3,i,k,j) ! ToDo: no -conjg? end do; end do; end do + call utilities_FFTvectorBackward() end subroutine utilities_fourierScalarGradient @@ -822,8 +822,10 @@ end subroutine utilities_fourierScalarGradient !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierVectorDivergence() + call utilities_FFTvectorForward() scalarField_fourier(1:cells1Red,1:cells(3),1:cells2) = sum(vectorField_fourier(1:3,1:cells1Red,1:cells(3),1:cells2) & *conjg(-xi1st),1) + call utilities_FFTscalarBackward() end subroutine utilities_fourierVectorDivergence