223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

next round in modularization

This commit is contained in:
Martin Diehl 2012-07-25 14:01:39 +00:00
parent 9c7a826f00
commit 4ed68bb4ae
5 changed files with 738 additions and 394 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
code/DAMASK_spectral_Driver.f90
View File

@ -78,7 +78,8 @@ program DAMASK_spectral_Driver
restartInc restartInc
use numerics, only: & use numerics, only: &
rotation_tol rotation_tol, &
myspectralsolver
use homogenization, only: & use homogenization, only: &
materialpoint_sizeResults, & materialpoint_sizeResults, &
@ -109,10 +110,8 @@ program DAMASK_spectral_Driver
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file ! variables related to information from load case and geom file
real(pReal), dimension(9) :: & real(pReal), dimension(9) :: temp_valueVector !> temporarily from loadcase file when reading in tensors
temp_valueVector !> temporarily from loadcase file when reading in tensors logical, dimension(9) :: temp_maskVector !> 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 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 + 1_pInt)*5_pInt +& ! time, (log)incs, temp, restartfrequency, and outputfrequency
1_pInt, & ! dropguessing 1_pInt, & ! dropguessing
@ -147,7 +146,7 @@ program DAMASK_spectral_Driver
call DAMASK_interface_init call DAMASK_interface_init
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') ' <<<+- DAMASK_spectral init -+>>>' write(6,'(a)') ' <<<+- DAMASK_spectral_Driver init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a)') ' Working Directory: ',trim(getSolverWorkingDirectoryName()) 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' if (debugGeneral) write(6,'(a)') 'Header of result file written out'
endif 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 ! Loop over loadcases defined in the currentLoadcase file
@ -390,7 +397,10 @@ print*, 'my Unit closed'
write(6,'(a)') '##################################################################' write(6,'(a)') '##################################################################'
write(6,'(A,I5.5,A,es12.5)') 'Increment ', totalIncsCounter, ' Time ',time 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, & guessmode,timeinc,timeinc_old, &
P_BC = bc(currentLoadcase)%stress, & P_BC = bc(currentLoadcase)%stress, &
F_BC = bc(currentLoadcase)%deformation, & F_BC = bc(currentLoadcase)%deformation, &
@ -399,6 +409,18 @@ print*, 'my Unit closed'
velgrad = bc(currentLoadcase)%velGradApplied, & velgrad = bc(currentLoadcase)%velGradApplied, &
rotation_BC = bc(currentLoadcase)%rotation) 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)') ''
write(6,'(a)') '==================================================================' write(6,'(a)') '=================================================================='
if(solres%converged) then if(solres%converged) then
@ -420,6 +442,15 @@ print*, 'my Unit closed'
enddo incLooping enddo incLooping
enddo loadCaseLooping 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,'(a)') '##################################################################' write(6,'(a)') '##################################################################'
write(6,'(i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & write(6,'(i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', &

259
code/DAMASK_spectral_SolverAL.f90
View File

@ -1,6 +1,263 @@
module DAMASK_spectral_SolverAL 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, 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
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 end module DAMASK_spectral_SolverAL

182
code/DAMASK_spectral_SolverBasic.f90
View File

@ -1,11 +1,112 @@
module DAMASK_spectral_SolverBasic 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 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 contains
subroutine Basic_Init() subroutine basic_init()
call Utilities_Init()
end 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) 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: & use FEsolving, only: &
restartWrite restartWrite
implicit none implicit none
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
real(pReal), intent(in) :: timeinc, timeinc_old real(pReal), intent(in) :: timeinc, timeinc_old
real(pReal), intent(in) :: guessmode 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 real(pReal), dimension(3,3), intent(in) :: P_BC,F_BC,rotation_BC
logical, dimension(9), intent(in) :: mask_stressVector 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
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! 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) :: iter
integer(pInt) :: i, j, k integer(pInt) :: i, j, k
logical :: ForwardResults logical :: ForwardData
real(pReal) :: defgradDet real(pReal) :: defgradDet
real(pReal) :: defgradDetMax, defgradDetMin real(pReal) :: defgradDetMax, defgradDetMin
mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3])) mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3]))
mask_defgrad = merge(zeroes,ones,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) close(777)
endif endif
ForwardResults = .True. ForwardData = .True.
if (velgrad) then ! calculate deltaF_aim from given L and current F if (velgrad) then ! calculate deltaF_aim from given L and current F
deltaF_aim = timeinc * mask_defgrad * math_mul33x33(F_BC, F_aim) deltaF_aim = timeinc * mask_defgrad * math_mul33x33(F_BC, F_aim)
else ! deltaF_aim = fDot *timeinc where applicable else ! deltaF_aim = fDot *timeinc where applicable
@ -102,7 +192,8 @@ real(pReal), dimension(3,3) :: &
1.0_pReal,F_lastInc,coordinates) 1.0_pReal,F_lastInc,coordinates)
iter = 0_pInt 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)))& convergenceLoop: do while((iter < itmax .and. (any([err_div ,err_stress] > 1.0_pReal)))&
.or. iter < itmin) .or. iter < itmin)
@ -118,23 +209,24 @@ S = S_lastInc(rotation_BC,mask_stressVector)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call constitutiveResponse(coordinates,F,F_lastInc,temperature,timeinc,&
P,C,P_av,ForwardData,rotation_BC)
ForwardData = .False.
call constitutiveResponse(ForwardResults,timeInc)
ForwardResults = .False.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
if(any(mask_stressVector)) then ! calculate stress BC if applied if(any(mask_stressVector)) then ! calculate stress BC if applied
err_stress = BCcorrection(mask_stressVector,P_BC,S) err_stress = BCcorrection(mask_stressVector,P_BC,P_av,F_aim,S)
else else
err_stress = 0.0_pReal err_stress = 0.0_pReal
endif endif
F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient ! updated deformation gradient
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) = 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) 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 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 function basic_solution
subroutine basic_destroy()
implicit none
call Utilities_destroy()
end subroutine basic_destroy
end module DAMASK_spectral_SolverBasic end module DAMASK_spectral_SolverBasic

608
code/DAMASK_spectral_Utilities.f90
View File

@ -47,20 +47,6 @@ module DAMASK_spectral_Utilities
implicit none 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 ! variables storing information for spectral method and FFTW
type(C_PTR) :: plan_forward, plan_backward ! plans for fftw type(C_PTR) :: plan_forward, plan_backward ! plans for fftw
@ -86,28 +72,22 @@ module DAMASK_spectral_Utilities
!variables controlling debugging !variables controlling debugging
logical :: debugGeneral, debugDivergence, debugRestart, debugFFTW 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 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 integer(pInt), dimension(3) :: res = 1_pInt
real(pReal) :: wgt real(pReal) :: wgt
integer(pInt) :: res1_red, Npoints integer(pInt) :: res1_red, Npoints
!--------------------------------------------------------------------------------------------------
! solution state
type solutionState
logical :: converged = .false.
logical :: regrid = .false.
logical :: term_ill = .false.
end type solutionState
contains contains
subroutine Utilities_init(F,P,F_...) subroutine Utilities_init(C_ref)
use DAMASK_interface, only: &
getSolverJobName
use mesh, only : & use mesh, only : &
mesh_spectral_getResolution, & mesh_spectral_getResolution, &
@ -127,23 +107,16 @@ subroutine Utilities_init(F,P,F_...)
debug_spectralRestart, & debug_spectralRestart, &
debug_spectralFFTW debug_spectralFFTW
use FEsolving, only: &
restartInc
use numerics, only: & use numerics, only: &
memory_efficient memory_efficient
use CPFEM, only: &
CPFEM_general
use IO, only: &
IO_read_JobBinaryFile, &
IO_write_JobBinaryFile
implicit none implicit none
real(pReal), dimension(3,3) :: temp33_Real, xiDyad real(pReal), dimension(3,3) :: temp33_Real, xiDyad
integer(pInt) :: i, j, k, l, m, n, q, ierr integer(pInt) :: i, j, k, l, m, n, q, ierr
integer(pInt), dimension(3) :: k_s 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) :: tensorField ! field in real and fourier space
type(C_PTR) :: scalarField_realC, scalarField_fourierC type(C_PTR) :: scalarField_realC, scalarField_fourierC
@ -151,7 +124,7 @@ subroutine Utilities_init(F,P,F_...)
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') ' <<<+- DAMASK_spectralSolver init -+>>>' write(6,'(a)') ' <<<+- DAMASK_spectralSolver Utilities init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a)') '' write(6,'(a)') ''
@ -172,12 +145,7 @@ subroutine Utilities_init(F,P,F_...)
Npoints = res(1)*res(2)*res(3) Npoints = res(1)*res(2)*res(3)
wgt = 1.0/real(Npoints,pReal) wgt = 1.0/real(Npoints,pReal)
allocate (F ( 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) ! start out isothermally
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
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) 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_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 call c_f_pointer(tensorField, field_fourier, [ res1_red, res(2),res(3),3,3]) ! place a pointer for a complex representation on tensorField
@ -233,37 +201,6 @@ subroutine Utilities_init(F,P,F_...)
if (debugGeneral) write(6,'(a)') 'FFTW initialized' 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) ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around)
if (divergence_correction) then if (divergence_correction) then
@ -285,19 +222,6 @@ subroutine Utilities_init(F,P,F_...)
xi(1:3,i,j,k) = real(k_s, pReal)/virt_dim xi(1:3,i,j,k) = real(k_s, pReal)/virt_dim
enddo; enddo; enddo 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 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(1,1,1,3,3,3,3), source = 0.0_pReal)
else ! precalculation of gamma_hat field else ! precalculation of gamma_hat field
@ -317,252 +241,250 @@ subroutine Utilities_init(F,P,F_...)
endif endif
end subroutine Utilities_init end subroutine Utilities_init
real(pReal) function convolution(calcDivergence, field_aim,) real(pReal) function convolution(calcDivergence,field_aim,C_ref)
use numerics, only: &
memory_efficient, & use numerics, only: &
err_div_tol memory_efficient, &
real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors err_div_tol
real(pReal) :: err_div = 0.0_pReal
real(pReal), dimension(3,3) :: temp33_Real real(pReal), dimension(3,3) :: xiDyad ! product of wave vectors
integer(pInt) :: i, j, k, l, m, n, q 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 !variables for additional output due to general debugging
real(pReal) :: maxCorrectionSym, maxCorrectionSkew real(pReal) :: maxCorrectionSym, maxCorrectionSkew
logical :: calcDivergence logical :: calcDivergence
real(pReal), dimension(3,3) :: field_avg, field_aim real(pReal), dimension(3,3) :: field_avg, field_aim
integer(pInt) :: row, column integer(pInt) :: row, column
real(pReal) :: field_av_L2, err_div_RMS, err_real_div_RMS, err_post_div_RMS,& real(pReal) :: field_av_L2, err_div_RMS, err_real_div_RMS, err_post_div_RMS,&
err_div_max, err_real_div_max err_div_max, err_real_div_max
complex(pReal), dimension(3) :: temp3_complex complex(pReal), dimension(3) :: temp3_complex
complex(pReal), dimension(3,3) :: temp33_complex complex(pReal), dimension(3,3) :: temp33_complex
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! actual spectral method ! actual spectral method
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '... doing convolution .................' write(6,'(a)') '... doing convolution .................'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! copy one component of the stress field to to a single FT and check for mismatch ! copy one component of the stress field to to a single FT and check for mismatch
if (debugFFTW) then 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 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 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 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) cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! call function to calculate divergence from math (for post processing) to check results ! call function to calculate divergence from math (for post processing) to check results
if (debugDivergence) & if (debugDivergence) &
call divergence_fft(res,virt_dim,3_pInt,& 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 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 ! 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 ! comparing 1 and 3x3 FT results
if (debugFFTW) then if (debugFFTW) then
call fftw_execute_dft(plan_scalarField_forth,scalarField_real,scalarField_fourier) 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,i1,1x,i1)') 'checking FT results of compontent ', row, column
write(6,'(a,2(es11.4,1x))') 'max FT relative error = ',& write(6,'(a,2(es11.4,1x))') 'max FT relative error = ',&
maxval( real((scalarField_fourier(1:res1_red,1:res(2),1:res(3))-& 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))/& field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), & scalarField_fourier(1:res1_red,1:res(2),1:res(3)))), &
maxval(aimag((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))/& field_fourier(1:res1_red,1:res(2),1:res(3),row,column))/&
scalarField_fourier(1:res1_red,1:res(2),1:res(3)))) scalarField_fourier(1:res1_red,1:res(2),1:res(3))))
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! removing highest frequencies ! removing highest frequencies
field_fourier ( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)& field_fourier ( res1_red,1:res(2) , 1:res(3) ,1:3,1:3)&
= cmplx(0.0_pReal,0.0_pReal,pReal) = 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)& 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) = cmplx(0.0_pReal,0.0_pReal,pReal)
if(res(3)>1_pInt) & if(res(3)>1_pInt) &
field_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& 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) = cmplx(0.0_pReal,0.0_pReal,pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space ! calculating RMS divergence criterion in Fourier space
if( calcDivergence) then if(calcDivergence) then
field_avg = real(field_fourier(1,1,1,1:3,1:3),pReal)*wgt 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) 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))))) math_transpose33(field_avg)))))
err_div_RMS = 0.0_pReal err_div_RMS = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2) 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. do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
err_div_RMS = err_div_RMS & 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 + 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 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),& +sum(aimag(math_mul33x3_complex(field_fourier(i,j,k,1:3,1:3),&
xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal)) xi(1:3,i,j,k))*TWOPIIMG)**2.0_pReal))
enddo enddo
err_div_RMS = err_div_RMS & ! Those two layers (DC and Nyquist) do not have a conjugate complex counterpart 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),& + sum( real(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)&
+ sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),& + sum(aimag(math_mul33x3_complex(field_fourier(1 ,j,k,1:3,1:3),&
xi(1:3,1 ,j,k))*TWOPIIMG)**2.0_pReal)& 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),& + 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)& 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),& + 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) xi(1:3,res1_red,j,k))*TWOPIIMG)**2.0_pReal)
enddo; enddo enddo; enddo
err_div_RMS = sqrt(err_div_RMS)*wgt ! RMS in real space calculated with Parsevals theorem from Fourier space 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 err_div = err_div_RMS/field_av_L2 ! criterion to stop iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate additional divergence criteria and report ! calculate additional divergence criteria and report
if (debugDivergence) then ! calculate divergence again if (debugDivergence) then ! calculate divergence again
err_div_max = 0.0_pReal err_div_max = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red 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 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 xi(1:3,i,j,k))*TWOPIIMG
err_div_max = max(err_div_max,sum(abs(temp3_Complex)**2.0_pReal)) 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 divergence_fourier(i,j,k,1:3) = temp3_Complex ! need divergence NOT squared
enddo; enddo; enddo enddo; enddo; enddo
call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real) ! already weighted call fftw_execute_dft_c2r(plan_divergence,divergence_fourier,divergence_real) ! already weighted
err_real_div_RMS = 0.0_pReal err_real_div_RMS = 0.0_pReal
err_post_div_RMS = 0.0_pReal err_post_div_RMS = 0.0_pReal
err_real_div_max = 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) 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_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_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 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 enddo; enddo; enddo
err_real_div_RMS = sqrt(wgt*err_real_div_RMS) ! RMS in real space 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_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_max = sqrt( err_real_div_max) ! max in real space
err_div_max = sqrt( err_div_max) ! max in Fourier 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 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 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 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 FT max = ',err_div_max
write(6,'(a,es11.4)') 'error divergence Real max = ',err_real_div_max write(6,'(a,es11.4)') 'error divergence Real max = ',err_real_div_max
endif endif
write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/err_div_tol,& write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/err_div_tol,&
' (',err_div,' N/m³)' ' (',err_div,' N/m³)'
end if end if
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! to the actual spectral method calculation (mechanical equilibrium) ! to the actual spectral method calculation (mechanical equilibrium)
if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat 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
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
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) &
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k)
xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
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(l,m) = sum(C_ref(l,m,1:3,1:3)*xiDyad) temp33_Real = math_inv33(temp33_Real)
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)&
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)
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) &
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) *&
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)) field_fourier(i,j,k,1:3,1:3) = temp33_Complex
field_fourier(i,j,k,1:3,1:3) = temp33_Complex endif
endif enddo; enddo; enddo
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
else ! use precalculated gamma-operator 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) *&
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt,res1_red field_fourier(i,j,k,1:3,1:3))
forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt) & field_fourier(i,j,k, 1:3,1:3) = temp33_Complex
temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) *& enddo; enddo; enddo
field_fourier(i,j,k,1:3,1:3)) endif
field_fourier(i,j,k, 1:3,1:3) = temp33_Complex 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
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 ! comparing 1 and 3x3 inverse FT results
if (debugFFTW) then if (debugFFTW) then
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red 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) scalarField_fourier(i,j,k) = field_fourier(i,j,k,row,column)
enddo; enddo; enddo enddo; enddo; enddo
do i = 0_pInt, res(1)/2_pInt-2_pInt ! unpack fft data for conj complex symmetric part do i = 0_pInt, res(1)/2_pInt-2_pInt ! unpack fft data for conj complex symmetric part
m = 1_pInt m = 1_pInt
do k = 1_pInt, res(3) do k = 1_pInt, res(3)
n = 1_pInt n = 1_pInt
do j = 1_pInt, res(2) do j = 1_pInt, res(2)
scalarField_fourier(res(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m)) scalarField_fourier(res(1)-i,j,k) = conjg(scalarField_fourier(2+i,n,m))
if(n == 1_pInt) n = res(2) + 1_pInt if(n == 1_pInt) n = res(2) + 1_pInt
n = n-1_pInt n = n-1_pInt
enddo enddo
if(m == 1_pInt) m = res(3) + 1_pInt if(m == 1_pInt) m = res(3) + 1_pInt
m = m -1_pInt m = m -1_pInt
enddo; enddo enddo; enddo
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing the inverse FT ! 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 ! comparing 1 and 3x3 inverse FT results
if (debugFFTW) then if (debugFFTW) then
write(6,'(a,i1,1x,i1)') 'checking iFT results of compontent ', row, column write(6,'(a,i1,1x,i1)') 'checking iFT results of compontent ', row, column
call fftw_execute_dft(plan_scalarField_back,scalarField_fourier,scalarField_real) call fftw_execute_dft(plan_scalarField_back,scalarField_fourier,scalarField_real)
write(6,'(a,es11.4)') 'max iFT relative error = ',& write(6,'(a,es11.4)') 'max iFT relative error = ',&
maxval((real(scalarField_real(1:res(1),1:res(2),1:res(3)))-& 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))/& 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)))) real(scalarField_real(1:res(1),1:res(2),1:res(3))))
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate some additional output ! calculate some additional output
if(debugGeneral) then if(debugGeneral) then
maxCorrectionSkew = 0.0_pReal maxCorrectionSkew = 0.0_pReal
maxCorrectionSym = 0.0_pReal maxCorrectionSym = 0.0_pReal
temp33_Real = 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) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
maxCorrectionSym = max(maxCorrectionSym,& maxCorrectionSym = max(maxCorrectionSym,&
maxval(math_symmetric33(field_real(i,j,k,1:3,1:3)))) maxval(math_symmetric33(field_real(i,j,k,1:3,1:3))))
maxCorrectionSkew = max(maxCorrectionSkew,& maxCorrectionSkew = max(maxCorrectionSkew,&
maxval(math_skew33(field_real(i,j,k,1:3,1:3)))) 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) temp33_Real = temp33_Real + field_real(i,j,k,1:3,1:3)
enddo; enddo; enddo enddo; enddo; enddo
write(6,'(a,1x,es11.4)') 'max symmetric correction of deformation =',& write(6,'(a,1x,es11.4)') 'max symmetric correction of deformation =',&
maxCorrectionSym*wgt maxCorrectionSym*wgt
write(6,'(a,1x,es11.4)') 'max skew correction of deformation =',& write(6,'(a,1x,es11.4)') 'max skew correction of deformation =',&
maxCorrectionSkew*wgt maxCorrectionSkew*wgt
write(6,'(a,1x,es11.4)') 'max sym/skew of avg correction = ',& write(6,'(a,1x,es11.4)') 'max sym/skew of avg correction = ',&
maxval(math_symmetric33(temp33_real))/& maxval(math_symmetric33(temp33_real))/&
maxval(math_skew33(temp33_real)) maxval(math_skew33(temp33_real))
endif endif
field_real = field_real * wgt field_real = field_real * wgt
convolution = err_div/err_div_tol convolution = err_div/err_div_tol
end function convolution 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 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 real(pReal), dimension(3,3), intent(in) :: rot_BC
logical, dimension(9), intent(in) :: mask_stressVector1 logical, dimension(9), intent(in) :: mask_stressVector1
real(pReal), dimension(3,3,3,3) :: C_lastInc 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 integer(pInt) :: size_reduced = 0_pInt
real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC) real(pReal), dimension(:,:), allocatable :: s_reduced, c_reduced ! reduced compliance and stiffness (only for stress BC)
logical :: errmatinv logical :: errmatinv
size_reduced = count(mask_stressVector1) 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 (c_reduced(size_reduced,size_reduced), source =0.0_pReal)
allocate (s_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 C_lastInc = math_rotate_forward3333(C,rot_BC) ! calculate stiffness from former inc
temp99_Real = math_Plain3333to99(C_lastInc) temp99_Real = math_Plain3333to99(C_lastInc)
k = 0_pInt ! build reduced stiffness k = 0_pInt ! build reduced stiffness
@ -588,9 +510,7 @@ function S_lastInc(rot_BC,mask_stressVector1)
j = j + 1_pInt j = j + 1_pInt
temp99_Real(n,m) = s_reduced(k,j) temp99_Real(n,m) = s_reduced(k,j)
endif; enddo; endif; enddo endif; enddo; endif; enddo
S_lastInc = (math_Plain99to3333(temp99_Real)) S_lastInc = math_Plain99to3333(temp99_Real)
if (allocated(c_reduced)) deallocate(c_reduced)
if (allocated(c_reduced)) deallocate(c_reduced)
end function S_lastInc end function S_lastInc
@ -598,76 +518,90 @@ end function S_lastInc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate reduced compliance ! calculate reduced compliance
real(pReal) function BCcorrection(mask_stressVector,P_BC,S_lastInc) 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 use numerics, only: err_stress_tolrel, err_stress_tolabs
real(pReal) :: err_stress, err_stress_tol
real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal logical, dimension(9) :: mask_stressVector
real(pReal), dimension(3,3,3,3) :: S_lastInc real(pReal) :: err_stress, err_stress_tol
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal
P_BC , & real(pReal), dimension(3,3,3,3) :: S_lastInc
mask_stress, & real(pReal), dimension(3,3) :: &
mask_defgrad P_BC , &
mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3])) P_av, &
mask_defgrad = merge(zeroes,ones,reshape(mask_stressVector,[3,3])) 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 ! stress BC handling
! calculate stress BC if applied ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (P_av - P_BC))) ! maximum deviaton (tensor norm not applicable) 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 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)') ''
write(6,'(a)') '... correcting deformation gradient to fulfill BCs ...............' 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, & write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, &
' (',err_stress,' Pa)' ' (',err_stress,' Pa)'
F_aim = F_aim - math_mul3333xx33(S_lastInc, ((P_av - P_BC))) ! residual on given stress components 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) write(6,'(a,1x,es11.4)')'determinant of new deformation = ',math_det33(F_aim)
BCcorrection = err_stress/err_stress_tol BCcorrection = err_stress/err_stress_tol
end function BCcorrection end function BCcorrection
subroutine constitutiveResponse(F,P,ForwardData,timeinc) subroutine constitutiveResponse(coordinates,F,F_lastInc,temperature,timeinc,&
use debug, only: & P,C,P_av,ForwardData,rotation_BC)
debug_reset, & use debug, only: &
debug_info debug_reset, &
use CPFEM, only: & debug_info
CPFEM_general use CPFEM, only: &
use FEsolving, only: restartWrite 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 real(pReal) :: timeinc
logical :: ForwardData logical :: ForwardData
integer(pInt) :: i, j, k, ielem integer(pInt) :: i, j, k, ielem
integer(pInt) :: CPFEM_mode integer(pInt) :: CPFEM_mode
real(pReal), dimension(3,3,3,3) :: dPdF real(pReal), dimension(3,3,3,3) :: dPdF, C
real(pReal), dimension(6) :: sigma ! cauchy stress real(pReal), dimension(6) :: sigma ! cauchy stress
real(pReal), dimension(6,6) :: dsde real(pReal), dimension(6,6) :: dsde
real(pReal), dimension(3,3) :: P_av_lab, rotation_BC real(pReal), dimension(3,3) :: P_av_lab, P_av, rotation_BC
if (ForwardData) then if (ForwardData) then
CPFEM_mode = 1_pInt CPFEM_mode = 1_pInt
else else
CPFEM_mode = 2_pInt CPFEM_mode = 2_pInt
endif endif
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '... update stress field P(F) .....................................' write(6,'(a)') '... update stress field P(F) .....................................'
ielem = 0_pInt ielem = 0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(3_pInt,& ! collect cycle 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), & 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) 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 P = 0.0_pReal ! needed because of the padding for FFTW
C = 0.0_pReal C = 0.0_pReal
P_av_lab = 0.0_pReal P_av_lab = 0.0_pReal
ielem = 0_pInt ielem = 0_pInt
call debug_reset() call debug_reset()
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
ielem = ielem + 1_pInt ielem = ielem + 1_pInt
call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, 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) 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) temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(i,j,k,1:3,1:3),dPdF)
CPFEM_mode = 2_pInt CPFEM_mode = 2_pInt
C = C + dPdF C = C + dPdF
P_av_lab = P_av_lab + P(i,j,k,1:3,1:3) P_av_lab = P_av_lab + P(i,j,k,1:3,1:3)
enddo; enddo; enddo enddo; enddo; enddo
call debug_info() call debug_info()
restartWrite = .false. restartWrite = .false.
@ -678,4 +612,22 @@ use FEsolving, only: restartWrite
C = C * wgt C = C * wgt
end subroutine constitutiveResponse 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 end module DAMASK_spectral_Utilities

8
code/numerics.f90
View File

@ -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 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 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 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 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 itmax = 20_pInt, & ! maximum number of iterations
itmin = 2_pInt ! minimum 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) fftw_timelimit = IO_floatValue(line,positions,2_pInt)
case ('fftw_plan_mode') case ('fftw_plan_mode')
fftw_plan_mode = IO_stringValue(line,positions,2_pInt) fftw_plan_mode = IO_stringValue(line,positions,2_pInt)
case ('myspectralsolver')
myspectralsolver = IO_stringValue(line,positions,2_pInt)
case ('rotation_tol') case ('rotation_tol')
rotation_tol = IO_floatValue(line,positions,2_pInt) rotation_tol = IO_floatValue(line,positions,2_pInt)
case ('divergence_correction') case ('divergence_correction')
@ -256,7 +259,7 @@ subroutine numerics_init
#endif #endif
#ifndef Spectral #ifndef Spectral
case ('err_div_tol','err_stress_tolrel','err_stress_tolabs',& 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') 'rotation_tol','divergence_correction','update_gamma')
call IO_warning(40_pInt,ext_msg=tag) call IO_warning(40_pInt,ext_msg=tag)
#endif #endif
@ -348,6 +351,7 @@ subroutine numerics_init
write(6,'(a24,1x,es8.1)') ' fftw_timelimit: ',fftw_timelimit write(6,'(a24,1x,es8.1)') ' fftw_timelimit: ',fftw_timelimit
endif endif
write(6,'(a24,1x,a)') ' fftw_plan_mode: ',trim(fftw_plan_mode) 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,i8)') ' fftw_planner_flag: ',fftw_planner_flag
write(6,'(a24,1x,es8.1)') ' rotation_tol: ',rotation_tol write(6,'(a24,1x,es8.1)') ' rotation_tol: ',rotation_tol
write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction write(6,'(a24,1x,L8,/)') ' divergence_correction: ',divergence_correction