From ce98cfdd5eb6d384405480d8988d5a5f757fdf47 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 07:55:22 +0100 Subject: [PATCH 01/14] padding is handled centrally in the FFT forward routines --- src/grid/grid_damage_spectral.f90 | 1 - src/grid/grid_mech_spectral_basic.f90 | 1 - src/grid/grid_mech_spectral_polarisation.f90 | 3 --- src/grid/grid_thermal_spectral.f90 | 1 - 4 files changed, 6 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 816cdd4b1..51804d52b 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -304,7 +304,6 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) phi_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_real = 0.0_pReal scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of damage field diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index bead280a7..ef958d718 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -545,7 +545,6 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme - tensorField_real = 0.0_pReal 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" err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 2b4ea364a..c82048d77 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -611,7 +611,6 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! - tensorField_real = 0.0_pReal do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) tensorField_real(1:3,1:3,i,j,k) = & num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& @@ -643,7 +642,6 @@ subroutine formResidual(in, FandF_tau, & P_av-P_aim, & params%stress_mask))) ! calculate divergence - tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r_F !< stress field in disguise call utilities_FFTtensorForward err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress @@ -662,7 +660,6 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! calculating curl - tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F call utilities_FFTtensorForward err_curl = utilities_curlRMS() diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 29ae07769..b34717027 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -326,7 +326,6 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) T_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_real = 0.0_pReal scalarField_real(1:cells(1),1:cells(2),1:cells3) = T_current call utilities_FFTscalarForward call utilities_fourierScalarGradient !< calculate gradient of temperature field From 20da5663c0b233d839c31b34f30a325faeadb955 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 08:14:38 +0100 Subject: [PATCH 02/14] simplified, avoid intermediate writes --- src/grid/grid_damage_spectral.f90 | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 51804d52b..fc1defac3 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -330,14 +330,10 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) call utilities_FFTscalarBackward - where(scalarField_real(1:cells(1),1:cells(2),1:cells3) > phi_lastInc) & - scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_lastInc - where(scalarField_real(1:cells(1),1:cells(2),1:cells3) < num%residualStiffness) & - scalarField_real(1:cells(1),1:cells(2),1:cells3) = num%residualStiffness - !-------------------------------------------------------------------------------------------------- ! constructing residual - r = scalarField_real(1:cells(1),1:cells(2),1:cells3) - phi_current + r = max(min(scalarField_real(1:cells(1),1:cells(2),1:cells3),phi_lastInc),num%residualStiffness) & + - phi_current err_PETSc = 0 end subroutine formResidual From cd2a21509a3294548d069eb64647a8163a801517 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 08:55:15 +0100 Subject: [PATCH 03/14] avoid depenencies on global state requires on extra forward FFT pre iteration for basic scheme --- src/grid/grid_mech_spectral_basic.f90 | 2 +- src/grid/grid_mech_spectral_polarisation.f90 | 12 ++---------- src/grid/spectral_utilities.f90 | 18 ++++++++++++++---- 3 files changed, 17 insertions(+), 15 deletions(-) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index ef958d718..4dca95e80 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -536,6 +536,7 @@ subroutine formResidual(in, F, & F,params%Delta_t,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + err_div = utilities_divergenceRMS(r) !-------------------------------------------------------------------------------------------------- ! stress BC handling @@ -547,7 +548,6 @@ subroutine formResidual(in, F, & ! 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" - err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use 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 diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index c82048d77..9c9a83766 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -634,6 +634,8 @@ subroutine formResidual(in, FandF_tau, & P_av,C_volAvg,C_minMaxAvg, & F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + err_div = utilities_divergenceRMS(r_F) + err_curl = utilities_curlRMS(F) !-------------------------------------------------------------------------------------------------- ! stress BC handling @@ -641,10 +643,6 @@ subroutine formResidual(in, FandF_tau, & err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), & P_av-P_aim, & params%stress_mask))) -! calculate divergence - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = r_F !< stress field in disguise - call utilities_FFTtensorForward - err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -658,12 +656,6 @@ subroutine formResidual(in, FandF_tau, & + r_F_tau(1:3,1:3,i,j,k) end do; end do; end do -!-------------------------------------------------------------------------------------------------- -! calculating curl - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F - call utilities_FFTtensorForward - err_curl = utilities_curlRMS() - end subroutine formResidual end module grid_mechanical_spectral_polarisation diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 6cb7edd30..db2efbfec 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -616,9 +616,11 @@ end subroutine utilities_fourierGreenConvolution !-------------------------------------------------------------------------------------------------- -!> @brief calculate root mean square of divergence of field_fourier +!> @brief Calculate root mean square of divergence. !-------------------------------------------------------------------------------------------------- -real(pReal) function utilities_divergenceRMS() +real(pReal) function utilities_divergenceRMS(tensorField) + + real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: tensorField integer :: i, j, k integer(MPI_INTEGER_KIND) :: err_MPI @@ -628,6 +630,9 @@ real(pReal) function utilities_divergenceRMS() print'(/,1x,a)', '... calculating divergence ................................................' flush(IO_STDOUT) + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = tensorField + call utilities_FFTtensorforward() + rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal) !-------------------------------------------------------------------------------------------------- @@ -660,9 +665,11 @@ end function utilities_divergenceRMS !-------------------------------------------------------------------------------------------------- -!> @brief calculate max of curl of field_fourier +!> @brief Calculate root mean square of curl. !-------------------------------------------------------------------------------------------------- -real(pReal) function utilities_curlRMS() +real(pReal) function utilities_curlRMS(tensorField) + + real(pReal), dimension(3,3,cells(1),cells(2),cells3), intent(in) :: tensorField integer :: i, j, k, l integer(MPI_INTEGER_KIND) :: err_MPI @@ -673,6 +680,9 @@ real(pReal) function utilities_curlRMS() print'(/,1x,a)', '... calculating curl ......................................................' flush(IO_STDOUT) + tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = tensorField + call utilities_FFTtensorforward() + rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal) !-------------------------------------------------------------------------------------------------- From 18b89239291d62ff2239107108854d0570772632 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 09:25:06 +0100 Subject: [PATCH 04/14] 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 From cb6df618feca7068d2933c116e9a53b36d736cfd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 10:50:15 +0100 Subject: [PATCH 05/14] avoid global variables --- src/grid/grid_mech_spectral_basic.f90 | 15 ++++++------- src/grid/grid_mech_spectral_polarisation.f90 | 15 +++++++------ src/grid/spectral_utilities.f90 | 22 +++++++++++++------- 3 files changed, 27 insertions(+), 25 deletions(-) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index a065633de..1ec6d334c 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -491,7 +491,7 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the residual vector +!> @brief Construct the residual vector. !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, F, & r, dummy, err_PETSc) @@ -501,15 +501,17 @@ subroutine formResidual(in, F, & intent(in) :: F !< deformation gradient field PetscScalar, dimension(3,3,X_RANGE,Y_RANGE,Z_RANGE), & intent(out) :: r !< residuum field + PetscObject :: dummy + PetscErrorCode :: err_PETSc + real(pReal), dimension(3,3) :: & deltaF_aim PetscInt :: & PETScIter, & nfuncs - PetscObject :: dummy - PetscErrorCode :: err_PETSc integer(MPI_INTEGER_KIND) :: err_MPI + call SNESGetNumberFunctionEvals(SNES_mechanical,nfuncs,err_PETSc) CHKERRQ(err_PETSc) call SNESGetIterationNumber(SNES_mechanical,PETScIter,err_PETSc) @@ -544,12 +546,7 @@ subroutine formResidual(in, F, & F_aim = F_aim - deltaF_aim err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask))) - 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_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier - -!-------------------------------------------------------------------------------------------------- -! constructing residual - r = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too + r = utilities_fourierGammaConvolution(r,params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier end subroutine formResidual diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index fea34d7d0..235bb873e 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -551,7 +551,7 @@ end subroutine converged !-------------------------------------------------------------------------------------------------- -!> @brief forms the residual vector +!> @brief Construct the residual vector. !-------------------------------------------------------------------------------------------------- subroutine formResidual(in, FandF_tau, & r, dummy,err_PETSc) @@ -561,6 +561,9 @@ subroutine formResidual(in, FandF_tau, & target, intent(in) :: FandF_tau PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& target, intent(out) :: r !< residuum field + PetscObject :: dummy + PetscErrorCode :: err_PETSc + PetscScalar, pointer, dimension(:,:,:,:,:) :: & F, & F_tau, & @@ -569,13 +572,10 @@ subroutine formResidual(in, FandF_tau, & PetscInt :: & PETScIter, & nfuncs - PetscObject :: dummy - PetscErrorCode :: err_PETSc integer(MPI_INTEGER_KIND) :: err_MPI integer :: & i, j, k, e -!--------------------------------------------------------------------------------------------------- F => FandF_tau(1:3,1:3,1,& XG_RANGE,YG_RANGE,ZG_RANGE) @@ -612,17 +612,16 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) - tensorField_real(1:3,1:3,i,j,k) = & + r_F_tau(1:3,1:3,i,j,k) = & num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& num%alpha*matmul(F(1:3,1:3,i,j,k), & 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 - call utilities_fourierGammaConvolution(params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) - !-------------------------------------------------------------------------------------------------- ! constructing residual - r_F_tau = num%beta*F - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) + r_F_tau = num%beta*F & + - utilities_fourierGammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 858fc20db..9e5282de0 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -42,7 +42,7 @@ module spectral_utilities !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW - real(C_DOUBLE), public, dimension(:,:,:,:,:), pointer :: tensorField_real !< tensor field in real space + real(C_DOUBLE), dimension(:,:,:,:,:), pointer :: tensorField_real !< tensor field in real space real(C_DOUBLE), public, dimension(:,:,:,:), pointer :: vectorField_real !< vector field in real space real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field in real space complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:,:), pointer :: tensorField_fourier !< tensor field in Fourier space @@ -505,12 +505,14 @@ end subroutine utilities_FFTvectorBackward !-------------------------------------------------------------------------------------------------- !> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim !-------------------------------------------------------------------------------------------------- -subroutine utilities_fourierGammaConvolution(fieldAim) +function utilities_fourierGammaConvolution(field, fieldAim) result(gammaField) + real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: field real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution - complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx - real(pReal), dimension(6,6) :: A, A_inv + real(pReal), dimension(3,3,cells(1),cells(2),cells3) :: gammaField + complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx + real(pReal), dimension(6,6) :: A, A_inv integer :: & i, j, k, & l, m, n, o @@ -520,8 +522,11 @@ subroutine utilities_fourierGammaConvolution(fieldAim) print'(/,1x,a)', '... doing gamma convolution ...............................................' flush(IO_STDOUT) - call utilities_FFTtensorForward() -!-------------------------------------------------------------------------------------------------- + tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = field + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) + + !-------------------------------------------------------------------------------------------------- ! do the actual spectral method calculation (mechanical equilibrium) memoryEfficient: if (num%memory_efficient) then !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat) @@ -582,9 +587,10 @@ 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() + call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) + gammaField = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3)*wgt -end subroutine utilities_fourierGammaConvolution +end function utilities_fourierGammaConvolution !-------------------------------------------------------------------------------------------------- From 7de3da50e7329a1bd961a3f3ea40c2a1adb6ecfb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 11:48:21 +0100 Subject: [PATCH 06/14] include weighting operation into Gamma operator avoids point-wise multiplication. --- src/grid/spectral_utilities.f90 | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 9e5282de0..30ea11d26 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -379,6 +379,7 @@ end subroutine spectral_utilities_init subroutine utilities_updateGamma(C) real(pReal), intent(in), dimension(3,3,3,3) :: C !< input stiffness to store as reference stiffness + complex(pReal), dimension(3,3) :: temp33_cmplx, xiDyad_cmplx real(pReal), dimension(6,6) :: A, A_inv integer :: & @@ -386,7 +387,8 @@ subroutine utilities_updateGamma(C) l, m, n, o logical :: err - C_ref = C + + C_ref = C/wgt if (.not. num%memory_efficient) then gamma_hat = cmplx(0.0_pReal,0.0_pReal,pReal) ! for the singular point and any non invertible A @@ -586,9 +588,9 @@ function utilities_fourierGammaConvolution(field, fieldAim) result(gammaField) !$OMP END PARALLEL DO end if memoryEfficient - if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim/wgt,0.0_pReal,pReal) + if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim,0.0_pReal,pReal) call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) - gammaField = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3)*wgt + gammaField = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) end function utilities_fourierGammaConvolution @@ -634,8 +636,9 @@ real(pReal) function utilities_divergenceRMS(tensorField) print'(/,1x,a)', '... calculating divergence ................................................' flush(IO_STDOUT) - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = tensorField - call utilities_FFTtensorforward() + tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = tensorField + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal) @@ -684,8 +687,9 @@ real(pReal) function utilities_curlRMS(tensorField) print'(/,1x,a)', '... calculating curl ......................................................' flush(IO_STDOUT) - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = tensorField - call utilities_FFTtensorforward() + tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = tensorField + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) rescaledGeom = cmplx(geomSize/scaledGeomSize,0.0_pReal,pReal) @@ -737,7 +741,7 @@ end function utilities_curlRMS !-------------------------------------------------------------------------------------------------- -!> @brief calculates mask compliance tensor used to adjust F to fullfill stress BC +!> @brief Calculate masked compliance tensor used to adjust F to fullfill stress BC. !-------------------------------------------------------------------------------------------------- function utilities_maskedCompliance(rot_BC,mask_stress,C) From ad3c18b29b25b19d5fc5ed50eb87b6e18cfc6179 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 12:24:16 +0100 Subject: [PATCH 07/14] avoid use of global variables --- src/grid/grid_damage_spectral.f90 | 10 +++---- src/grid/grid_mech_spectral_basic.f90 | 2 +- src/grid/grid_mech_spectral_polarisation.f90 | 2 +- src/grid/grid_thermal_spectral.f90 | 10 +++---- src/grid/spectral_utilities.f90 | 29 ++++++++++++-------- 5 files changed, 27 insertions(+), 26 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index f09992f55..6e88ea55a 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -318,16 +318,14 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & - + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & - + mu_ref*phi_current(i,j,k) + r(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + + mu_ref*phi_current(i,j,k) end do; end do; end do - call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) - !-------------------------------------------------------------------------------------------------- ! constructing residual - r = max(min(scalarField_real(1:cells(1),1:cells(2),1:cells3),phi_lastInc),num%residualStiffness) & + r = max(min(utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t),phi_lastInc),num%residualStiffness) & - phi_current err_PETSc = 0 diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 1ec6d334c..cab3d4c7b 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -546,7 +546,7 @@ subroutine formResidual(in, F, & F_aim = F_aim - deltaF_aim err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask))) - r = utilities_fourierGammaConvolution(r,params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier + r = utilities_GammaConvolution(r,params%rotation_BC%rotate(deltaF_aim,active=.true.)) end subroutine formResidual diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 235bb873e..90e6d346f 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -621,7 +621,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! constructing residual r_F_tau = num%beta*F & - - utilities_fourierGammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) + - utilities_GammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index c6b190254..0f970c8da 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -341,17 +341,15 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & - + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & - + mu_ref*T_current(i,j,k) + r(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & + + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + + mu_ref*T_current(i,j,k) end do; end do; end do - call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t) - !-------------------------------------------------------------------------------------------------- ! constructing residual r = T_current & - - scalarField_real(1:cells(1),1:cells(2),1:cells3) + - utilities_GreenConvolution(r, K_ref, mu_ref, params%Delta_t) err_PETSc = 0 end subroutine formResidual diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 30ea11d26..7bccfc0b0 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -116,8 +116,8 @@ module spectral_utilities public :: & spectral_utilities_init, & utilities_updateGamma, & - utilities_fourierGammaConvolution, & - utilities_fourierGreenConvolution, & + utilities_GammaConvolution, & + utilities_GreenConvolution, & utilities_divergenceRMS, & utilities_curlRMS, & utilities_fourierScalarGradient, & @@ -507,7 +507,7 @@ end subroutine utilities_FFTvectorBackward !-------------------------------------------------------------------------------------------------- !> @brief doing convolution gamma_hat * field_real, ensuring that average value = fieldAim !-------------------------------------------------------------------------------------------------- -function utilities_fourierGammaConvolution(field, fieldAim) result(gammaField) +function utilities_GammaConvolution(field, fieldAim) result(gammaField) real(pReal), intent(in), dimension(3,3,cells(1),cells(2),cells3) :: field real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution @@ -528,8 +528,6 @@ function utilities_fourierGammaConvolution(field, fieldAim) result(gammaField) tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = field call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) - !-------------------------------------------------------------------------------------------------- -! do the actual spectral method calculation (mechanical equilibrium) memoryEfficient: if (num%memory_efficient) then !$OMP PARALLEL DO PRIVATE(l,m,n,o,temp33_cmplx,xiDyad_cmplx,A,A_inv,err,gamma_hat) do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red @@ -592,22 +590,27 @@ function utilities_fourierGammaConvolution(field, fieldAim) result(gammaField) call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) gammaField = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) -end function utilities_fourierGammaConvolution +end function utilities_GammaConvolution !-------------------------------------------------------------------------------------------------- -!> @brief doing convolution DamageGreenOp_hat * field_real +!> @brief Convolution of Greens' operator for damage/thermal. !-------------------------------------------------------------------------------------------------- -subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t) +function utilities_GreenConvolution(field, D_ref, mu_ref, Delta_t) result(greenField) + real(pReal), intent(in), dimension(cells(1),cells(2),cells3) :: field real(pReal), dimension(3,3), intent(in) :: D_ref real(pReal), intent(in) :: mu_ref, Delta_t + real(pReal), dimension(cells(1),cells(2),cells3) :: greenField + complex(pReal) :: GreenOp_hat integer :: i, j, k -!-------------------------------------------------------------------------------------------------- -! do the actual spectral method calculation - call utilities_FFTscalarForward() + + scalarField_real(cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + scalarField_real(1:cells(1), 1:cells(2),1:cells3) = field + call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) + !$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) & @@ -618,7 +621,9 @@ subroutine utilities_fourierGreenConvolution(D_ref, mu_ref, Delta_t) !$OMP END PARALLEL DO call utilities_FFTscalarBackward() -end subroutine utilities_fourierGreenConvolution + greenField = scalarField_real(1:cells(1),1:cells(2),1:cells3) + +end function utilities_GreenConvolution !-------------------------------------------------------------------------------------------------- From f22ff8fa257866cc181bc2a6210309fca4e397b5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 12:36:56 +0100 Subject: [PATCH 08/14] avoid state-changing functions requires explicit padding, i.e. a little bit of code duplication --- src/grid/spectral_utilities.f90 | 103 ++++++++------------------------ 1 file changed, 24 insertions(+), 79 deletions(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 7bccfc0b0..e353b7fb6 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -430,68 +430,6 @@ subroutine utilities_updateGamma(C) end subroutine utilities_updateGamma -!-------------------------------------------------------------------------------------------------- -!> @brief forward FFT of data in field_real to field_fourier -!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set -! to 0.0 -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTtensorForward() - - tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) - -end subroutine utilities_FFTtensorForward - - -!-------------------------------------------------------------------------------------------------- -!> @brief backward FFT of data in field_fourier to field_real -!> @details Does an weighted inverse FFT transform from complex to real -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTtensorBackward() - - call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) - tensorField_real = tensorField_real * wgt ! normalize the result by number of elements - -end subroutine utilities_FFTtensorBackward - -!-------------------------------------------------------------------------------------------------- -!> @brief forward FFT of data in scalarField_real to scalarField_fourier -!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set -! to 0.0 -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTscalarForward() - - scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) - -end subroutine utilities_FFTscalarForward - - -!-------------------------------------------------------------------------------------------------- -!> @brief backward FFT of data in scalarField_fourier to scalarField_real -!> @details Does an weighted inverse FFT transform from complex to real -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTscalarBackward() - - call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) - scalarField_real = scalarField_real * wgt ! normalize the result by number of elements - -end subroutine utilities_FFTscalarBackward - - -!-------------------------------------------------------------------------------------------------- -!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed -!> @details Does an unweighted FFT transform from real to complex. Extra padding entries are set -! to 0.0 -!-------------------------------------------------------------------------------------------------- -subroutine utilities_FFTvectorForward() - - vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) - -end subroutine utilities_FFTvectorForward - - !-------------------------------------------------------------------------------------------------- !> @brief backward FFT of data in field_fourier to field_real !> @details Does an weighted inverse FFT transform from complex to real @@ -619,7 +557,8 @@ function utilities_GreenConvolution(field, D_ref, mu_ref, Delta_t) result(greenF 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() + call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) + scalarField_real = scalarField_real * wgt ! normalize the result by number of elements greenField = scalarField_real(1:cells(1),1:cells(2),1:cells3) @@ -823,11 +762,13 @@ subroutine utilities_fourierScalarGradient() integer :: i, j, k - call utilities_FFTscalarForward() + scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal + call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) 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() + call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) + vectorField_real = vectorField_real * wgt ! normalize the result by number of elements end subroutine utilities_fourierScalarGradient @@ -837,10 +778,12 @@ end subroutine utilities_fourierScalarGradient !-------------------------------------------------------------------------------------------------- subroutine utilities_fourierVectorDivergence() - call utilities_FFTvectorForward() + vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal + call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) 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() + call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) + scalarField_real = scalarField_real * wgt ! normalize the result by number of elements end subroutine utilities_fourierVectorDivergence @@ -1075,8 +1018,9 @@ subroutine utilities_updateCoords(F) step = geomSize/real(cells, pReal) !-------------------------------------------------------------------------------------------------- ! integration in Fourier space to get fluctuations of cell center discplacements - tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) = F - call utilities_FFTtensorForward() + tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = F + tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) !$OMP PARALLEL DO do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red @@ -1089,7 +1033,8 @@ subroutine utilities_updateCoords(F) end do; end do; end do !$OMP END PARALLEL DO - call utilities_FFTvectorBackward() + call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) + vectorField_real = vectorField_real * wgt ! normalize the result by number of elements !-------------------------------------------------------------------------------------------------- ! average F @@ -1183,38 +1128,38 @@ subroutine selfTest() call random_number(tensorField_real) tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal tensorField_real_ = tensorField_real - call utilities_FFTtensorForward() + call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) if (worldsize==1) then if (any(dNeq(sum(sum(sum(tensorField_real_,dim=5),dim=4),dim=3)/tensorField_fourier(:,:,1,1,1)%re,1.0_pReal,1.0e-12_pReal))) & error stop 'tensorField avg' endif - call utilities_FFTtensorBackward() + call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - if (maxval(abs(tensorField_real_ - tensorField_real))>5.0e-15_pReal) error stop 'tensorField' + if (maxval(abs(tensorField_real_ - tensorField_real*wgt))>5.0e-15_pReal) error stop 'tensorField' call random_number(vectorField_real) vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal vectorField_real_ = vectorField_real - call utilities_FFTvectorForward() + call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) if (worldsize==1) then if (any(dNeq(sum(sum(sum(vectorField_real_,dim=4),dim=3),dim=2)/vectorField_fourier(:,1,1,1)%re,1.0_pReal,1.0e-12_pReal))) & error stop 'vector avg' endif - call utilities_FFTvectorBackward() + call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - if (maxval(abs(vectorField_real_ - vectorField_real))>5.0e-15_pReal) error stop 'vectorField' + if (maxval(abs(vectorField_real_ - vectorField_real*wgt))>5.0e-15_pReal) error stop 'vectorField' call random_number(scalarField_real) scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal scalarField_real_ = scalarField_real - call utilities_FFTscalarForward() + call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) if (worldsize==1) then if (dNeq(sum(sum(sum(scalarField_real_,dim=3),dim=2),dim=1)/scalarField_fourier(1,1,1)%re,1.0_pReal,1.0e-12_pReal)) & error stop 'scalar avg' endif - call utilities_FFTscalarBackward() + call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal - if (maxval(abs(scalarField_real_ - scalarField_real))>5.0e-15_pReal) error stop 'scalarField' + if (maxval(abs(scalarField_real_ - scalarField_real*wgt))>5.0e-15_pReal) error stop 'scalarField' end subroutine selfTest From 6db3b72c89c3affd859599eccb5d0349da62e8a4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 13:05:12 +0100 Subject: [PATCH 09/14] avoid global variables extra memory (one vector field) required --- src/grid/grid_damage_spectral.f90 | 10 +++--- src/grid/grid_thermal_spectral.f90 | 10 +++--- src/grid/spectral_utilities.f90 | 53 +++++++++++++++++------------- 3 files changed, 41 insertions(+), 32 deletions(-) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 6e88ea55a..3e514b960 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -302,23 +302,23 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) PetscErrorCode, intent(out) :: err_PETSc integer :: i, j, k, ce + real(pReal), dimension(3,cells(1),cells(2),cells3) :: vectorField phi_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_real(1:cells(1),1:cells(2),1:cells3) = phi_current - call utilities_fourierScalarGradient() + vectorField = utilities_ScalarGradient(phi_current) 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)) + vectorField(1:3,i,j,k) = matmul(homogenization_K_phi(ce) - K_ref, vectorField(1:3,i,j,k)) end do; end do; end do - call utilities_fourierVectorDivergence() + r = utilities_VectorDivergence(vectorField) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - r(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + mu_ref*phi_current(i,j,k) end do; end do; end do diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 0f970c8da..b1dde568f 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -325,23 +325,23 @@ subroutine formResidual(in,x_scal,r,dummy,err_PETSc) PetscErrorCode, intent(out) :: err_PETSc integer :: i, j, k, ce + real(pReal), dimension(3,cells(1),cells(2),cells3) :: vectorField T_current = x_scal !-------------------------------------------------------------------------------------------------- ! evaluate polarization field - scalarField_real(1:cells(1),1:cells(2),1:cells3) = T_current - call utilities_fourierScalarGradient() + vectorField = utilities_ScalarGradient(T_current) 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)) + vectorField(1:3,i,j,k) = matmul(homogenization_K_T(ce) - K_ref, vectorField(1:3,i,j,k)) end do; end do; end do - call utilities_fourierVectorDivergence() + r = utilities_VectorDivergence(vectorField) ce = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1,cells(1) ce = ce + 1 - r(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & + r(i,j,k) = params%Delta_t*(r(i,j,k) + homogenization_f_T(ce)) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + mu_ref*T_current(i,j,k) end do; end do; end do diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index e353b7fb6..aebb7fa15 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -42,16 +42,16 @@ module spectral_utilities !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW - real(C_DOUBLE), dimension(:,:,:,:,:), pointer :: tensorField_real !< tensor field in real space - real(C_DOUBLE), public, dimension(:,:,:,:), pointer :: vectorField_real !< vector field in real space - real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field in real space - complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:,:), pointer :: tensorField_fourier !< tensor field in Fourier space - complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field in Fourier space - complex(C_DOUBLE_COMPLEX), dimension(:,:,:), pointer :: scalarField_fourier !< scalar field in Fourier space - complex(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method - complex(pReal), dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives - complex(pReal), dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives - real(pReal), dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness + real(C_DOUBLE), dimension(:,:,:,:,:), pointer :: tensorField_real !< tensor field in real space + real(C_DOUBLE), dimension(:,:,:,:), pointer :: vectorField_real !< vector field in real space + real(C_DOUBLE), dimension(:,:,:), pointer :: scalarField_real !< scalar field in real space + complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:,:), pointer :: tensorField_fourier !< tensor field in Fourier space + complex(C_DOUBLE_COMPLEX), dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field in Fourier space + complex(C_DOUBLE_COMPLEX), dimension(:,:,:), pointer :: scalarField_fourier !< scalar field in Fourier space + complex(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method + complex(pReal), dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives + complex(pReal), dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives + real(pReal), dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness !-------------------------------------------------------------------------------------------------- @@ -120,8 +120,8 @@ module spectral_utilities utilities_GreenConvolution, & utilities_divergenceRMS, & utilities_curlRMS, & - utilities_fourierScalarGradient, & - utilities_fourierVectorDivergence, & + utilities_ScalarGradient, & + utilities_VectorDivergence, & utilities_maskedCompliance, & utilities_constitutiveResponse, & utilities_calculateRate, & @@ -755,37 +755,46 @@ end function utilities_maskedCompliance !-------------------------------------------------------------------------------------------------- -!> @brief calculate scalar gradient in fourier field +!> @brief Calculate gradient of scalar field. !-------------------------------------------------------------------------------------------------- -subroutine utilities_fourierScalarGradient() +function utilities_ScalarGradient(field) result(grad) + + real(pReal), intent(in), dimension( cells(1),cells(2),cells3) :: field + real(pReal), dimension(3,cells(1),cells(2),cells3) :: grad integer :: i, j, k - scalarField_real(cells(1)+1:cells1Red*2,:,:) = 0.0_pReal + scalarField_real(cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + scalarField_real(1:cells(1), 1:cells(2),1:cells3) = field call fftw_mpi_execute_dft_r2c(planScalarForth,scalarField_real,scalarField_fourier) 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 fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) - vectorField_real = vectorField_real * wgt ! normalize the result by number of elements + grad = vectorField_real(1:3,1:cells(1),1:cells(2),1:cells3)*wgt -end subroutine utilities_fourierScalarGradient +end function utilities_ScalarGradient !-------------------------------------------------------------------------------------------------- -!> @brief calculate vector divergence in fourier field +!> @brief Calculate divergence of vector field. !-------------------------------------------------------------------------------------------------- -subroutine utilities_fourierVectorDivergence() +function utilities_VectorDivergence(field) result(div) - vectorField_real(1:3,cells(1)+1:cells1Red*2,:,:) = 0.0_pReal + real(pReal), intent(in), dimension(3,cells(1),cells(2),cells3) :: field + real(pReal), dimension( cells(1),cells(2),cells3) :: div + + + vectorField_real(1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal + vectorField_real(1:3,1:cells(1), 1:cells(2),1:cells3) = field call fftw_mpi_execute_dft_r2c(planVectorForth,vectorField_real,vectorField_fourier) 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 fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) - scalarField_real = scalarField_real * wgt ! normalize the result by number of elements + div = scalarField_real(1:cells(1),1:cells(2),1:cells3)*wgt -end subroutine utilities_fourierVectorDivergence +end function utilities_VectorDivergence !-------------------------------------------------------------------------------------------------- From eb226d237f2070858aaec46b7da24c605e9057e0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 13:25:21 +0100 Subject: [PATCH 10/14] better readable --- src/grid/spectral_utilities.f90 | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index aebb7fa15..3a3582f7d 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -300,8 +300,8 @@ subroutine spectral_utilities_init() !-------------------------------------------------------------------------------------------------- ! allocation - allocate (xi1st (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension - allocate (xi2nd (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension + allocate (xi1st (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for first derivatives, only half the size for first dimension + allocate (xi2nd (3,cells1Red,cells(3),cells2),source = cmplx(0.0_pReal,0.0_pReal,pReal)) ! frequencies for second derivatives, only half the size for first dimension !-------------------------------------------------------------------------------------------------- ! tensor MPI fftw plans @@ -525,6 +525,7 @@ function utilities_GammaConvolution(field, fieldAim) result(gammaField) end if memoryEfficient if (cells3Offset == 0) tensorField_fourier(1:3,1:3,1,1,1) = cmplx(fieldAim,0.0_pReal,pReal) + call fftw_mpi_execute_dft_c2r(planTensorBack,tensorField_fourier,tensorField_real) gammaField = tensorField_real(1:3,1:3,1:cells(1),1:cells(2),1:cells3) @@ -557,10 +558,9 @@ function utilities_GreenConvolution(field, D_ref, mu_ref, Delta_t) result(greenF scalarField_fourier(i,k,j) = scalarField_fourier(i,k,j)*GreenOp_hat end do; end do; end do !$OMP END PARALLEL DO - call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) - scalarField_real = scalarField_real * wgt ! normalize the result by number of elements - greenField = scalarField_real(1:cells(1),1:cells(2),1:cells3) + call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) + greenField = scalarField_real(1:cells(1),1:cells(2),1:cells3) * wgt end function utilities_GreenConvolution @@ -1025,12 +1025,19 @@ subroutine utilities_updateCoords(F) step = geomSize/real(cells, pReal) - !-------------------------------------------------------------------------------------------------- - ! integration in Fourier space to get fluctuations of cell center discplacements + tensorField_real(1:3,1:3,1:cells(1), 1:cells(2),1:cells3) = F tensorField_real(1:3,1:3,cells(1)+1:cells1Red*2,1:cells(2),1:cells3) = 0.0_pReal call fftw_mpi_execute_dft_r2c(planTensorForth,tensorField_real,tensorField_fourier) + !-------------------------------------------------------------------------------------------------- + ! average F + if (cells3Offset == 0) Favg = tensorField_fourier(1:3,1:3,1,1,1)%re*wgt + call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + + !-------------------------------------------------------------------------------------------------- + ! integration in Fourier space to get fluctuations of cell center discplacements !$OMP PARALLEL DO do j = 1, cells2; do k = 1, cells(3); do i = 1, cells1Red if (any([i,j+cells2Offset,k] /= 1)) then @@ -1045,12 +1052,6 @@ subroutine utilities_updateCoords(F) call fftw_mpi_execute_dft_c2r(planVectorBack,vectorField_fourier,vectorField_real) vectorField_real = vectorField_real * wgt ! normalize the result by number of elements - !-------------------------------------------------------------------------------------------------- - ! average F - if (cells3Offset == 0) Favg = tensorField_fourier(1:3,1:3,1,1,1)%re*wgt - call MPI_Bcast(Favg,9_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - !-------------------------------------------------------------------------------------------------- ! pad cell center fluctuations along z-direction (needed when running MPI simulation) IPfluct_padded(1:3,1:cells(1),1:cells(2),2:cells3+1) = vectorField_real(1:3,1:cells(1),1:cells(2),1:cells3) From 9b80ff623b687a3be84b2876982ede91bd20004c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Nov 2022 13:35:55 +0100 Subject: [PATCH 11/14] faster operation explicit weighting not needed --- src/grid/spectral_utilities.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 3a3582f7d..d5ec36702 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -552,7 +552,7 @@ function utilities_GreenConvolution(field, D_ref, mu_ref, Delta_t) result(greenF !$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) & + GreenOp_hat = cmplx(wgt,0.0_pReal,pReal) & / (cmplx(mu_ref,0.0_pReal,pReal) + cmplx(Delta_t,0.0_pReal,pReal) & * sum(conjg(xi1st(1:3,i,k,j))* matmul(cmplx(D_ref,0.0_pReal,pReal),xi1st(1:3,i,k,j)))) scalarField_fourier(i,k,j) = scalarField_fourier(i,k,j)*GreenOp_hat @@ -560,7 +560,7 @@ function utilities_GreenConvolution(field, D_ref, mu_ref, Delta_t) result(greenF !$OMP END PARALLEL DO call fftw_mpi_execute_dft_c2r(planScalarBack,scalarField_fourier,scalarField_real) - greenField = scalarField_real(1:cells(1),1:cells(2),1:cells3) * wgt + greenField = scalarField_real(1:cells(1),1:cells(2),1:cells3) end function utilities_GreenConvolution From 34fb7e921a4754fa0eb8b5d1953e7c1adaff7504 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Nov 2022 12:58:50 +0100 Subject: [PATCH 12/14] use self-documenting code the comments did not anything that was not clear from the variable/function names --- src/grid/grid_mech_spectral_basic.f90 | 20 ++++++--------- src/grid/grid_mech_spectral_polarisation.f90 | 27 ++++++-------------- 2 files changed, 16 insertions(+), 31 deletions(-) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index cab3d4c7b..c570e94e7 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -519,8 +519,6 @@ subroutine formResidual(in, F, & if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment -!-------------------------------------------------------------------------------------------------- -! begin of new iteration newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax @@ -531,17 +529,15 @@ subroutine formResidual(in, F, & flush(IO_STDOUT) end if newIteration -!-------------------------------------------------------------------------------------------------- -! evaluate constitutive response - call utilities_constitutiveResponse(r, & ! residuum gets field of first PK stress (to save memory) - P_av,C_volAvg,C_minMaxAvg, & - F,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' - err_div = utilities_divergenceRMS(r) + associate (P => r) + call utilities_constitutiveResponse(P, & + P_av,C_volAvg,C_minMaxAvg, & + F,params%Delta_t,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' + err_div = utilities_divergenceRMS(P) + end associate -!-------------------------------------------------------------------------------------------------- -! stress BC handling deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - deltaF_aim err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask))) diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 90e6d346f..e2ec8488c 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -597,8 +597,6 @@ subroutine formResidual(in, FandF_tau, & if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment -!-------------------------------------------------------------------------------------------------- -! begin of new iteration newIteration: if (totalIter <= PETScIter) then totalIter = totalIter + 1 print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax @@ -609,38 +607,29 @@ subroutine formResidual(in, FandF_tau, & flush(IO_STDOUT) end if newIteration -!-------------------------------------------------------------------------------------------------- -! do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) r_F_tau(1:3,1:3,i,j,k) = & num%beta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& num%alpha*matmul(F(1:3,1:3,i,j,k), & 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 - -!-------------------------------------------------------------------------------------------------- -! constructing residual r_F_tau = num%beta*F & - utilities_GammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) -!-------------------------------------------------------------------------------------------------- -! evaluate constitutive response - call utilities_constitutiveResponse(r_F, & ! "residuum" gets field of first PK stress (to save memory) - P_av,C_volAvg,C_minMaxAvg, & - F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC) - call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) - err_div = utilities_divergenceRMS(r_F) - err_curl = utilities_curlRMS(F) + associate (P => r_F) + call utilities_constitutiveResponse(P, & + P_av,C_volAvg,C_minMaxAvg, & + F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC) + call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) + err_div = utilities_divergenceRMS(P) + err_curl = utilities_curlRMS(F) + end associate -!-------------------------------------------------------------------------------------------------- -! stress BC handling F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), & P_av-P_aim, & params%stress_mask))) -!-------------------------------------------------------------------------------------------------- -! constructing residual e = 0 do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) e = e + 1 From 2173c9e4992c88ecc2999d077086e2ef506b022b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Nov 2022 13:36:03 +0100 Subject: [PATCH 13/14] undo weighting needed for restart --- src/grid/spectral_utilities.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index d5ec36702..54dbf8c43 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -1118,7 +1118,7 @@ subroutine utilities_saveReferenceStiffness() open(newunit=fileUnit, file=getSolverJobName()//'.C_ref',& status='replace',access='stream',action='write',iostat=ierr) if (ierr /=0) call IO_error(100,ext_msg='could not open file '//getSolverJobName()//'.C_ref') - write(fileUnit) C_ref + write(fileUnit) C_ref*wgt close(fileUnit) end if From cad4cbc5d24ba749749f47339f9ae2e5867d7c62 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Nov 2022 21:04:42 +0100 Subject: [PATCH 14/14] circument bug in gfortran associate to strided pointer seems to cause trouble --- src/grid/DAMASK_grid.f90 | 2 -- src/grid/grid_mech_spectral_polarisation.f90 | 35 ++++++++++++++------ 2 files changed, 25 insertions(+), 12 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 1b785ca3d..c9dea0166 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -106,8 +106,6 @@ program DAMASK_grid external :: & quit - class(tNode), pointer :: & - tmp type(tDict), pointer :: & config_load, & num_grid, & diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index e2ec8488c..46ff04adb 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -616,29 +616,44 @@ subroutine formResidual(in, FandF_tau, & r_F_tau = num%beta*F & - utilities_GammaConvolution(r_F_tau,params%rotation_BC%rotate(num%beta*F_aim,active=.true.)) + err_curl = utilities_curlRMS(F) + +#ifdef __GFORTRAN__ + call utilities_constitutiveResponse(r_F, & +#else associate (P => r_F) call utilities_constitutiveResponse(P, & +#endif P_av,C_volAvg,C_minMaxAvg, & F - r_F_tau/num%beta,params%Delta_t,params%rotation_BC) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) +#ifdef __GFORTRAN__ + err_div = utilities_divergenceRMS(r_F) +#else err_div = utilities_divergenceRMS(P) - err_curl = utilities_curlRMS(F) +#endif + e = 0 + do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) + e = e + 1 + r_F(1:3,1:3,i,j,k) = & + math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & +#ifdef __GFORTRAN__ + r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & +#else + P(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & +#endif + math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & + + r_F_tau(1:3,1:3,i,j,k) + end do; end do; end do +#ifndef __GFORTRAN__ end associate +#endif F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), & P_av-P_aim, & params%stress_mask))) - e = 0 - do k = 1, cells3; do j = 1, cells(2); do i = 1, cells(1) - e = e + 1 - r_F(1:3,1:3,i,j,k) = & - math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,e) + C_scale), & - r_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & - math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & - + r_F_tau(1:3,1:3,i,j,k) - end do; end do; end do end subroutine formResidual