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 !--------------------------------------------------------------------------------------------------