diff --git a/code/DAMASK_spectral_Driver.f90 b/code/DAMASK_spectral_Driver.f90 index f49e3a435..235e74ece 100644 --- a/code/DAMASK_spectral_Driver.f90 +++ b/code/DAMASK_spectral_Driver.f90 @@ -78,7 +78,8 @@ program DAMASK_spectral_Driver restartInc use numerics, only: & - rotation_tol + rotation_tol, & + myspectralsolver use homogenization, only: & materialpoint_sizeResults, & @@ -109,10 +110,8 @@ program DAMASK_spectral_Driver !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file - real(pReal), dimension(9) :: & - temp_valueVector !> temporarily from loadcase file when reading in tensors - logical, dimension(9) :: & - temp_maskVector !> temporarily from loadcase file when reading in tensors + real(pReal), dimension(9) :: temp_valueVector !> temporarily from loadcase file when reading in tensors + logical, dimension(9) :: temp_maskVector !> temporarily from loadcase file when reading in tensors integer(pInt), parameter :: maxNchunksLoadcase = (1_pInt + 9_pInt)*3_pInt +& ! deformation, rotation, and stress (1_pInt + 1_pInt)*5_pInt +& ! time, (log)incs, temp, restartfrequency, and outputfrequency 1_pInt, & ! dropguessing @@ -147,7 +146,7 @@ program DAMASK_spectral_Driver call DAMASK_interface_init write(6,'(a)') '' - write(6,'(a)') ' <<<+- DAMASK_spectral init -+>>>' + write(6,'(a)') ' <<<+- DAMASK_spectral_Driver init -+>>>' write(6,'(a)') ' $Id$' #include "compilation_info.f90" write(6,'(a)') ' Working Directory: ',trim(getSolverWorkingDirectoryName()) @@ -343,7 +342,15 @@ print*, 'my Unit closed' if (debugGeneral) write(6,'(a)') 'Header of result file written out' endif - call Basic_init() + select case (myspectralsolver) + + case (DAMASK_spectral_SolverBasic_label) + call basic_init() + + case (DAMASK_spectral_SolverAL_label) + call AL_init() + + end select !################################################################################################## ! Loop over loadcases defined in the currentLoadcase file @@ -390,7 +397,10 @@ print*, 'my Unit closed' write(6,'(a)') '##################################################################' write(6,'(A,I5.5,A,es12.5)') 'Increment ', totalIncsCounter, ' Time ',time - solres =basic_solution (& + select case (myspectralsolver) + + case (DAMASK_spectral_SolverBasic_label) + solres = basic_solution (& guessmode,timeinc,timeinc_old, & P_BC = bc(currentLoadcase)%stress, & F_BC = bc(currentLoadcase)%deformation, & @@ -398,6 +408,18 @@ print*, 'my Unit closed' mask_stressVector = bc(currentLoadcase)%maskStressVector, & velgrad = bc(currentLoadcase)%velGradApplied, & rotation_BC = bc(currentLoadcase)%rotation) + + case (DAMASK_spectral_SolverAL_label) + solres = AL_solution (& + guessmode,timeinc,timeinc_old, & + P_BC = bc(currentLoadcase)%stress, & + F_BC = bc(currentLoadcase)%deformation, & + ! temperature_bc = bc(currentLoadcase)%temperature, & + mask_stressVector = bc(currentLoadcase)%maskStressVector, & + velgrad = bc(currentLoadcase)%velGradApplied, & + rotation_BC = bc(currentLoadcase)%rotation) + + end select write(6,'(a)') '' write(6,'(a)') '==================================================================' @@ -420,6 +442,15 @@ print*, 'my Unit closed' enddo incLooping enddo loadCaseLooping + select case (myspectralsolver) + + case (DAMASK_spectral_SolverBasic_label) + call basic_destroy() + + case (DAMASK_spectral_SolverAL_label) + call AL_destroy() + + end select write(6,'(a)') '' write(6,'(a)') '##################################################################' write(6,'(i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & diff --git a/code/DAMASK_spectral_SolverAL.f90 b/code/DAMASK_spectral_SolverAL.f90 index ca3e947be..03923f7df 100644 --- a/code/DAMASK_spectral_SolverAL.f90 +++ b/code/DAMASK_spectral_SolverAL.f90 @@ -1,6 +1,263 @@ module DAMASK_spectral_SolverAL - use DAMASK_spectral_Utilities + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + + use DAMASK_spectral_Utilities + + use math + + use mesh, only : & + mesh_spectral_getResolution, & + mesh_spectral_getDimension + + implicit none + + character (len=*), parameter, public :: & + DAMASK_spectral_SolverAL_label = 'AL' +!-------------------------------------------------------------------------------------------------- +! common pointwise data + real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc, F_lambda, F_lambda_lastInc, P + real(pReal), dimension(:,:,:,:), allocatable :: coordinates + real(pReal), dimension(:,:,:), allocatable :: temperature + +!-------------------------------------------------------------------------------------------------- +! stress, stiffness and compliance average etc. + real(pReal), dimension(3,3) :: & + F_aim = math_I3, & + F_aim_lastInc = math_I3, & + P_av + real(pReal), dimension(3,3,3,3) :: & + C_ref = 0.0_pReal, & + C = 0.0_pReal + +!-------------------------------------------------------------------------------------------------- +! solution state + + contains + + subroutine AL_init() + + use IO, only: & + IO_read_JobBinaryFile, & + IO_write_JobBinaryFile + + use FEsolving, only: & + restartInc -end module DAMASK_spectral_SolverAL \ No newline at end of file + use DAMASK_interface, only: & + getSolverJobName + + implicit none + + integer(pInt) :: i, j, k + res = mesh_spectral_getResolution() + geomdim = mesh_spectral_getDimension() + + 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 (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) + 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 + 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 + 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,'F_aim',trim(getSolverJobName()),size(F_aim)) + read (777,rec=1) F_aim + close (777) + call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) + read (777,rec=1) F_aim_lastInc + close (777) + + coordinates = 0.0 ! change it later!!! + endif + + call constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& + P,C,P_av,.false.,math_I3) + +!-------------------------------------------------------------------------------------------------- +! reference stiffness + if (restartInc == 1_pInt) then + C_ref = C + call IO_write_jobBinaryFile(777,'C_ref',size(C_ref)) + write (777,rec=1) C_ref + close(777) + elseif (restartInc > 1_pInt) then + call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C_ref)) + read (777,rec=1) C_ref + close (777) + endif + + call Utilities_Init(C_ref) + + end subroutine AL_init + +type(solutionState) function AL_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,mask_stressVector,velgrad,rotation_BC) + + use numerics, only: & + itmax,& + itmin + + use IO, only: & + IO_write_JobBinaryFile + + use FEsolving, only: & + restartWrite + + implicit none + +!-------------------------------------------------------------------------------------------------- +! input data for solution + + real(pReal), intent(in) :: timeinc, timeinc_old + real(pReal), intent(in) :: guessmode + logical, intent(in) :: velgrad + real(pReal), dimension(3,3), intent(in) :: P_BC,F_BC,rotation_BC + logical, dimension(9), intent(in) :: mask_stressVector + +!-------------------------------------------------------------------------------------------------- +! loop variables, convergence etc. + + real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal + real(pReal), dimension(3,3) :: temp33_Real + real(pReal), dimension(3,3,3,3) :: S + real(pReal), dimension(3,3) :: mask_stress, & + mask_defgrad, & + deltaF_aim, & + F_aim_lab, & + F_aim_lab_lastIter + real(pReal) :: err_div, err_stress + integer(pInt) :: iter + integer(pInt) :: i, j, k + logical :: ForwardData + real(pReal) :: defgradDet + real(pReal) :: defgradDetMax, defgradDetMin + + mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3])) + mask_defgrad = merge(zeroes,ones,reshape(mask_stressVector,[3,3])) + + if (restartWrite) then + write(6,'(a)') 'writing converged results for restart' + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) ! writing deformation gradient field to file + write (777,rec=1) F_LastInc + close (777) + call IO_write_jobBinaryFile(777,'C',size(C)) + write (777,rec=1) C + close(777) + endif + + ForwardData = .True. + if (velgrad) then ! calculate deltaF_aim from given L and current F + deltaF_aim = timeinc * mask_defgrad * math_mul33x33(F_BC, F_aim) + else ! deltaF_aim = fDot *timeinc where applicable + deltaF_aim = timeinc * mask_defgrad * F_BC + endif + +!-------------------------------------------------------------------------------------------------- +! winding forward of deformation aim in loadcase system + temp33_Real = F_aim + F_aim = F_aim & + + guessmode * mask_stress * (F_aim - F_aim_lastInc)*timeinc/timeinc_old & + + 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) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + temp33_Real = F(i,j,k,1:3,1:3) + F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon + + guessmode * (F(i,j,k,1:3,1:3) - F_lastInc(i,j,k,1:3,1:3))*timeinc/timeinc_old& ! guessing... + + (1.0_pReal-guessmode) * deltaF_aim ! if not guessing, use prescribed average deformation where applicable + F_lastInc(i,j,k,1:3,1:3) = temp33_Real + enddo; enddo; enddo + call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),& ! calculate current coordinates + 1.0_pReal,F_lastInc,coordinates) + + iter = 0_pInt + S = S_lastInc(rotation_BC,mask_stressVector,C) + + convergenceLoop: do while((iter < itmax .and. (any([err_div ,err_stress] > 1.0_pReal)))& + .or. iter < itmin) + + iter = iter + 1_pInt +!-------------------------------------------------------------------------------------------------- +! report begin of new iteration + write(6,'(a)') '' + write(6,'(a)') '==================================================================' + write(6,'(3(a,i6.6))') ' @ Iter. ',itmin,' < ',iter,' < ',itmax + write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& + math_transpose33(F_aim) + F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC) + +!-------------------------------------------------------------------------------------------------- +! evaluate constitutive response + call constitutiveResponse(coordinates,F,F_lastInc,temperature,timeinc,& + P,C,P_av,ForwardData,rotation_BC) + ForwardData = .False. + +!-------------------------------------------------------------------------------------------------- +! stress BC handling + if(any(mask_stressVector)) then ! calculate stress BC if applied + err_stress = BCcorrection(mask_stressVector,P_BC,P_av,F_aim,S) + else + err_stress = 0.0_pReal + endif + + F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame + +!-------------------------------------------------------------------------------------------------- +! updated deformation gradient + field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = P + err_div = convolution(.True.,F_aim_lab_lastIter - F_aim_lab, C_ref) + + 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 + +!-------------------------------------------------------------------------------------------------- +! calculate bounds of det(F) and report + if(debugGeneral) then + defgradDetMax = -huge(1.0_pReal) + defgradDetMin = +huge(1.0_pReal) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + defgradDet = math_det33(F(i,j,k,1:3,1:3)) + defgradDetMax = max(defgradDetMax,defgradDet) + defgradDetMin = min(defgradDetMin,defgradDet) + enddo; enddo; enddo + + write(6,'(a,1x,es11.4)') 'max determinant of deformation =', defgradDetMax + write(6,'(a,1x,es11.4)') 'min determinant of deformation =', defgradDetMin + endif + enddo convergenceLoop + +end function AL_solution + +subroutine AL_destroy() + +implicit none + +call Utilities_destroy() + +end subroutine AL_destroy + +end module DAMASK_spectral_SolverAL diff --git a/code/DAMASK_spectral_SolverBasic.f90 b/code/DAMASK_spectral_SolverBasic.f90 index 82e630681..51fb14a8d 100644 --- a/code/DAMASK_spectral_SolverBasic.f90 +++ b/code/DAMASK_spectral_SolverBasic.f90 @@ -1,11 +1,112 @@ module DAMASK_spectral_SolverBasic -use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + use DAMASK_spectral_Utilities + + use math + + use mesh, only : & + mesh_spectral_getResolution, & + mesh_spectral_getDimension + + implicit none + + character (len=*), parameter, public :: & + DAMASK_spectral_SolverBasic_label = 'basic' + +!-------------------------------------------------------------------------------------------------- +! common pointwise data + real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc, P + real(pReal), dimension(:,:,:,:), allocatable :: coordinates + real(pReal), dimension(:,:,:), allocatable :: temperature + +!-------------------------------------------------------------------------------------------------- +! stress, stiffness and compliance average etc. + real(pReal), dimension(3,3) :: & + F_aim = math_I3, & + F_aim_lastInc = math_I3, & + P_av + real(pReal), dimension(3,3,3,3) :: & + C_ref = 0.0_pReal, & + C = 0.0_pReal + + contains -subroutine Basic_Init() - call Utilities_Init() -end subroutine basic_Init + subroutine basic_init() + + use IO, only: & + IO_read_JobBinaryFile, & + IO_write_JobBinaryFile + + use FEsolving, only: & + restartInc + + use DAMASK_interface, only: & + getSolverJobName + + implicit none + integer(pInt) :: i,j,k + + res = mesh_spectral_getResolution() + geomdim = mesh_spectral_getDimension() + + 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 (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) + 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 + 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 + 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,'F_aim',trim(getSolverJobName()),size(F_aim)) + read (777,rec=1) F_aim + close (777) + call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) + read (777,rec=1) F_aim_lastInc + close (777) + + coordinates = 0.0 ! change it later!!! + endif + + call constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& + P,C,P_av,.false.,math_I3) + +!-------------------------------------------------------------------------------------------------- +! reference stiffness + if (restartInc == 1_pInt) then + C_ref = C + call IO_write_jobBinaryFile(777,'C_ref',size(C_ref)) + write (777,rec=1) C_ref + close(777) + elseif (restartInc > 1_pInt) then + call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C_ref)) + read (777,rec=1) C_ref + close (777) + endif + + call Utilities_Init(C_ref) + + end subroutine basic_init type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,mask_stressVector,velgrad,rotation_BC) @@ -19,45 +120,34 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F use FEsolving, only: & restartWrite - - implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution + real(pReal), intent(in) :: timeinc, timeinc_old real(pReal), intent(in) :: guessmode - logical, intent(in) :: velgrad + logical, intent(in) :: velgrad real(pReal), dimension(3,3), intent(in) :: P_BC,F_BC,rotation_BC - logical, dimension(9), intent(in) :: mask_stressVector - - - real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal - - - - - - real(pReal), dimension(3,3) :: temp33_Real ! compliance and stiffness in matrix notation - - - -real(pReal), dimension(3,3,3,3) :: S -real(pReal), dimension(3,3) :: & - mask_stress, & - mask_defgrad, & - deltaF_aim, & - F_aim_lab, & - F_aim_lab_lastIter + logical, dimension(9), intent(in) :: mask_stressVector !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. - real(pReal) :: err_div, err_stress + + real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal + real(pReal), dimension(3,3) :: temp33_Real + real(pReal), dimension(3,3,3,3) :: S + real(pReal), dimension(3,3) :: mask_stress, & + mask_defgrad, & + deltaF_aim, & + F_aim_lab, & + F_aim_lab_lastIter + real(pReal) :: err_div, err_stress integer(pInt) :: iter integer(pInt) :: i, j, k - logical :: ForwardResults - real(pReal) :: defgradDet - real(pReal) :: defgradDetMax, defgradDetMin + logical :: ForwardData + real(pReal) :: defgradDet + real(pReal) :: defgradDetMax, defgradDetMin mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3])) mask_defgrad = merge(zeroes,ones,reshape(mask_stressVector,[3,3])) @@ -72,7 +162,7 @@ real(pReal), dimension(3,3) :: & close(777) endif - ForwardResults = .True. + ForwardData = .True. if (velgrad) then ! calculate deltaF_aim from given L and current F deltaF_aim = timeinc * mask_defgrad * math_mul33x33(F_BC, F_aim) else ! deltaF_aim = fDot *timeinc where applicable @@ -102,7 +192,8 @@ real(pReal), dimension(3,3) :: & 1.0_pReal,F_lastInc,coordinates) iter = 0_pInt -S = S_lastInc(rotation_BC,mask_stressVector) + S = S_lastInc(rotation_BC,mask_stressVector,C) + convergenceLoop: do while((iter < itmax .and. (any([err_div ,err_stress] > 1.0_pReal)))& .or. iter < itmin) @@ -118,23 +209,24 @@ S = S_lastInc(rotation_BC,mask_stressVector) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response - - call constitutiveResponse(ForwardResults,timeInc) - ForwardResults = .False. + call constitutiveResponse(coordinates,F,F_lastInc,temperature,timeinc,& + P,C,P_av,ForwardData,rotation_BC) + ForwardData = .False. + !-------------------------------------------------------------------------------------------------- ! stress BC handling - if(any(mask_stressVector)) then ! calculate stress BC if applied - err_stress = BCcorrection(mask_stressVector,P_BC,S) - else - err_stress = 0.0_pReal - endif + if(any(mask_stressVector)) then ! calculate stress BC if applied + err_stress = BCcorrection(mask_stressVector,P_BC,P_av,F_aim,S) + else + err_stress = 0.0_pReal + endif F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame !-------------------------------------------------------------------------------------------------- ! updated deformation gradient field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = P - err_div = convolution(.True.,F_aim_lab_lastIter - F_aim_lab) + err_div = convolution(.True.,F_aim_lab_lastIter - F_aim_lab, C_ref) 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 @@ -158,4 +250,12 @@ S = S_lastInc(rotation_BC,mask_stressVector) end function basic_solution -end module DAMASK_spectral_SolverBasic \ No newline at end of file +subroutine basic_destroy() + +implicit none + +call Utilities_destroy() + +end subroutine basic_destroy + +end module DAMASK_spectral_SolverBasic diff --git a/code/DAMASK_spectral_Utilities.f90 b/code/DAMASK_spectral_Utilities.f90 index 6e900122c..735498286 100644 --- a/code/DAMASK_spectral_Utilities.f90 +++ b/code/DAMASK_spectral_Utilities.f90 @@ -46,20 +46,6 @@ module DAMASK_spectral_Utilities IO_error implicit none - - type solutionState ! mask of stress boundary conditions - logical :: converged = .false. - logical :: regrid = .false. - logical :: term_ill = .false. - end type solutionState - - character(len=5) :: solverType, parameter = 'basic' -!-------------------------------------------------------------------------------------------------- -! common pointwise data - real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc, P - real(pReal), dimension(:,:,:,:), allocatable :: coordinates - real(pReal), dimension(:,:,:), allocatable :: temperature - !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW @@ -86,28 +72,22 @@ module DAMASK_spectral_Utilities !variables controlling debugging logical :: debugGeneral, debugDivergence, debugRestart, debugFFTW -!-------------------------------------------------------------------------------------------------- -! stress, stiffness and compliance average etc. - real(pReal), dimension(3,3) :: & - F_aim = math_I3, & - F_aim_lastInc = math_I3, & - P_av - - real(pReal), dimension(3,3,3,3) :: & - C_ref = 0.0_pReal, & - C = 0.0_pReal real(pReal), dimension(3) :: geomdim = 0.0_pReal, virt_dim = 0.0_pReal ! physical dimension of volume element per direction integer(pInt), dimension(3) :: res = 1_pInt real(pReal) :: wgt integer(pInt) :: res1_red, Npoints +!-------------------------------------------------------------------------------------------------- +! solution state + type solutionState + logical :: converged = .false. + logical :: regrid = .false. + logical :: term_ill = .false. + end type solutionState contains -subroutine Utilities_init(F,P,F_...) - - use DAMASK_interface, only: & - getSolverJobName +subroutine Utilities_init(C_ref) use mesh, only : & mesh_spectral_getResolution, & @@ -126,24 +106,17 @@ subroutine Utilities_init(F,P,F_...) debug_spectralDivergence, & debug_spectralRestart, & debug_spectralFFTW - - use FEsolving, only: & - restartInc + use numerics, only: & memory_efficient - use CPFEM, only: & - CPFEM_general - - use IO, only: & - IO_read_JobBinaryFile, & - IO_write_JobBinaryFile - implicit none real(pReal), dimension(3,3) :: temp33_Real, xiDyad integer(pInt) :: i, j, k, l, m, n, q, ierr integer(pInt), dimension(3) :: k_s + real(pReal), dimension(3,3,3,3) :: & + C_ref type(C_PTR) :: tensorField ! field in real and fourier space type(C_PTR) :: scalarField_realC, scalarField_fourierC @@ -151,7 +124,7 @@ subroutine Utilities_init(F,P,F_...) write(6,'(a)') '' - write(6,'(a)') ' <<<+- DAMASK_spectralSolver init -+>>>' + write(6,'(a)') ' <<<+- DAMASK_spectralSolver Utilities init -+>>>' write(6,'(a)') ' $Id$' #include "compilation_info.f90" write(6,'(a)') '' @@ -172,12 +145,7 @@ subroutine Utilities_init(F,P,F_...) Npoints = res(1)*res(2)*res(3) wgt = 1.0/real(Npoints,pReal) - 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 (xi (3,res1_red,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) ! start out isothermally + allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) ! start out isothermally tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8) call c_f_pointer(tensorField, field_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField call c_f_pointer(tensorField, field_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField @@ -232,37 +200,6 @@ subroutine Utilities_init(F,P,F_...) endif if (debugGeneral) write(6,'(a)') 'FFTW initialized' - -!-------------------------------------------------------------------------------------------------- -! init fields - if (restartInc == 1_pInt) then ! no deformation (no restart) - 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 - 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 - 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,'F_aim',trim(getSolverJobName()),size(F_aim)) - read (777,rec=1) F_aim - close (777) - call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) - read (777,rec=1) F_aim_lastInc - close (777) - - coordinates = 0.0 ! change it later!!! - endif - call constitutiveResponse(.FALSE.,0.0_pReal) !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) @@ -284,19 +221,6 @@ subroutine Utilities_init(F,P,F_...) k_s(1) = i - 1_pInt xi(1:3,i,j,k) = real(k_s, pReal)/virt_dim enddo; enddo; enddo - -!-------------------------------------------------------------------------------------------------- -! calculate the gamma operator - if (restartInc == 1_pInt) then - C_ref = C - call IO_write_jobBinaryFile(777,'C_ref',size(C_ref)) - write (777,rec=1) C_ref - close(777) - elseif (restartInc > 1_pInt) then - call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C_ref)) - read (777,rec=1) C_ref - close (777) - endif if(memory_efficient) then ! allocate just single fourth order tensor allocate (gamma_hat(1,1,1,3,3,3,3), source = 0.0_pReal) @@ -317,252 +241,250 @@ subroutine Utilities_init(F,P,F_...) endif end subroutine Utilities_init -real(pReal) function convolution(calcDivergence, field_aim,) - use numerics, only: & - memory_efficient, & - err_div_tol - real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors - real(pReal) :: err_div = 0.0_pReal - real(pReal), dimension(3,3) :: temp33_Real - integer(pInt) :: i, j, k, l, m, n, q - +real(pReal) function convolution(calcDivergence,field_aim,C_ref) + + use numerics, only: & + memory_efficient, & + err_div_tol + + real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors + real(pReal) :: err_div = 0.0_pReal + real(pReal), dimension(3,3) :: temp33_Real + integer(pInt) :: i, j, k, l, m, n, q + real(pReal), dimension(3,3,3,3) :: C_ref + !-------------------------------------------------------------------------------------------------- !variables for additional output due to general debugging - real(pReal) :: maxCorrectionSym, maxCorrectionSkew - logical :: calcDivergence - real(pReal), dimension(3,3) :: field_avg, field_aim + real(pReal) :: maxCorrectionSym, maxCorrectionSkew + logical :: calcDivergence + real(pReal), dimension(3,3) :: field_avg, field_aim integer(pInt) :: row, column - real(pReal) :: field_av_L2, err_div_RMS, err_real_div_RMS, err_post_div_RMS,& - err_div_max, err_real_div_max - complex(pReal), dimension(3) :: temp3_complex - complex(pReal), dimension(3,3) :: temp33_complex + real(pReal) :: field_av_L2, err_div_RMS, err_real_div_RMS, err_post_div_RMS,& + err_div_max, err_real_div_max + complex(pReal), dimension(3) :: temp3_complex + complex(pReal), dimension(3,3) :: temp33_complex - !-------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- ! actual spectral method - write(6,'(a)') '' - write(6,'(a)') '... doing convolution .................' + write(6,'(a)') '' + write(6,'(a)') '... doing convolution .................' !-------------------------------------------------------------------------------------------------- ! copy one component of the stress field to to a single FT and check for mismatch - if (debugFFTW) then - row = 3 ! (mod(totalIncsCounter+iter-2_pInt,9_pInt))/3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1 - column = 3 !(mod(totalIncsCounter+iter-2_pInt,3_pInt)) + 1_pInt - scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component - cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) - endif + if (debugFFTW) then + row = 3 ! (mod(totalIncsCounter+iter-2_pInt,9_pInt))/3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1 + column = 3 !(mod(totalIncsCounter+iter-2_pInt,3_pInt)) + 1_pInt + scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component + cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) + endif !-------------------------------------------------------------------------------------------------- ! call function to calculate divergence from math (for post processing) to check results - if (debugDivergence) & - call divergence_fft(res,virt_dim,3_pInt,& - field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),divergence_post) ! padding - + if (debugDivergence) & + call divergence_fft(res,virt_dim,3_pInt,& + field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),divergence_post) ! padding + !-------------------------------------------------------------------------------------------------- ! doing the FT because it simplifies calculation of average stress in real space also - call fftw_execute_dft_r2c(plan_forward,field_real,field_fourier) + call fftw_execute_dft_r2c(plan_forward,field_real,field_fourier) !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 FT results - if (debugFFTW) then - call fftw_execute_dft(plan_scalarField_forth,scalarField_real,scalarField_fourier) - write(6,'(a,i1,1x,i1)') 'checking FT results of compontent ', row, column - write(6,'(a,2(es11.4,1x))') 'max FT relative error = ',& - maxval( real((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& - field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& - scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), & - maxval(aimag((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& - field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& - scalarField_fourier(1:res1_red,1:res(2),1:res(3)))) - endif + if (debugFFTW) then + call fftw_execute_dft(plan_scalarField_forth,scalarField_real,scalarField_fourier) + write(6,'(a,i1,1x,i1)') 'checking FT results of compontent ', row, column + write(6,'(a,2(es11.4,1x))') 'max FT relative error = ',& + maxval( real((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& + field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& + scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), & + maxval(aimag((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& + field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/& + scalarField_fourier(1:res1_red,1:res(2),1:res(3)))) + endif !-------------------------------------------------------------------------------------------------- ! removing highest frequencies - field_fourier ( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& - = cmplx(0.0_pReal,0.0_pReal,pReal) - field_fourier (1:res1_red, res(2)/2_pInt+1_pInt,1:res(3) ,1:3,1:3)& - = cmplx(0.0_pReal,0.0_pReal,pReal) - if(res(3)>1_pInt) & - field_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& - = cmplx(0.0_pReal,0.0_pReal,pReal) + field_fourier ( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& + = cmplx(0.0_pReal,0.0_pReal,pReal) + field_fourier (1:res1_red, res(2)/2_pInt+1_pInt,1:res(3) ,1:3,1:3)& + = cmplx(0.0_pReal,0.0_pReal,pReal) + if(res(3)>1_pInt) & + field_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& + = cmplx(0.0_pReal,0.0_pReal,pReal) !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space - if( calcDivergence) then - field_avg = real(field_fourier(1,1,1,1:3,1:3),pReal)*wgt - - field_av_L2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(field_avg,& ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html) - math_transpose33(field_avg))))) - err_div_RMS = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2) - do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. - err_div_RMS = err_div_RMS & - + 2.0_pReal*(sum (real(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again - xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector - +sum(aimag(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& - xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) - enddo - err_div_RMS = err_div_RMS & ! Those two layers (DC and Nyquist) do not have a conjugate complex counterpart - + sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& - xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& - + sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& - xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& - + sum( real(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& - xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)& - + sum(aimag(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& - xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) - enddo; enddo - - err_div_RMS = sqrt(err_div_RMS)*wgt ! RMS in real space calculated with Parsevals theorem from Fourier space - err_div = err_div_RMS/field_av_L2 ! criterion to stop iterations - - + if(calcDivergence) then + field_avg = real(field_fourier(1,1,1,1:3,1:3),pReal)*wgt + + field_av_L2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(field_avg,& ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html) + math_transpose33(field_avg))))) + err_div_RMS = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2) + do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. + err_div_RMS = err_div_RMS & + + 2.0_pReal*(sum (real(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& ! (sqrt(real(a)**2 + aimag(a)**2))**2 = real(a)**2 + aimag(a)**2. do not take square root and square again + xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)& ! --> sum squared L_2 norm of vector + +sum(aimag(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),& + xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) + enddo + err_div_RMS = err_div_RMS & ! Those two layers (DC and Nyquist) do not have a conjugate complex counterpart + + sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& + xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& + + sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& + xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& + + sum( real(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& + xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)& + + sum(aimag(math_mul33x3_complex(field_fourier(res1_red,j,k,1:3,1:3),& + xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal) + enddo; enddo + + err_div_RMS = sqrt(err_div_RMS)*wgt ! RMS in real space calculated with Parsevals theorem from Fourier space + err_div = err_div_RMS/field_av_L2 ! criterion to stop iterations + + !-------------------------------------------------------------------------------------------------- ! calculate additional divergence criteria and report - if (debugDivergence) then ! calculate divergence again - err_div_max = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red - temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier - xi(1:3,i,j,k))*TWOPIIMG - err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal)) - divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared - enddo; enddo; enddo - - 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_div_max = sqrt( err_div_max) ! max in Fourier space - - write(6,'(a,es11.4)') 'error divergence FT RMS = ',err_div_RMS - write(6,'(a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS - write(6,'(a,es11.4)') 'error divergence post RMS = ',err_post_div_RMS - write(6,'(a,es11.4)') 'error divergence FT max = ',err_div_max - write(6,'(a,es11.4)') 'error divergence Real max = ',err_real_div_max - endif - write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/err_div_tol,& - ' (',err_div,' N/m³)' - end if + if (debugDivergence) then ! calculate divergence again + err_div_max = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red + temp3_Complex = math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3)*wgt,& ! weighting P_fourier + xi(1:3,i,j,k))*TWOPIIMG + err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal)) + divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared + enddo; enddo; enddo + + 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_div_max = sqrt( err_div_max) ! max in Fourier space + + write(6,'(a,es11.4)') 'error divergence FT RMS = ',err_div_RMS + write(6,'(a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS + write(6,'(a,es11.4)') 'error divergence post RMS = ',err_post_div_RMS + write(6,'(a,es11.4)') 'error divergence FT max = ',err_div_max + write(6,'(a,es11.4)') 'error divergence Real max = ',err_real_div_max + endif + write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/err_div_tol,& + ' (',err_div,' N/m³)' + end if !-------------------------------------------------------------------------------------------------- ! to the actual spectral method calculation (mechanical equilibrium) - if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat - - do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red - if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 - forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) - forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - temp33_Real(l,m) = sum(C_ref(l,m,1:3,1:3)*xiDyad) - temp33_Real = math_inv33(temp33_Real) - forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, q=1_pInt:3_pInt)& - gamma_hat(1,1,1, l,m,n,q) = temp33_Real(l,n)*xiDyad(m,q) - 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)) - 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)) - field_fourier(i,j,k, 1:3,1:3) = temp33_Complex - enddo; enddo; enddo - - endif - field_fourier(1,1,1,1:3,1:3) = cmplx(field_aim,0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat + do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red + if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + temp33_Real(l,m) = sum(C_ref(l,m,1:3,1:3)*xiDyad) + temp33_Real = math_inv33(temp33_Real) + forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, q=1_pInt:3_pInt)& + gamma_hat(1,1,1, l,m,n,q) = temp33_Real(l,n)*xiDyad(m,q) + 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)) + 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)) + field_fourier(i,j,k, 1:3,1:3) = temp33_Complex + enddo; enddo; enddo + endif + field_fourier(1,1,1,1:3,1:3) = cmplx(field_aim,0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 !-------------------------------------------------------------------------------------------------- ! 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 - 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) - n = 1_pInt - do j = 1_pInt, res(2) - scalarField_fourier(res(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m)) - if(n == 1_pInt) n = res(2) + 1_pInt - n = n-1_pInt - enddo - if(m == 1_pInt) m = res(3) + 1_pInt - m = m -1_pInt - enddo; enddo - endif + 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 + 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) + n = 1_pInt + do j = 1_pInt, res(2) + scalarField_fourier(res(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m)) + if(n == 1_pInt) n = res(2) + 1_pInt + n = n-1_pInt + enddo + if(m == 1_pInt) m = res(3) + 1_pInt + m = m -1_pInt + enddo; enddo + endif !-------------------------------------------------------------------------------------------------- ! doing the inverse FT - call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient + call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 inverse FT results - if (debugFFTW) then - write(6,'(a,i1,1x,i1)') 'checking iFT results of compontent ', row, column - call fftw_execute_dft(plan_scalarField_back,scalarField_fourier,scalarField_real) - write(6,'(a,es11.4)') 'max iFT relative error = ',& - maxval((real(scalarField_real(1:res(1),1:res(2),1:res(3)))-& - field_real(1:res(1),1:res(2),1:res(3),row,column))/& - real(scalarField_real(1:res(1),1:res(2),1:res(3)))) - endif + if (debugFFTW) then + write(6,'(a,i1,1x,i1)') 'checking iFT results of compontent ', row, column + call fftw_execute_dft(plan_scalarField_back,scalarField_fourier,scalarField_real) + write(6,'(a,es11.4)') 'max iFT relative error = ',& + maxval((real(scalarField_real(1:res(1),1:res(2),1:res(3)))-& + field_real(1:res(1),1:res(2),1:res(3),row,column))/& + real(scalarField_real(1:res(1),1:res(2),1:res(3)))) + endif !-------------------------------------------------------------------------------------------------- ! calculate some additional output - if(debugGeneral) then - maxCorrectionSkew = 0.0_pReal - maxCorrectionSym = 0.0_pReal - temp33_Real = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - maxCorrectionSym = max(maxCorrectionSym,& - maxval(math_symmetric33(field_real(i,j,k,1:3,1:3)))) - maxCorrectionSkew = max(maxCorrectionSkew,& - maxval(math_skew33(field_real(i,j,k,1:3,1:3)))) - temp33_Real = temp33_Real + field_real(i,j,k,1:3,1:3) - enddo; enddo; enddo - write(6,'(a,1x,es11.4)') 'max symmetric correction of deformation =',& - maxCorrectionSym*wgt - write(6,'(a,1x,es11.4)') 'max skew correction of deformation =',& - maxCorrectionSkew*wgt - write(6,'(a,1x,es11.4)') 'max sym/skew of avg correction = ',& - maxval(math_symmetric33(temp33_real))/& - maxval(math_skew33(temp33_real)) - endif - field_real = field_real * wgt - convolution = err_div/err_div_tol + if(debugGeneral) then + maxCorrectionSkew = 0.0_pReal + maxCorrectionSym = 0.0_pReal + temp33_Real = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + maxCorrectionSym = max(maxCorrectionSym,& + maxval(math_symmetric33(field_real(i,j,k,1:3,1:3)))) + maxCorrectionSkew = max(maxCorrectionSkew,& + maxval(math_skew33(field_real(i,j,k,1:3,1:3)))) + temp33_Real = temp33_Real + field_real(i,j,k,1:3,1:3) + enddo; enddo; enddo + write(6,'(a,1x,es11.4)') 'max symmetric correction of deformation =',& + maxCorrectionSym*wgt + write(6,'(a,1x,es11.4)') 'max skew correction of deformation =',& + maxCorrectionSkew*wgt + write(6,'(a,1x,es11.4)') 'max sym/skew of avg correction = ',& + maxval(math_symmetric33(temp33_real))/& + maxval(math_skew33(temp33_real)) + endif + field_real = field_real * wgt + convolution = err_div/err_div_tol + end function convolution -function S_lastInc(rot_BC,mask_stressVector1) +function S_lastInc(rot_BC,mask_stressVector1,C) + real(pReal), dimension(3,3,3,3) :: S_lastInc - integer(pInt) :: i, j, k, m,n + real(pReal), dimension(3,3,3,3), intent(in) :: C + integer(pInt) :: i, j, k, m, n real(pReal), dimension(3,3), intent(in) :: rot_BC logical, dimension(9), intent(in) :: mask_stressVector1 real(pReal), dimension(3,3,3,3) :: C_lastInc - real(pReal), dimension(9,9) :: temp99_Real + real(pReal), dimension(9,9) :: temp99_Real integer(pInt) :: size_reduced = 0_pInt real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC) logical :: errmatinv size_reduced = count(mask_stressVector1) - if (allocated(c_reduced)) deallocate(c_reduced) - if (allocated(c_reduced)) deallocate(c_reduced) allocate (c_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal) - - C_lastInc = math_rotate_forward3333(C,rot_BC) ! calculate stiffness from former inc temp99_Real = math_Plain3333to99(C_lastInc) k = 0_pInt ! build reduced stiffness @@ -588,9 +510,7 @@ function S_lastInc(rot_BC,mask_stressVector1) j = j + 1_pInt temp99_Real(n,m) = s_reduced(k,j) endif; enddo; endif; enddo - S_lastInc = (math_Plain99to3333(temp99_Real)) - if (allocated(c_reduced)) deallocate(c_reduced) - if (allocated(c_reduced)) deallocate(c_reduced) + S_lastInc = math_Plain99to3333(temp99_Real) end function S_lastInc @@ -598,76 +518,90 @@ end function S_lastInc !-------------------------------------------------------------------------------------------------- ! calculate reduced compliance -real(pReal) function BCcorrection(mask_stressVector,P_BC,S_lastInc) -use numerics, only: err_stress_tolrel, err_stress_tolabs - logical, dimension(9) :: mask_stressVector - real(pReal) :: err_stress, err_stress_tol - real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal - real(pReal), dimension(3,3,3,3) :: S_lastInc - real(pReal), dimension(3,3) :: & - P_BC , & - mask_stress, & - mask_defgrad - mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3])) - mask_defgrad = merge(zeroes,ones,reshape(mask_stressVector,[3,3])) + real(pReal) function BCcorrection(mask_stressVector,P_BC,P_av,F_aim,S_lastInc) + + use numerics, only: err_stress_tolrel, err_stress_tolabs + + logical, dimension(9) :: mask_stressVector + real(pReal) :: err_stress, err_stress_tol + real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal + real(pReal), dimension(3,3,3,3) :: S_lastInc + real(pReal), dimension(3,3) :: & + P_BC , & + P_av, & + F_aim, & + mask_stress, & + mask_defgrad + mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3])) + mask_defgrad = merge(zeroes,ones,reshape(mask_stressVector,[3,3])) !-------------------------------------------------------------------------------------------------- ! stress BC handling - ! calculate stress BC if applied - err_stress = maxval(abs(mask_stress * (P_av - P_BC))) ! maximum deviaton (tensor norm not applicable) - err_stress_tol = min(maxval(abs(P_av)) * err_stress_tolrel,err_stress_tolabs) ! don't use any tensor norm for the relative criterion because the comparison should be coherent - write(6,'(a)') '' - write(6,'(a)') '... correcting deformation gradient to fulfill BCs ...............' - write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, & - ' (',err_stress,' Pa)' - F_aim = F_aim - math_mul3333xx33(S_lastInc, ((P_av - P_BC))) ! residual on given stress components - write(6,'(a,1x,es11.4)')'determinant of new deformation = ',math_det33(F_aim) - BCcorrection = err_stress/err_stress_tol +! calculate stress BC if applied + err_stress = maxval(abs(mask_stress * (P_av - P_BC))) ! maximum deviaton (tensor norm not applicable) + err_stress_tol = min(maxval(abs(P_av)) * err_stress_tolrel,err_stress_tolabs) ! don't use any tensor norm for the relative criterion because the comparison should be coherent + write(6,'(a)') '' + write(6,'(a)') '... correcting deformation gradient to fulfill BCs ...............' + write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, & + ' (',err_stress,' Pa)' + F_aim = F_aim - math_mul3333xx33(S_lastInc, ((P_av - P_BC))) ! residual on given stress components + write(6,'(a,1x,es11.4)')'determinant of new deformation = ',math_det33(F_aim) + BCcorrection = err_stress/err_stress_tol + end function BCcorrection -subroutine constitutiveResponse(F,P,ForwardData,timeinc) - use debug, only: & - debug_reset, & - debug_info - use CPFEM, only: & - CPFEM_general -use FEsolving, only: restartWrite +subroutine constitutiveResponse(coordinates,F,F_lastInc,temperature,timeinc,& + P,C,P_av,ForwardData,rotation_BC) + use debug, only: & + debug_reset, & + debug_info + use CPFEM, only: & + CPFEM_general + use FEsolving, only: restartWrite + + implicit none + + 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) :: timeinc logical :: ForwardData integer(pInt) :: i, j, k, ielem integer(pInt) :: CPFEM_mode - real(pReal), dimension(3,3,3,3) :: dPdF - real(pReal), dimension(6) :: sigma ! cauchy stress - real(pReal), dimension(6,6) :: dsde - real(pReal), dimension(3,3) :: P_av_lab, rotation_BC + 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 + if (ForwardData) then CPFEM_mode = 1_pInt else CPFEM_mode = 2_pInt endif - write(6,'(a)') '' - write(6,'(a)') '... update stress field P(F) .....................................' - ielem = 0_pInt - 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 + write(6,'(a)') '' + write(6,'(a)') '... update stress field P(F) .....................................' + ielem = 0_pInt + 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) - enddo; enddo; enddo + 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) - CPFEM_mode = 2_pInt - C = C + dPdF - P_av_lab = P_av_lab + P(i,j,k,1:3,1:3) + 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) + 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() restartWrite = .false. @@ -677,5 +611,23 @@ use FEsolving, only: restartWrite math_transpose33(P_av)/1.e6_pReal C = C * wgt end subroutine constitutiveResponse + +subroutine Utilities_destroy + +implicit none + +if (debugDivergence) then + call fftw_destroy_plan(plan_divergence) +endif + +if (debugFFTW) then + call fftw_destroy_plan(plan_scalarField_forth) + call fftw_destroy_plan(plan_scalarField_back) +endif + +call fftw_destroy_plan(plan_forward) +call fftw_destroy_plan(plan_backward) + +end subroutine Utilities_destroy end module DAMASK_spectral_Utilities diff --git a/code/numerics.f90 b/code/numerics.f90 index c4bac1339..dd849d07a 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -79,7 +79,8 @@ real(pReal) :: err_div_tol = 0.1_pReal, & err_stress_tolabs = huge(1.0_pReal), & ! absolute tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress fftw_timelimit = -1.0_pReal, & ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit rotation_tol = 1.0e-12_pReal ! tolerance of rotation specified in loadcase, Default 1.0e-12: first guess -character(len=64) :: fftw_plan_mode = 'FFTW_PATIENT' ! reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag +character(len=64) :: fftw_plan_mode = 'FFTW_PATIENT', & ! reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag + myspectralsolver = 'basic' ! spectral solution method integer(pInt) :: fftw_planner_flag = 32_pInt, & ! conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw itmax = 20_pInt, & ! maximum number of iterations itmin = 2_pInt ! minimum number of iterations @@ -247,6 +248,8 @@ subroutine numerics_init fftw_timelimit = IO_floatValue(line,positions,2_pInt) case ('fftw_plan_mode') fftw_plan_mode = IO_stringValue(line,positions,2_pInt) + case ('myspectralsolver') + myspectralsolver = IO_stringValue(line,positions,2_pInt) case ('rotation_tol') rotation_tol = IO_floatValue(line,positions,2_pInt) case ('divergence_correction') @@ -256,7 +259,7 @@ subroutine numerics_init #endif #ifndef Spectral case ('err_div_tol','err_stress_tolrel','err_stress_tolabs',& - 'itmax', 'itmin','memory_efficient','fftw_timelimit','fftw_plan_mode', & + 'itmax', 'itmin','memory_efficient','fftw_timelimit','fftw_plan_mode','myspectralsolver', & 'rotation_tol','divergence_correction','update_gamma') call IO_warning(40_pInt,ext_msg=tag) #endif @@ -348,6 +351,7 @@ subroutine numerics_init write(6,'(a24,1x,es8.1)') ' fftw_timelimit: ',fftw_timelimit endif write(6,'(a24,1x,a)') ' fftw_plan_mode: ',trim(fftw_plan_mode) + write(6,'(a24,1x,a)') ' myspectralsolver: ',trim(myspectralsolver) write(6,'(a24,1x,i8)') ' fftw_planner_flag: ',fftw_planner_flag write(6,'(a24,1x,es8.1)') ' rotation_tol: ',rotation_tol write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction