diff --git a/code/DAMASK_spectral_SolverAL.f90 b/code/DAMASK_spectral_SolverAL.f90 index da90d726c..5bbe9130e 100644 --- a/code/DAMASK_spectral_SolverAL.f90 +++ b/code/DAMASK_spectral_SolverAL.f90 @@ -31,6 +31,7 @@ module DAMASK_spectral_SolverAL #include #include #include +#include #include #endif @@ -48,12 +49,13 @@ module DAMASK_spectral_SolverAL !-------------------------------------------------------------------------------------------------- ! PETSc data + DM, private :: da SNES, private :: snes Vec, private :: solution_vec !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, F_lastInc, F_lambda, F_lambda_lastInc, P + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc, P real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates real(pReal), private, dimension(:,:,:), allocatable :: temperature @@ -72,8 +74,7 @@ module DAMASK_spectral_SolverAL real(pReal), private :: err_stress, err_f, err_p logical, private :: ForwardData - real(pReal), private, dimension(3,3) :: & - mask_stress = 0.0_pReal + real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal contains @@ -103,7 +104,8 @@ module DAMASK_spectral_SolverAL use mesh, only: & res, & - geomdim + geomdim, & + mesh_NcpElems use math, only: & math_invSym3333 @@ -115,7 +117,7 @@ module DAMASK_spectral_SolverAL PetscErrorCode :: ierr_psc PetscObject :: dummy PetscMPIInt :: rank - DM :: da + PetscScalar, pointer :: xx_psc(:,:,:,:) call Utilities_init() @@ -125,39 +127,61 @@ module DAMASK_spectral_SolverAL #include "compilation_info.f90" write(6,'(a)') '' - allocate (F ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (F_lambda ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (F_lambda_lastInc(res(1),res(2),res(3),3,3), source = 0.0_pReal) - allocate (P ( res(1), res(2),res(3),3,3), source = 0.0_pReal) + allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) + allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal) + allocate (P (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) - + + !-------------------------------------------------------------------------------------------------- + ! PETSc Init + call PetscInitialize(PETSC_NULL_CHARACTER,ierr_psc) + call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr_psc) + call SNESCreate(PETSC_COMM_WORLD,snes,ierr_psc) + call DMDACreate3d(PETSC_COMM_WORLD, & + DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & + DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & + 18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr_psc) + + call DMCreateGlobalVector(da,solution_vec,ierr_psc) + call DMDASetLocalFunction(da,AL_formResidual,ierr_psc) + call SNESSetDM(snes,da,ierr_psc) + call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr_psc) + call PetscOptionsInsertString(petsc_options,ierr_psc) + call SNESSetFromOptions(snes,ierr_psc) + !-------------------------------------------------------------------------------------------------- ! init fields + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc) + 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_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) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - F(i,j,k,1:3,1:3) = math_I3 - F_lastInc(i,j,k,1:3,1:3) = math_I3 - F_lambda(i,j,k,1:3,1:3) = math_I3 - F_lambda_lastInc(i,j,k,1:3,1:3) = math_I3 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)) - read (777,rec=1) F + trim(getSolverJobName()),size(F_lastInc)) + read (777,rec=1) xx_psc(0:8,:,:,:) close (777) call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& trim(getSolverJobName()),size(F_lastInc)) read (777,rec=1) F_lastInc close (777) call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',& - trim(getSolverJobName()),size(F_lambda)) - read (777,rec=1) F + trim(getSolverJobName()),size(F_lambda_lastInc)) + read (777,rec=1) xx_psc(9:17,:,:,:) close (777) call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',& trim(getSolverJobName()),size(F_lambda_lastInc)) @@ -173,9 +197,12 @@ module DAMASK_spectral_SolverAL coordinates = 0.0 ! change it later!!! endif - call Utilities_constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& - P,C,P_av,.false.,math_I3) - + 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 DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc) + print*, 'restored' !-------------------------------------------------------------------------------------------------- ! reference stiffness if (restartInc == 1_pInt) then @@ -191,25 +218,6 @@ module DAMASK_spectral_SolverAL call Utilities_updateGamma(C) C_scale = C S_scale = math_invSym3333(C) - - !-------------------------------------------------------------------------------------------------- - ! PETSc Init - call PetscInitialize(PETSC_NULL_CHARACTER,ierr_psc) - call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr_psc) - - call SNESCreate(PETSC_COMM_WORLD,snes,ierr_psc) - call DMDACreate3d(PETSC_COMM_WORLD, & - DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & - DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & - 18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr_psc) - call DMCreateGlobalVector(da,solution_vec,ierr_psc) - call DMDASetLocalFunction(da,AL_formResidual,ierr_psc) - - call SNESSetDM(snes,da,ierr_psc) - call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr_psc) - call PetscOptionsInsertString(petsc_options,ierr_psc) - call SNESSetFromOptions(snes,ierr_psc) - call DMDestroy(da,ierr_psc) end subroutine AL_init @@ -255,7 +263,7 @@ module DAMASK_spectral_SolverAL !-------------------------------------------------------------------------------------------------- ! - PetscScalar, pointer :: xx_psc(:) + PetscScalar, pointer :: xx_psc(:,:,:,:) PetscErrorCode ierr_psc !-------------------------------------------------------------------------------------------------- @@ -283,12 +291,16 @@ module DAMASK_spectral_SolverAL + deltaF_aim F_aim_lastInc = temp33_Real F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame - + !-------------------------------------------------------------------------------------------------- ! update local deformation gradient and coordinates deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) - 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 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))) + 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) !-------------------------------------------------------------------------------------------------- @@ -301,15 +313,7 @@ module DAMASK_spectral_SolverAL params%P_BC = P_BC%values params%rotation_BC = rotation_BC params%timeinc = timeinc - - ctr = 1_pInt - call VecGetArrayF90(solution_vec,xx_psc,ierr_psc) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - xx_psc(ctr:ctr+8_pInt) = reshape(F(i,j,k,1:3,1:3),[9]) - xx_psc(ctr+9_pInt:ctr+17_pInt) = reshape(F_lambda(i,j,k,1:3,1:3),[9]) - ctr = ctr + 18_pInt - enddo; enddo; enddo - call VecRestoreArrayF90(solution_vec,xx_psc,ierr_psc) + call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr_psc) call SNESGetConvergedReason(snes,reason,ierr_psc) if (reason > 0 ) AL_solution%converged = .true. @@ -340,7 +344,7 @@ module DAMASK_spectral_SolverAL implicit none - integer(pInt) :: i,j,k + integer(pInt) :: i,j,k,l, ctr real(pReal), dimension (3,3) :: temp33_real DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE) @@ -361,16 +365,10 @@ module DAMASK_spectral_SolverAL write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& math_transpose33(F_aim) - do k=in(DMDA_LOCAL_INFO_ZS)+1,in(DMDA_LOCAL_INFO_ZS)+in(DMDA_LOCAL_INFO_ZM) - do j=in(DMDA_LOCAL_INFO_YS)+1,in(DMDA_LOCAL_INFO_YS)+in(DMDA_LOCAL_INFO_YM) - do i=in(DMDA_LOCAL_INFO_XS)+1,in(DMDA_LOCAL_INFO_XS)+in(DMDA_LOCAL_INFO_XM) - F(i,j,k,1:3,1:3) = reshape(x_scal(1:9,i,j,k),[3,3]) - F_lambda(i,j,k,1:3,1:3) = reshape(x_scal(10:18,i,j,k),[3,3]) - enddo; enddo; enddo - !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc,& + 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) ForwardData = .False. @@ -383,37 +381,38 @@ 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) - field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_scale,F_lambda(i,j,k,1:3,1:3)-F(i,j,k,1:3,1:3)) + P(1:3,1:3,i,j,k) = math_mul3333xx33(S_scale,P(1:3,1:3,i,j,k)) + math_I3 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]) + + 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 = 0.0_pReal - err_p = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - temp33_real = field_real(i,j,k,1:3,1:3) - F(i,j,k,1:3,1:3) - err_f = err_f + sum(temp33_real*temp33_real) - - temp33_real = F_lambda(i,j,k,1:3,1:3) - & - (math_mul3333xx33(S_scale,P(i,j,k,1:3,1:3)) + math_I3) - err_p = err_p + sum(temp33_real*temp33_real) - enddo; enddo; enddo - err_f = wgt*sqrt(err_f) - err_p = wgt*sqrt(err_p) - - do k=in(DMDA_LOCAL_INFO_ZS)+1,in(DMDA_LOCAL_INFO_ZS)+in(DMDA_LOCAL_INFO_ZM) - do j=in(DMDA_LOCAL_INFO_YS)+1,in(DMDA_LOCAL_INFO_YS)+in(DMDA_LOCAL_INFO_YM) - do i=in(DMDA_LOCAL_INFO_XS)+1,in(DMDA_LOCAL_INFO_XS)+in(DMDA_LOCAL_INFO_XM) - temp33_real = math_mul3333xx33(S_scale,P(i,j,k,1:3,1:3)) + math_I3 - F_lambda(i,j,k,1:3,1:3) & - + F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3) - f_scal(1:9,i,j,k) = reshape(temp33_real,[9]) - f_scal(10:18,i,j,k) = reshape(F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3),[9]) - enddo; enddo; enddo - - return end subroutine AL_formResidual !-------------------------------------------------------------------------------------------------- diff --git a/code/DAMASK_spectral_SolverBasic.f90 b/code/DAMASK_spectral_SolverBasic.f90 index 6967a327e..15eefd1f6 100644 --- a/code/DAMASK_spectral_SolverBasic.f90 +++ b/code/DAMASK_spectral_SolverBasic.f90 @@ -79,18 +79,18 @@ subroutine basic_init() write(6,'(a)') '' - allocate (F ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (P ( res(1), res(2),res(3),3,3), source = 0.0_pReal) + allocate (F ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) + allocate (F_lastInc ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) + allocate (P ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) !-------------------------------------------------------------------------------------------------- ! init fields if (restartInc == 1_pInt) then ! no deformation (no restart) + F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity + F_lastInc = F do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - F(i,j,k,1:3,1:3) = math_I3 - F_lastInc(i,j,k,1:3,1:3) = math_I3 coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - geomdim/real(2_pInt*res,pReal) enddo; enddo; enddo @@ -114,7 +114,7 @@ subroutine basic_init() coordinates = 0.0 ! change it later!!! endif - + !no rotation bc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) call Utilities_constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& P,C,temp33_Real,.false.,math_I3) @@ -259,17 +259,14 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme field_real = 0.0_pReal - field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = P + field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(P,[res(1),res(2),res(3),3,3],& + order=[4,5,1,2,3]) ! field real has a different order call Utilities_forwardFFT() err_div = Utilities_divergenceRMS() call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab) call Utilities_backwardFFT() - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization - enddo; enddo; enddo - + F = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av) - if (basic_solution%converged .and. iter > itmin) exit enddo convergenceLoop diff --git a/code/DAMASK_spectral_Utilities.f90 b/code/DAMASK_spectral_Utilities.f90 index 9d58fd8ac..404a7cbeb 100644 --- a/code/DAMASK_spectral_Utilities.f90 +++ b/code/DAMASK_spectral_Utilities.f90 @@ -71,6 +71,11 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, sets debug flags, create plans for fftw +!> @details Sets the debug levels for general, divergence, restart and fftw from the biwise coding +!> provided by the debug module to logicals. +!> Allocates all fields used by FFTW and create the corresponding plans depending on the debug +!> level chosen. +!> Initializes FFTW. !-------------------------------------------------------------------------------------------------- subroutine Utilities_init() @@ -185,15 +190,19 @@ subroutine Utilities_init() enddo; enddo; enddo if(memory_efficient) then ! allocate just single fourth order tensor - allocate (gamma_hat(1,1,1,3,3,3,3), source = 0.0_pReal) + allocate (gamma_hat(3,3,3,3,1,1,1), source = 0.0_pReal) else ! precalculation of gamma_hat field - allocate (gamma_hat(res1_red ,res(2),res(3),3,3,3,3), source =0.0_pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + allocate (gamma_hat(3,3,3,3,res1_red ,res(2),res(3)), source =0.0_pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 endif end subroutine Utilities_init !-------------------------------------------------------------------------------------------------- !> @brief updates references stiffness and potentially precalculated gamma operator +!> @details Sets the current reference stiffness to the stiffness given as an argument. +!> If the gamma operator is precalculated, it is calculated with this stiffness. +!> In case of a on-the-fly calculation, only the reference stiffness is updated. +!> The gamma operator is filtered depening on the filter selected in numerics !-------------------------------------------------------------------------------------------------- subroutine Utilities_updateGamma(C) @@ -218,15 +227,18 @@ subroutine Utilities_updateGamma(C) temp33_Real = math_inv33(temp33_Real) filter = Utilities_getFilter(xi(1:3,i,j,k)) forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)& - gamma_hat(i,j,k, l,m,n,o) = filter*temp33_Real(l,n)*xiDyad(m,o) + gamma_hat(l,m,n,o, i,j,k) = filter*temp33_Real(l,n)*xiDyad(m,o) endif enddo; enddo; enddo - gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + gamma_hat(1:3,1:3,1:3,1:3, 1,1,1) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 endif end subroutine Utilities_updateGamma !-------------------------------------------------------------------------------------------------- !> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed +!> Does an unweighted FFT transform from real to complex. +!> In case of debugging the FFT, also one component of the tensor (specified by row and column) +!> is independetly transformed complex to complex and compared to the whole tensor transform !-------------------------------------------------------------------------------------------------- subroutine Utilities_forwardFFT(row,column) use mesh, only : & @@ -276,6 +288,14 @@ subroutine Utilities_forwardFFT(row,column) = cmplx(0.0_pReal,0.0_pReal,pReal) end subroutine Utilities_forwardFFT + +!-------------------------------------------------------------------------------------------------- +!> @brief backward FFT of data in field_fourier to field_real +!> Does an inverse FFT transform from complex to real +!> In case of debugging the FFT, also one component of the tensor (specified by row and column) +!> is independetly transformed complex to complex and compared to the whole tensor transform +!> results is weighted by number of points stored in wgt +!-------------------------------------------------------------------------------------------------- subroutine Utilities_backwardFFT(row,column) implicit none @@ -286,9 +306,7 @@ subroutine Utilities_backwardFFT(row,column) !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 inverse FT results if (debugFFTW) then - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red - scalarField_fourier(i,j,k) = field_fourier(i,j,k,row,column) - enddo; enddo; enddo + scalarField_fourier = field_fourier(1:res1_red,1:res(2),1:res(3),row,column) do i = 0_pInt, res(1)/2_pInt-2_pInt ! unpack fft data for conj complex symmetric part m = 1_pInt do k = 1_pInt, res(3) @@ -303,8 +321,7 @@ subroutine Utilities_backwardFFT(row,column) enddo; enddo endif -!-------------------------------------------------------------------------------------------------- -! doing the inverse FT + call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient !-------------------------------------------------------------------------------------------------- @@ -353,16 +370,16 @@ subroutine Utilities_fourierConvolution(fieldAim) temp33_Real = math_inv33(temp33_Real) filter = Utilities_getFilter(xi(1:3,i,j,k)) forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)& - gamma_hat(1,1,1, l,m,n,o) = filter*temp33_Real(l,n)*xiDyad(m,o) + gamma_hat(l,m,n,o, 1,1,1) = filter*temp33_Real(l,n)*xiDyad(m,o) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - temp33_Complex(l,m) = sum(gamma_hat(1,1,1, l,m, 1:3,1:3) * field_fourier(i,j,k,1:3,1:3)) + temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3, 1,1,1) * field_fourier(i,j,k,1:3,1:3)) field_fourier(i,j,k,1:3,1:3) = temp33_Complex endif enddo; enddo; enddo else ! use precalculated gamma-operator do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt,res1_red forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt) & - temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) * field_fourier(i,j,k,1:3,1:3)) + temp33_Complex(m,n) = sum(gamma_hat(m,n,1:3,1:3, i,j,k) * field_fourier(i,j,k,1:3,1:3)) field_fourier(i,j,k, 1:3,1:3) = temp33_Complex enddo; enddo; enddo endif @@ -422,18 +439,9 @@ real(pReal) function Utilities_divergenceRMS() call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real) ! already weighted - err_real_div_RMS = 0.0_pReal - err_post_div_RMS = 0.0_pReal - err_real_div_max = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - err_real_div_RMS = err_real_div_RMS + sum(divergence_real(i,j,k,1:3)**2.0_pReal) ! avg of squared L_2 norm of div(stress) in real space - err_post_div_RMS = err_post_div_RMS + sum(divergence_post(i,j,k,1:3)**2.0_pReal) ! avg of squared L_2 norm of div(stress) in real space - err_real_div_max = max(err_real_div_max,sum(divergence_real(i,j,k,1:3)**2.0_pReal)) ! max of squared L_2 norm of div(stress) in real space - enddo; enddo; enddo - - err_real_div_RMS = sqrt(wgt*err_real_div_RMS) ! RMS in real space - err_post_div_RMS = sqrt(wgt*err_post_div_RMS) ! RMS in real space - err_real_div_max = sqrt( err_real_div_max) ! max in real space + err_real_div_RMS = sqrt(wgt*sum(divergence_real**2.0_pReal)) ! RMS in real space + err_post_div_RMS = sqrt(wgt*sum(divergence_post**2.0_pReal)) ! RMS in real space + err_real_div_max = sqrt(maxval(sum(divergence_real**2.0_pReal,dim=4))) ! max in real space err_div_max = sqrt( err_div_max) ! max in Fourier space write(6,'(1x,a,es11.4)') 'error divergence FT RMS = ',err_div_RMS @@ -516,7 +524,7 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti real(pReal), dimension(res(1),res(2),res(3)) :: temperature real(pReal), dimension(res(1),res(2),res(3),3) :: coordinates - real(pReal), dimension(res(1),res(2),res(3),3,3) :: F,F_lastInc, P + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: F,F_lastInc, P real(pReal) :: timeinc logical :: ForwardData integer(pInt) :: i, j, k, ielem @@ -524,7 +532,7 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti real(pReal), dimension(3,3,3,3) :: dPdF, C real(pReal), dimension(6) :: sigma ! cauchy stress real(pReal), dimension(6,6) :: dsde - real(pReal), dimension(3,3) :: P_av_lab, P_av, rotation_BC + real(pReal), dimension(3,3) :: P_av, rotation_BC write(6,'(a)') '' write(6,'(a)') '... evaluating constitutive response .................' @@ -539,28 +547,27 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt call CPFEM_general(3_pInt,& ! collect cycle - coordinates(i,j,k,1:3), F_lastInc(i,j,k,1:3,1:3),F(i,j,k,1:3,1:3), & - temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(i,j,k,1:3,1:3),dPdF) + coordinates(i,j,k,1:3), F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), & + temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) enddo; enddo; enddo P = 0.0_pReal ! needed because of the padding for FFTW C = 0.0_pReal - P_av_lab = 0.0_pReal ielem = 0_pInt call debug_reset() do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, - coordinates(i,j,k,1:3),F_lastInc(i,j,k,1:3,1:3), F(i,j,k,1:3,1:3), & ! others get 2 (saves winding forward effort) - temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(i,j,k,1:3,1:3),dPdF) + coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort) + temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) CPFEM_mode = 2_pInt C = C + dPdF - P_av_lab = P_av_lab + P(i,j,k,1:3,1:3) enddo; enddo; enddo call debug_info() + + P_av = math_rotate_forward33(sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt,rotation_BC) !average of P rotated restartWrite = .false. - P_av_lab = P_av_lab * wgt - P_av = math_rotate_forward33(P_av_lab,rotation_BC) + write (6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') ' Piola-Kirchhoff stress / MPa =',& math_transpose33(P_av)/1.e6_pReal @@ -574,18 +581,17 @@ subroutine Utilities_forwardField(delta_aim,timeinc,timeinc_old,guessmode,field_ real(pReal), intent(in), dimension(3,3) :: delta_aim real(pReal), intent(in) :: timeinc, timeinc_old, guessmode - real(pReal), intent(inout), dimension(res(1),res(2),res(3),3,3) :: field_lastInc,field - integer(pInt) :: i,j,k - real(pReal), dimension(3,3) :: temp33_real + real(pReal), intent(inout), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,field - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - temp33_Real = Field(i,j,k,1:3,1:3) - Field(i,j,k,1:3,1:3) = Field(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon - + guessmode * (field(i,j,k,1:3,1:3) - Field_lastInc(i,j,k,1:3,1:3))*timeinc/timeinc_old& ! guessing... - + (1.0_pReal-guessmode) * delta_aim ! if not guessing, use prescribed average deformation where applicable - Field_lastInc(i,j,k,1:3,1:3) = temp33_Real - enddo; enddo; enddo - end subroutine Utilities_forwardField + if (guessmode == 1.0_pReal) then + field = field + (field-field_lastInc) * timeinc/timeinc_old + field_lastInc = (field + field_lastInc * timeinc/timeinc_old) /(1.0_pReal + timeinc/timeinc_old) + else + field_lastInc = field + field = field + spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3)) + endif + +end subroutine Utilities_forwardField real(pReal) function Utilities_getFilter(k) @@ -614,9 +620,7 @@ subroutine Utilities_destroy() implicit none - if (debugDivergence) then - call fftw_destroy_plan(plan_divergence) - endif + if (debugDivergence) call fftw_destroy_plan(plan_divergence) if (debugFFTW) then call fftw_destroy_plan(plan_scalarField_forth) diff --git a/code/DAMASK_spectral_interface.f90 b/code/DAMASK_spectral_interface.f90 index 7e09892db..0f8911dde 100644 --- a/code/DAMASK_spectral_interface.f90 +++ b/code/DAMASK_spectral_interface.f90 @@ -17,7 +17,7 @@ ! along with DAMASK. If not, see . ! !-------------------------------------------------------------------------------------------------- -!* $Id$ +! $Id$ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH @@ -84,8 +84,8 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn) gotGeometry = .false. write(6,'(a)') '' - write(6,'(a)') '<<<+- DAMASK_spectral_interface init -+>>>' - write(6,'(a)') '$Id$' + write(6,'(a)') ' <<<+- DAMASK_spectral_interface init -+>>>' + write(6,'(a)') ' $Id$' #include "compilation_info.f90" if ( present(loadcaseParameterIn) .and. present(geometryParameterIn)) then ! both mandatory parameters given in function call