diff --git a/code/DAMASK_spectral_SolverAL.f90 b/code/DAMASK_spectral_SolverAL.f90 index 5bbe9130e..68274fb4e 100644 --- a/code/DAMASK_spectral_SolverAL.f90 +++ b/code/DAMASK_spectral_SolverAL.f90 @@ -55,7 +55,7 @@ module DAMASK_spectral_SolverAL !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc, P + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates real(pReal), private, dimension(:,:,:), allocatable :: temperature @@ -113,11 +113,12 @@ module DAMASK_spectral_SolverAL implicit none integer(pInt) :: i,j,k + real(pReal), dimension(:,:,:,:,:), allocatable :: P PetscErrorCode :: ierr_psc PetscObject :: dummy PetscMPIInt :: rank - PetscScalar, pointer :: xx_psc(:,:,:,:) + PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:), F_lambda(:,:,:,:) call Utilities_init() @@ -153,27 +154,23 @@ module DAMASK_spectral_SolverAL !-------------------------------------------------------------------------------------------------- ! init fields call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc) - + F => xx_psc(0:8,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) if (restartInc == 1_pInt) then ! no deformation (no restart) - F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity + F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity F_lambda_lastInc = F_lastInc - - xx_psc(0:8,:,:,:) = reshape(F_lastInc,[9,res(1),res(2),res(3)]) - - xx_psc(9:17,:,:,:) = xx_psc(0:8,:,:,:) - - call flush(6) + F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + F_lambda = F do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - geomdim/real(2_pInt*res,pReal) enddo; enddo; enddo - print*, 'init done12' elseif (restartInc > 1_pInt) then ! using old values from file if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& restartInc - 1_pInt,' from file' call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) xx_psc(0:8,:,:,:) + read (777,rec=1) F close (777) call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& trim(getSolverJobName()),size(F_lastInc)) @@ -181,7 +178,7 @@ module DAMASK_spectral_SolverAL close (777) call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',& trim(getSolverJobName()),size(F_lambda_lastInc)) - read (777,rec=1) xx_psc(9:17,:,:,:) + read (777,rec=1) F_lambda close (777) call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',& trim(getSolverJobName()),size(F_lambda_lastInc)) @@ -197,12 +194,9 @@ module DAMASK_spectral_SolverAL coordinates = 0.0 ! change it later!!! endif - call Utilities_constitutiveResponse(coordinates,reshape(xx_psc(0:8,:,:,:),shape(F_lastInc)),& - reshape(xx_psc(0:8,:,:,:),shape(F_lastInc)),& - temperature,0.0_pReal,P,C,P_av,.false.,math_I3) - print*, 'const response' + call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc) - print*, 'restored' + !-------------------------------------------------------------------------------------------------- ! reference stiffness if (restartInc == 1_pInt) then @@ -253,19 +247,18 @@ module DAMASK_spectral_SolverAL real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode type(boundaryCondition), intent(in) :: P_BC,F_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC - SNESConvergedReason reason real(pReal), dimension(3,3) :: deltaF_aim, & F_aim_lab !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. real(pReal), dimension(3,3) :: temp33_Real - integer(pInt) :: ctr, i, j, k !-------------------------------------------------------------------------------------------------- ! - PetscScalar, pointer :: xx_psc(:,:,:,:) - PetscErrorCode ierr_psc - + PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:), F_lambda(:,:,:,:) + PetscErrorCode ierr_psc + SNESConvergedReason reason + !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver if (restartWrite) then @@ -296,10 +289,10 @@ module DAMASK_spectral_SolverAL ! update local deformation gradient and coordinates deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc) - call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc, & - xx_psc(1:9,1:res(1),1:res(2),1:res(3))) - !call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lambda_lastInc,& - ! reshape(xx_psc(10:18,1:res(1),1:res(2),1:res(3)),shape(F_lambda_lastInc))) + F => xx_psc(0:8,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) + call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc,F) + call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lambda_lastInc,F_lambda) call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc) call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) @@ -345,15 +338,22 @@ module DAMASK_spectral_SolverAL implicit none integer(pInt) :: i,j,k,l, ctr - real(pReal), dimension (3,3) :: temp33_real + real(pReal), dimension(3,3) :: temp33_Real DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE) - PetscScalar :: x_scal(in(DMDA_LOCAL_INFO_DOF),XG_RANGE,YG_RANGE,ZG_RANGE) - PetscScalar :: f_scal(in(DMDA_LOCAL_INFO_DOF),X_RANGE,Y_RANGE,Z_RANGE) + PetscScalar, target :: x_scal(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE) + PetscScalar, target :: f_scal(3,3,2,X_RANGE,Y_RANGE,Z_RANGE) + PetscScalar, pointer :: F(:,:,:,:,:), F_lambda(:,:,:,:,:) + PetscScalar, pointer :: residual_F(:,:,:,:,:), residual_F_lambda(:,:,:,:,:) PetscInt :: iter, nfuncs PetscObject :: dummy PetscErrorCode :: ierr_psc + F => x_scal(:,:,1,:,:,:) + F_lambda => x_scal(:,:,2,:,:,:) + residual_F => f_scal(:,:,1,:,:,:) + residual_F_lambda => f_scal(:,:,2,:,:,:) + call SNESGetNumberFunctionEvals(snes,nfuncs,ierr_psc) call SNESGetIterationNumber(snes,iter,ierr_psc) @@ -367,9 +367,8 @@ module DAMASK_spectral_SolverAL !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(coordinates,F_lastInc,reshape(x_scal(1:9,1:res(1),1:res(2),1:res(3)),shape(F_lastInc)),& - temperature,params%timeinc,& - P,C,P_av,ForwardData,params%rotation_BC) + call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, & + residual_F,C,P_av,ForwardData,params%rotation_BC) ForwardData = .False. !-------------------------------------------------------------------------------------------------- @@ -381,37 +380,22 @@ module DAMASK_spectral_SolverAL !-------------------------------------------------------------------------------------------------- ! doing Fourier transform field_real = 0.0_pReal - - do j = 1_pInt, 3_pInt; do i = 1_pInt, 3_pInt - ctr = 1_pInt - do k = 1_pInt, 3_pInt; do l = 1_pInt, 3_pInt - field_real(1:res(1),1:res(2),1:res(3),i,j) = field_real(1:res(1),1:res(2),1:res(3),i,j) + & - C_scale(i,j,l,k)*(x_scal(9+ctr,1:res(1),1:res(2),1:res(3)) - & - x_scal(ctr,1:res(1),1:res(2),1:res(3))) - ! P_temp(i,j,1:res(1),1:res(2),1:res(3)) = P_temp(i,j,1:res(1),1:res(2),1:res(3)) + & - ! S_scale(i,j,l,k)*P(l,k,1:res(1),1:res(2),1:res(3)) - ctr = ctr + 1_pInt - enddo; enddo - enddo; enddo - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - P(1:3,1:3,i,j,k) = math_mul3333xx33(S_scale,P(1:3,1:3,i,j,k)) + math_I3 + temp33_Real = math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) + math_I3 + residual_F(1:3,1:3,i,j,k) = temp33_Real + field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k)-F(1:3,1:3,i,j,k)) enddo; enddo; enddo call Utilities_forwardFFT() call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC)) call Utilities_backwardFFT() - f_scal(1:9,1:res(1),1:res(2),1:res(3)) = reshape(P,[9,res(1),res(2),res(3)]) - x_scal(10:18,1:res(1),1:res(2),1:res(3)) +& - x_scal(1:9,1:res(1),1:res(2),1:res(3)) - & - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),[9,res(1),res(2),res(3)],order=[2,3,4,1]) - f_scal(10:18,1:res(1),1:res(2),1:res(3)) = x_scal(1:9,1:res(1),1:res(2),1:res(3)) - & - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),& - [9,res(1),res(2),res(3)],order=[2,3,4,1]) + residual_F_lambda = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),& + [3,3,res(1),res(2),res(3)],order=[3,4,5,1,2]) + residual_F = residual_F - F_lambda + residual_F_lambda - err_f = wgt*sqrt(sum(f_scal(10:18,1:res(1),1:res(2),1:res(3))**2.0_pReal)) - err_p = wgt*sqrt(sum((f_scal( 1:9 ,1:res(1),1:res(2),1:res(3)) & - -f_scal(10:18,1:res(1),1:res(2),1:res(3)))**2.0_pReal)) + err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal)) + err_p = wgt*sqrt(sum((residual_F - residual_F_lambda)**2.0_pReal)) end subroutine AL_formResidual @@ -470,6 +454,7 @@ module DAMASK_spectral_SolverAL call VecDestroy(solution_vec,ierr_psc) call SNESDestroy(snes,ierr_psc) + call DMDestroy(da,ierr_psc) call PetscFinalize(ierr_psc) call Utilities_destroy() diff --git a/code/DAMASK_spectral_Utilities.f90 b/code/DAMASK_spectral_Utilities.f90 index 404a7cbeb..8b3093d98 100644 --- a/code/DAMASK_spectral_Utilities.f90 +++ b/code/DAMASK_spectral_Utilities.f90 @@ -357,6 +357,7 @@ subroutine Utilities_fourierConvolution(fieldAim) write(6,'(a)') '' write(6,'(a)') '... doing convolution .................' + write(6,'(a)') '' !-------------------------------------------------------------------------------------------------- ! to the actual spectral method calculation (mechanical equilibrium) @@ -536,6 +537,7 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti write(6,'(a)') '' write(6,'(a)') '... evaluating constitutive response .................' + write(6,'(a)') '' if (ForwardData) then CPFEM_mode = 1_pInt