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

added comments and structured the code. Temperature is not longer stored for each point but only to simplify restart behavior as long as it is not fully supported by the constitutive models

This commit is contained in:
Martin Diehl 2012-11-12 14:14:39 +00:00
parent 7974001c9d
commit 70c4e11742
4 changed files with 774 additions and 751 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

740
code/DAMASK_spectral_solverAL.f90
View File

@ -1,20 +1,18 @@
!--------------------------------------------------------------------------------------------------
! $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $
! $Id: DAMASK_spectral_solverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief AL scheme solver
!--------------------------------------------------------------------------------------------------
module DAMASK_spectral_SolverAL
module DAMASK_spectral_solverAL
use prec, only: &
pInt, &
pReal
use math, only: &
math_I3
use DAMASK_spectral_Utilities, only: &
use DAMASK_spectral_utilities, only: &
tSolutionState
implicit none
@ -26,13 +24,14 @@ module DAMASK_spectral_SolverAL
DAMASK_spectral_SolverAL_label = 'al'
!--------------------------------------------------------------------------------------------------
! derived types
! derived types ToDo: use here the type definition for a full loadcase including mask
type tSolutionParams
real(pReal), dimension(3,3) :: P_BC, rotation_BC
real(pReal) :: timeinc
end type tSolutionParams
type(tSolutionParams), private :: params
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! PETSc data
@ -42,27 +41,34 @@ module DAMASK_spectral_SolverAL
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc, F_lambdaDot, Fdot
real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, & !< field of previous compatible deformation gradients
F_lambda_lastInc, & !< field of previous incompatible deformation gradient
Fdot, & !< field of assumed rate of compatible deformation gradient
F_lambdaDot !< field of assumed rate of incopatible deformation gradient
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private, dimension(:,:,:), allocatable :: temperature
real(pReal), private :: temperature !< temperature, no spatial quantity at the moment
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: &
F_aimDot, &
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
P_av
character(len=1024), private :: incInfo
F_aimDot, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
C = 0.0_pReal, C_lastInc = 0.0_pReal, &
S = 0.0_pReal, &
C_scale = 0.0_pReal, &
C = 0.0_pReal, & !< current average stiffness
C_lastInc = 0.0_pReal, & !< previous average stiffness
S = 0.0_pReal, & !< current compliance (filled up with zeros)
C_scale = 0.0_pReal, &
S_scale = 0.0_pReal
real(pReal), private :: err_stress, err_f, err_p
logical, private :: ForwardData
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
real(pReal), private :: &
err_stress, & !< deviation from stress BC
err_f, & !< difference between compatible and incompatible deformation gradient
err_p !< difference of stress resulting from compatible and incompatible F
logical, private :: ForwardData
integer(pInt) :: reportIter = 0_pInt
contains
@ -70,215 +76,209 @@ contains
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
subroutine AL_init()
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO, only: &
IO_read_JobBinaryFile, &
IO_write_JobBinaryFile
use FEsolving, only: &
restartInc
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment)
use IO, only: &
IO_read_JobBinaryFile, &
IO_write_JobBinaryFile
use FEsolving, only: &
restartInc
use DAMASK_interface, only: &
getSolverJobName
use DAMASK_spectral_Utilities, only: &
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugRestart
use numerics, only: &
petsc_options
use mesh, only: &
res, &
geomdim, &
mesh_NcpElems
use math, only: &
math_invSym3333
use DAMASK_interface, only: &
getSolverJobName
use DAMASK_spectral_Utilities, only: &
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugRestart
use numerics, only: &
petsc_options
use mesh, only: &
res, &
geomdim, &
mesh_NcpElems
use math, only: &
math_invSym3333
implicit none
implicit none
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90>
integer(pInt) :: i,j,k
real(pReal), dimension(:,:,:,:,:), allocatable :: P
PetscErrorCode :: ierr
PetscObject :: dummy
PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda
call Utilities_init()
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>'
write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $'
#include "compilation_info.f90"
write(6,'(a)') ''
allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
! allocate (Fdot,source = F_lastInc) somethin like that should be possible
allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lambdaDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (P (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! PETSc Init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
call DMCreateGlobalVector(da,solution_vec,ierr)
call DMDASetLocalFunction(da,AL_formResidual,ierr)
call SNESSetDM(snes,da,ierr)
call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr)
call SNESSetFromOptions(snes,ierr)
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F_lambda_lastInc = F_lastInc
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
F_lambda = F
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo
elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
restartInc - 1_pInt,' from file'
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) 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,'convergedSpectralDefgradLambda',&
trim(getSolverJobName()),size(F_lambda_lastInc))
read (777,rec=1) F_lambda
close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',&
trim(getSolverJobName()),size(F_lambda_lastInc))
read (777,rec=1) F_lambda_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
!print*, shape(F)
call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr)
!--------------------------------------------------------------------------------------------------
! reference stiffness
if (restartInc == 1_pInt) then
call IO_write_jobBinaryFile(777,'C_ref',size(C))
write (777,rec=1) C
close(777)
elseif (restartInc > 1_pInt) then
call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C))
read (777,rec=1) C
close (777)
endif
call Utilities_updateGamma(C,.True.)
C_scale = C
S_scale = math_invSym3333(C)
integer(pInt) :: i,j,k
real(pReal), dimension(3,3, res(1), res(2),res(3)) :: P
end subroutine AL_init
PetscErrorCode :: ierr
PetscObject :: dummy
PetscScalar, pointer, dimension(:,:,:,:) :: xx_psc, F, F_lambda
call Utilities_init()
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>'
write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $'
#include "compilation_info.f90"
write(6,'(a)') ''
allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
! allocate (Fdot,source = F_lastInc) somethin like that should be possible
allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lambdaDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! PETSc Init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
CHKERRQ(ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr)
CHKERRQ(ierr)
call DMDASetLocalFunction(da,AL_formResidual,ierr)
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr)
CHKERRQ(ierr)
call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr)
call SNESSetFromOptions(snes,ierr)
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr) ! places pointer xx_psc on PETSc data
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F_lambda_lastInc = F_lastInc
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
F_lambda = F
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo
elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
restartInc - 1_pInt,' from file'
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) 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,'convergedSpectralDefgradLambda',&
trim(getSolverJobName()),size(F_lambda_lastInc))
read (777,rec=1) F_lambda
close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',&
trim(getSolverJobName()),size(F_lambda_lastInc))
read (777,rec=1) F_lambda_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 Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,P_av,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr)
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! reference stiffness
if (restartInc == 1_pInt) then
call IO_write_jobBinaryFile(777,'C_ref',size(C))
write (777,rec=1) C
close(777)
elseif (restartInc > 1_pInt) then
call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C))
read (777,rec=1) C
close (777)
endif
call Utilities_updateGamma(C,.True.)
C_scale = C
S_scale = math_invSym3333(C)
end subroutine AL_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the AL scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function &
type(tSolutionState) function &
AL_solution(incInfoIn,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
use numerics, only: &
update_gamma
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use mesh, only: &
res,&
geomdim,&
deformed_fft
use IO, only: &
IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: &
tBoundaryCondition, &
Utilities_forwardField, &
Utilities_calculateRate, &
Utilities_maskedCompliance, &
Utilities_updateGamma, &
cutBack
use FEsolving, only: &
restartWrite, &
terminallyIll
implicit none
use numerics, only: &
update_gamma
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use mesh, only: &
res,&
geomdim,&
deformed_fft
use IO, only: &
IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: &
tBoundaryCondition, &
Utilities_forwardField, &
Utilities_calculateRate, &
Utilities_maskedCompliance, &
Utilities_updateGamma, &
cutBack
use FEsolving, only: &
restartWrite, &
terminallyIll
implicit none
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90>
!--------------------------------------------------------------------------------------------------
! input data for solution
real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc
logical, intent(in) :: guess
type(tBoundaryCondition), intent(in) :: P_BC,F_BC
character(len=*), intent(in) :: incInfoIn
real(pReal), dimension(3,3), intent(in) :: rotation_BC
real(pReal), dimension(3,3) ,save :: F_aimDot
real(pReal), dimension(3,3) :: F_aim_lab
real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc
logical, intent(in) :: guess
type(tBoundaryCondition), intent(in) :: P_BC,F_BC
character(len=*), intent(in) :: incInfoIn
real(pReal), dimension(3,3), intent(in) :: rotation_BC
real(pReal), dimension(3,3) ,save :: F_aimDot
real(pReal), dimension(3,3) :: F_aim_lab
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
real(pReal), dimension(3,3) :: temp33_Real
real(pReal), dimension(3,3) :: temp33_Real
!--------------------------------------------------------------------------------------------------
!
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda
PetscErrorCode :: ierr
SNESConvergedReason ::reason
! PETSc Data
PetscScalar, dimension(:,:,:,:), pointer :: xx_psc, F, F_lambda
PetscErrorCode :: ierr
SNESConvergedReason ::reason
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
incInfo = incInfoIn
if (restartWrite) then
write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc))
write (777,rec=1) F_LastInc
close (777)
call IO_write_jobBinaryFile(777,'C',size(C))
write (777,rec=1) C
close(777)
endif
AL_solution%converged =.false.
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
incInfo = incInfoIn
if (restartWrite) then
write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc))
write (777,rec=1) F_LastInc
close (777)
call IO_write_jobBinaryFile(777,'C',size(C))
write (777,rec=1) C
close(777)
endif
AL_solution%converged =.false.
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr)
F => xx_psc(0:8,:,:,:)
F_lambda => xx_psc(9:17,:,:,:)
if ( cutBack) then
F_aim = F_aim_lastInc
F_lambda= reshape(F_lambda_lastInc,[9,res(1),res(2),res(3)])
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
C = C_lastInc
else
!--------------------------------------------------------------------------------------------------
if ( cutBack) then
F_aim = F_aim_lastInc
F_lambda= reshape(F_lambda_lastInc,[9,res(1),res(2),res(3)])
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
C = C_lastInc
else
C_lastInc = C
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
@ -287,7 +287,7 @@ else
elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed
f_aimDot = F_BC%maskFloat * F_BC%values
endif
if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
@ -300,94 +300,97 @@ else
F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)])
F_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)])
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), &
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), &
1.0_pReal,F_lastInc,coordinates)
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
! update local deformation gradient and coordinates
! deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC)
! deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC)
F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)])
F_lambda = reshape(Utilities_forwardField(timeinc,F_aim,F_lambda_lastInc,F_lambdadot),[9,res(1),res(2),res(3)])
F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)])
F_lambda = reshape(Utilities_forwardField(timeinc,F_aim,F_lambda_lastInc,F_lambdadot),[9,res(1),res(2),res(3)])
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr)
CHKERRQ(ierr)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr)
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C)
if (update_gamma) call Utilities_updateGamma(C,restartWrite)
ForwardData = .True.
mask_stress = P_BC%maskFloat
params%P_BC = P_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C)
if (update_gamma) call Utilities_updateGamma(C,restartWrite)
ForwardData = .True.
mask_stress = P_BC%maskFloat
params%P_BC = P_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr)
call SNESGetConvergedReason(snes,reason,ierr)
AL_solution%termIll = terminallyIll
terminallyIll = .false.
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr)
CHKERRQ(ierr)
call SNESGetConvergedReason(snes,reason,ierr)
CHKERRQ(ierr)
AL_solution%termIll = terminallyIll
terminallyIll = .false.
if (reason > 0 ) then
AL_solution%converged = .true.
AL_solution%iterationsNeeded = reportIter
endif
end function AL_solution
if (reason > 0 ) then
AL_solution%converged = .true.
AL_solution%iterationsNeeded = reportIter
endif
end function AL_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the AL residual vector
!--------------------------------------------------------------------------------------------------
subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: &
itmax, &
itmin
use math, only: &
math_rotate_backward33, &
math_transpose33, &
math_mul3333xx33
use mesh, only: &
res, &
wgt
use DAMASK_spectral_Utilities, only: &
field_real, &
Utilities_FFTforward, &
Utilities_fourierConvolution, &
Utilities_FFTbackward, &
Utilities_constitutiveResponse
use IO, only: IO_intOut
implicit none
subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
integer(pInt) :: i,j,k
integer(pInt), save :: callNo = 3_pInt
real(pReal), dimension(3,3) :: temp33_Real
logical :: report
DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE)
PetscScalar, target :: x_scal(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE)
PetscScalar, target :: f_scal(3,3,2,X_RANGE,Y_RANGE,Z_RANGE)
PetscScalar, pointer :: F(:,:,:,:,:), F_lambda(:,:,:,:,:)
PetscScalar, pointer :: residual_F(:,:,:,:,:), residual_F_lambda(:,:,:,:,:)
PetscInt :: iter, nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
use numerics, only: &
itmax, &
itmin
use math, only: &
math_rotate_backward33, &
math_transpose33, &
math_mul3333xx33
use mesh, only: &
res, &
wgt
use DAMASK_spectral_Utilities, only: &
field_real, &
Utilities_FFTforward, &
Utilities_fourierConvolution, &
Utilities_FFTbackward, &
Utilities_constitutiveResponse
use IO, only: IO_intOut
F => x_scal(:,:,1,:,:,:)
F_lambda => x_scal(:,:,2,:,:,:)
residual_F => f_scal(:,:,1,:,:,:)
residual_F_lambda => f_scal(:,:,2,:,:,:)
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr)
call SNESGetIterationNumber(snes,iter,ierr)
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
implicit none
integer(pInt) :: i,j,k
integer(pInt), save :: callNo = 3_pInt
real(pReal), dimension(3,3) :: temp33_Real
logical :: report
DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE)
PetscScalar, target :: x_scal(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE)
PetscScalar, target :: f_scal(3,3,2,X_RANGE,Y_RANGE,Z_RANGE)
PetscScalar, pointer :: F(:,:,:,:,:), F_lambda(:,:,:,:,:)
PetscScalar, pointer :: residual_F(:,:,:,:,:), residual_F_lambda(:,:,:,:,:)
PetscInt :: iter, nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
F => x_scal(:,:,1,:,:,:)
F_lambda => x_scal(:,:,2,:,:,:)
residual_F => f_scal(:,:,1,:,:,:)
residual_F_lambda => f_scal(:,:,2,:,:,:)
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr)
CHKERRQ(ierr)
call SNESGetIterationNumber(snes,iter,ierr)
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
if (iter == 0 .and. callNo>2) then
callNo = 0_pInt
reportIter = 0_pInt
@ -400,105 +403,112 @@ else
reportIter = reportIter + 1_pInt
endif
callNo = callNo +1_pInt
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, &
residual_F,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .False.
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, &
residual_F,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .False.
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc
err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc
err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc
!--------------------------------------------------------------------------------------------------
!
field_real = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
temp33_Real = math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) + math_I3
residual_F(1:3,1:3,i,j,k) = temp33_Real
field_real(i,j,k,1:3,1:3) = -math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k)-F(1:3,1:3,i,j,k))
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space
call Utilities_FFTforward()
call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC))
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_lambda = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),&
[3,3,res(1),res(2),res(3)],order=[3,4,5,1,2])
residual_F = residual_F - F_lambda + residual_F_lambda
!--------------------------------------------------------------------------------------------------
! calculating errors
err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal))
err_p = wgt*sqrt(sum((residual_F - residual_F_lambda)**2.0_pReal))
end subroutine AL_formResidual
!--------------------------------------------------------------------------------------------------
! doing Fourier transform
field_real = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
temp33_Real = math_mul3333xx33(S_scale,residual_F(1:3,1:3,i,j,k)) + math_I3
residual_F(1:3,1:3,i,j,k) = temp33_Real
field_real(i,j,k,1:3,1:3) = -math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k)-F(1:3,1:3,i,j,k))
enddo; enddo; enddo
!--------------------------------------------------------------------------------------------------
! doing Fourier transform
call Utilities_FFTforward()
call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC))
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_lambda = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),&
[3,3,res(1),res(2),res(3)],order=[3,4,5,1,2])
residual_F = residual_F - F_lambda + residual_F_lambda
!--------------------------------------------------------------------------------------------------
! calculating errors
err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal))
err_p = wgt*sqrt(sum((residual_F - residual_F_lambda)**2.0_pReal))
end subroutine AL_formResidual
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
err_f_tol, &
err_p_tol, &
err_stress_tolrel, &
err_stress_tolabs
implicit none
subroutine AL_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr)
SNES :: snes_local
PetscInt :: it
PetscReal :: xnorm, snorm, fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode ::ierr
logical :: Converged
Converged = (it > itmin .and. &
all([ err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol, &
err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol, &
err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal))
if (Converged) then
reason = 1
elseif (it > itmax) then
reason = -1
else
reason = 0
endif
write(6,'(a,es14.7)') 'error stress BC = ', err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)
write(6,'(a,es14.7)') 'error F = ', err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol
write(6,'(a,es14.7)') 'error P = ', err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol
write(6,'(/,a)') '=========================================================================='
end subroutine AL_converged
use numerics, only: &
itmax, &
itmin, &
err_f_tol, &
err_p_tol, &
err_stress_tolrel, &
err_stress_tolabs
implicit none
SNES :: snes_local
PetscInt :: it
PetscReal :: xnorm, snorm, fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode ::ierr
logical :: Converged
Converged = (it > itmin .and. &
all([ err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol, &
err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol, &
err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal))
if (Converged) then
reason = 1
elseif (it > itmax) then
reason = -1
else
reason = 0
endif
write(6,'(a,es14.7)') 'error stress BC = ', &
err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)
write(6,'(a,es14.7)') 'error F = ',&
err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol
write(6,'(a,es14.7)') 'error P = ', &
err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol
write(6,'(/,a)') '=========================================================================='
end subroutine AL_converged
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
subroutine AL_destroy()
use DAMASK_spectral_Utilities, only: &
Utilities_destroy
implicit none
PetscErrorCode :: ierr
call VecDestroy(solution_vec,ierr)
call SNESDestroy(snes,ierr)
call DMDestroy(da,ierr)
call PetscFinalize(ierr)
call Utilities_destroy()
subroutine AL_destroy()
end subroutine AL_destroy
use DAMASK_spectral_Utilities, only: &
Utilities_destroy
implicit none
PetscErrorCode :: ierr
end module DAMASK_spectral_SolverAL
call VecDestroy(solution_vec,ierr)
CHKERRQ(ierr)
call SNESDestroy(snes,ierr)
CHKERRQ(ierr)
call DMDestroy(da,ierr)
CHKERRQ(ierr)
call PetscFinalize(ierr)
CHKERRQ(ierr)
call Utilities_destroy()
end subroutine AL_destroy
end module DAMASK_spectral_SolverAL

156
code/DAMASK_spectral_solverBasic.f90
View File

@ -12,10 +12,8 @@ module DAMASK_spectral_SolverBasic
use prec, only: &
pInt, &
pReal
use math, only: &
math_I3
use DAMASK_spectral_Utilities, only: &
tSolutionState
@ -24,19 +22,23 @@ module DAMASK_spectral_SolverBasic
DAMASK_spectral_SolverBasic_label = 'basic'
!--------------------------------------------------------------------------------------------------
! pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, F_lastInc, Fdot, P
! pointwise global data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F, & !< deformation gradient field
F_lastInc, & !< deformation gradient field last increment
Fdot !< assumed rate for F n to F n+1
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private, dimension(:,:,:), allocatable :: temperature
real(pReal), private :: temperature !< not pointwise
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: &
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
F_aimDot = 0.0_pReal
F_aim = math_I3, & !< deformation gradient aim
F_aim_lastInc = math_I3, & !< deformation gradient aim last increment
F_aimDot = 0.0_pReal !< assumed rate
real(pReal), private,dimension(3,3,3,3) :: &
C = 0.0_pReal, C_lastInc = 0.0_pReal
C = 0.0_pReal, & !< average stiffness
C_lastInc = 0.0_pReal !< average stiffness last increment
contains
@ -50,27 +52,27 @@ subroutine basic_init()
IO_read_JobBinaryFile, &
IO_write_JobBinaryFile, &
IO_intOut
use FEsolving, only: &
restartInc
use DAMASK_interface, only: &
getSolverJobName
use DAMASK_spectral_Utilities, only: &
Utilities_init, &
Utilities_constitutiveResponse, &
Utilities_updateGamma, &
debugRestart
use mesh, only: &
res, &
geomdim
implicit none
integer(pInt) :: i,j,k
real(pReal), dimension(3,3) :: temp33_Real
real(pReal), dimension(3,3,3,3) :: temp3333_Real
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P
integer(pInt) :: &
i, j, k
real(pReal), dimension(3,3) :: &
temp33_Real
real(pReal), dimension(3,3,3,3) :: &
temp3333_Real
call Utilities_Init()
@ -79,15 +81,15 @@ subroutine basic_init()
#include "compilation_info.f90"
write(6,'(a)') ''
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lastInc ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (Fdot ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (P ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3),source = 0.0_pReal)
allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! init fields
! init fields and average quantities
if (restartInc == 1_pInt) then ! no deformation (no restart)
F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F_lastInc = F
@ -127,13 +129,9 @@ subroutine basic_init()
close (777)
coordinates = 0.0 ! change it later!!!
endif
call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,&
P,C,temp33_Real,.false.,math_I3)
call Utilities_constitutiveResponse(coordinates,F,F,temperature,0.0_pReal,P,C,temp33_Real,.false.,math_I3) ! constitutive response with no deformation in no time to get reference stiffness
!no rotation bc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
!--------------------------------------------------------------------------------------------------
! reference stiffness
if (restartInc == 1_pInt) then ! use initial stiffness as reference stiffness
temp3333_Real = C
endif
@ -148,7 +146,6 @@ end subroutine basic_init
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function &
basic_solution(incInfo,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
use numerics, only: &
itmax, &
itmin, &
@ -162,11 +159,10 @@ type(tSolutionState) function &
res,&
geomdim, &
deformed_fft, &
wgt
wgt
use IO, only: &
IO_write_JobBinaryFile, &
IO_intOut
use DAMASK_spectral_Utilities, only: &
tBoundaryCondition, &
field_real, &
@ -179,7 +175,6 @@ type(tSolutionState) function &
Utilities_updateGamma, &
Utilities_constitutiveResponse, &
Utilities_calculateRate
use FEsolving, only: &
restartWrite, &
restartRead, &
@ -190,26 +185,40 @@ type(tSolutionState) function &
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc
logical, intent(in) :: guess
type(tBoundaryCondition), intent(in) :: P_BC,F_BC
character(len=*), intent(in) :: incInfo
real(pReal), dimension(3,3), intent(in) :: rotation_BC
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
temperature_bc !< temperature (I wonder, why it is intent(in)
logical, intent(in) :: &
guess !< if .false., assume homogeneous addon
type(tBoundaryCondition), intent(in) :: &
P_BC,& !< stress boundary conditions
F_BC !< deformation boundary conditions
character(len=*), intent(in) :: &
incInfo !< string with information of current increment for output to screen
real(pReal), dimension(3,3), intent(in) :: &
rotation_BC !< rotation load to lab
real(pReal), dimension(3,3,3,3) :: S
real(pReal), dimension(3,3) :: F_aim_lab, &
F_aim_lab_lastIter, &
P_av
real(pReal), dimension(3,3,3,3) :: &
S !< current average compliance
real(pReal), dimension(3,3) :: &
F_aim_lab, &
F_aim_lab_lastIter, & !< aim of last iteration
P_av
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
real(pReal) :: err_div, err_stress
integer(pInt) :: iter, row, column, i, j, k
logical :: ForwardData
real(pReal) :: defgradDet, defgradDetMax, defgradDetMin
real(pReal) :: &
defgradDet, &
defgradDetMax, &
defgradDetMin
real(pReal), dimension(3,3) :: temp33_Real
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
! write restart information for spectral solver
if (restartWrite) then
write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F)) ! writing deformation gradient field to file
@ -235,14 +244,14 @@ type(tSolutionState) function &
close(777)
endif
if ( cutBack) then
!--------------------------------------------------------------------------------------------------
! forward data
if (cutBack) then
F_aim = F_aim_lastInc
F = F_lastInc
C = C_lastInc
else
C_lastInc = C
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
C_lastInc = C
if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F
f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim)
elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed
@ -251,9 +260,7 @@ type(tSolutionState) function &
if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), &
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & ! update coordinates and rate and forward last inc
1.0_pReal,F_lastInc,coordinates)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc,timeinc_old,guess,F_lastInc,F)
@ -265,13 +272,13 @@ type(tSolutionState) function &
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C)
if (update_gamma) call Utilities_updateGamma(C,restartWrite)
!--------------------------------------------------------------------------------------------------
! iteration till converged
if (.not. restartRead) ForwardData = .True.
iter = 0_pInt
convergenceLoop: do while(iter < itmax)
iter = iter + 1_pInt
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
@ -279,16 +286,16 @@ type(tSolutionState) function &
' @ 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
F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC)
basic_solution%termIll = .false.
call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,&
P,C,P_av,ForwardData,rotation_BC)
basic_solution%termIll = terminallyIll
terminallyIll = .False.
ForwardData = .False.
terminallyIll = .false.
ForwardData = .false.
!--------------------------------------------------------------------------------------------------
! stress BC handling
@ -326,7 +333,6 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress)
err_div_tol, &
err_stress_tolrel, &
err_stress_tolabs
use math, only: &
math_mul33x33, &
math_eigenvalues33, &
@ -344,8 +350,6 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress)
real(pReal) :: &
err_stress_tol, &
pAvgDivL2
pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(pAvgDiv,math_transpose33(pAvgDiv))))) ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html)
err_stress_tol = min(maxval(abs(pAvgStress))*err_stress_tolrel,err_stress_tolabs)
@ -363,8 +367,8 @@ end function basic_Converged
subroutine basic_destroy()
use DAMASK_spectral_Utilities, only: &
Utilities_destroy
use DAMASK_spectral_Utilities, only: &
Utilities_destroy
implicit none
call Utilities_destroy()
@ -372,44 +376,4 @@ subroutine basic_destroy()
end subroutine basic_destroy
end module DAMASK_spectral_SolverBasic
!--------------------------------------------------------------------------------------------------
! 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
!--------------------------------------------------------------------------------------------------
! 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

546
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -37,36 +37,37 @@ module DAMASK_spectral_SolverBasicPETSc
!--------------------------------------------------------------------------------------------------
! PETSc data
DM, private :: da
SNES, private :: snes
Vec, private :: solution_vec
DM, private :: da
SNES, private :: snes
Vec, private :: solution_vec
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private, dimension(:,:,:), allocatable :: temperature
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal) :: temperature
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: &
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
P_av
character(len=1024), private :: incInfo
real(pReal), private, dimension(3,3,3,3) :: &
C = 0.0_pReal, C_lastInc= 0.0_pReal, &
S = 0.0_pReal
real(pReal), private :: err_stress, err_div
logical, private :: ForwardData
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
real(pReal), private, dimension(3,3) :: &
F_aim = math_I3, &
F_aim_lastInc = math_I3, &
P_av
character(len=1024), private :: incInfo
real(pReal), private, dimension(3,3,3,3) :: &
C = 0.0_pReal, C_lastInc= 0.0_pReal, &
S = 0.0_pReal
real(pReal), private :: err_stress, err_div
logical, private :: ForwardData
real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
public :: basicPETSc_init, &
basicPETSc_solution ,&
basicPETSc_destroy
contains
public :: basicPETSc_init, &
basicPETSc_solution ,&
basicPETSc_destroy
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!--------------------------------------------------------------------------------------------------
@ -111,17 +112,19 @@ subroutine basicPETSc_init()
PetscObject :: dummy
call Utilities_init()
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSc init -+>>>'
write(6,'(a)') ' $Id: DAMASK_spectral_SolverBasicPETSC.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $'
#include "compilation_info.f90"
write(6,'(a)') ''
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr)
CHKERRQ(ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
@ -140,133 +143,130 @@ subroutine basicPETSc_init()
call SNESSetFromOptions(snes,ierr)
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
CHKERRQ(ierr)
if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo
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'
flush(6)
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
!--------------------------------------------------------------------------------------------------
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr) ! get the data out of PETSc to work with
CHKERRQ(ierr)
if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo
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'
flush(6)
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)
call Utilities_constitutiveResponse(coordinates,&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
temperature,0.0_pReal,P,C,P_av,.false.,math_I3)
coordinates = 0.0 ! change it later!!!
endif
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr)
call Utilities_constitutiveResponse(coordinates,&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
reshape(F(0:8,0:res(1)-1_pInt,0:res(2)-1_pInt,0:res(3)-1_pInt),[3,3,res(1),res(2),res(3)]),&
temperature,0.0_pReal,P,C,P_av,.false.,math_I3)
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr) ! write data back into PETSc
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! reference stiffness and Gamma operator
if (restartInc == 1_pInt) then
call IO_write_jobBinaryFile(777,'C_ref',size(C))
write (777,rec=1) C
close(777)
elseif (restartInc > 1_pInt) then
call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C))
read (777,rec=1) C
close (777)
endif
call Utilities_updateGamma(C,.True.)
!--------------------------------------------------------------------------------------------------
! reference stiffness
if (restartInc == 1_pInt) then
call IO_write_jobBinaryFile(777,'C_ref',size(C))
write (777,rec=1) C
close(777)
elseif (restartInc > 1_pInt) then
call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C))
read (777,rec=1) C
close (777)
endif
call Utilities_updateGamma(C,.True.)
end subroutine basicPETSc_init
end subroutine basicPETSc_init
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Basic PETSC scheme with internal iterations
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function &
basicPETSc_solution(incInfoIn,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
use numerics, only: &
update_gamma
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use mesh, only: &
res,&
geomdim,&
deformed_fft
use IO, only: &
IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: &
tBoundaryCondition, &
Utilities_calculateRate, &
Utilities_forwardField, &
Utilities_maskedCompliance, &
Utilities_updateGamma, &
cutBack
use FEsolving, only: &
restartWrite, &
terminallyIll
implicit none
type(tSolutionState) function &
basicPETSc_solution(incInfoIn,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC)
use numerics, only: &
update_gamma
use math, only: &
math_mul33x33 ,&
math_rotate_backward33
use mesh, only: &
res,&
geomdim,&
deformed_fft
use IO, only: &
IO_write_JobBinaryFile
use DAMASK_spectral_Utilities, only: &
tBoundaryCondition, &
Utilities_calculateRate, &
Utilities_forwardField, &
Utilities_maskedCompliance, &
Utilities_updateGamma, &
cutBack
use FEsolving, only: &
restartWrite, &
terminallyIll
implicit none
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90>
!--------------------------------------------------------------------------------------------------
! input data for solution
real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc
logical, intent(in):: guess
type(tBoundaryCondition), intent(in) :: P_BC,F_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC
real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc
logical, intent(in):: guess
type(tBoundaryCondition), intent(in) :: P_BC,F_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC
character(len=*), intent(in) :: incInfoIn
real(pReal), dimension(3,3),save :: F_aimDot=0.0_pReal
real(pReal), dimension(3,3) :: &
F_aim_lab
real(pReal), dimension(3,3),save :: F_aimDot=0.0_pReal
real(pReal), dimension(3,3) :: &
F_aim_lab
!--------------------------------------------------------------------------------------------------
! loop variables, convergence etc.
real(pReal), dimension(3,3) :: temp33_Real
real(pReal), dimension(3,3) :: temp33_Real
!--------------------------------------------------------------------------------------------------
!
PetscScalar, pointer :: F(:,:,:,:)
PetscErrorCode :: ierr
SNESConvergedReason :: reason
PetscScalar, pointer :: F(:,:,:,:)
PetscErrorCode :: ierr
SNESConvergedReason :: reason
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
incInfo = incInfoIn
if (restartWrite) then
write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc))
write (777,rec=1) F_LastInc
close (777)
call IO_write_jobBinaryFile(777,'C',size(C))
write (777,rec=1) C
close(777)
endif
call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
incInfo = incInfoIn
if (restartWrite) then
write(6,'(a)') 'writing converged results for restart'
call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc))
write (777,rec=1) F_LastInc
close (777)
call IO_write_jobBinaryFile(777,'C',size(C))
write (777,rec=1) C
close(777)
endif
call DMDAVecGetArrayF90(da,solution_vec,F,ierr)
if ( cutBack) then
F_aim = F_aim_lastInc
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
C = C_lastInc
else
C_lastInc = C
if ( cutBack) then
F_aim = F_aim_lastInc
F = reshape(F_lastInc,[9,res(1),res(2),res(3)])
C = C_lastInc
else
C_lastInc = C
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
@ -277,7 +277,7 @@ else
endif
if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old
F_aim_lastInc = F_aim
!--------------------------------------------------------------------------------------------------
! update coordinates and rate and forward last inc
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), &
@ -291,167 +291,177 @@ else
F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr)
CHKERRQ(ierr)
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C)
if (update_gamma) call Utilities_updateGamma(C,restartWrite)
ForwardData = .True.
mask_stress = P_BC%maskFloat
params%P_BC = P_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C)
if (update_gamma) call Utilities_updateGamma(C,restartWrite)
ForwardData = .True.
mask_stress = P_BC%maskFloat
params%P_BC = P_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr)
call SNESGetConvergedReason(snes,reason,ierr)
basicPETSc_solution%termIll = terminallyIll
terminallyIll = .false.
BasicPETSC_solution%converged =.false.
if (reason > 0 ) then
BasicPETSC_solution%converged = .true.
endif
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr)
CHKERRQ(ierr)
call SNESGetConvergedReason(snes,reason,ierr)
CHKERRQ(ierr)
basicPETSc_solution%termIll = terminallyIll
terminallyIll = .false.
BasicPETSC_solution%converged =.false.
if (reason > 0 ) then
BasicPETSC_solution%converged = .true.
endif
end function BasicPETSc_solution
end function BasicPETSc_solution
!--------------------------------------------------------------------------------------------------
!> @brief forms the AL residual vector
!--------------------------------------------------------------------------------------------------
subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr)
use numerics, only: &
itmax, &
itmin
use math, only: &
math_rotate_backward33, &
math_transpose33, &
math_mul3333xx33
use mesh, only: &
res, &
wgt
use DAMASK_spectral_Utilities, only: &
field_real, &
Utilities_FFTforward, &
Utilities_FFTbackward, &
Utilities_fourierConvolution, &
Utilities_constitutiveResponse, &
Utilities_divergenceRMS
use IO, only : IO_intOut
implicit none
subroutine BasicPETSC_formResidual(myIn,x_scal,f_scal,dummy,ierr)
real(pReal), dimension(3,3) :: F_aim_lab_lastIter, F_aim_lab
DMDALocalInfo, dimension(*) :: myIn
PetscScalar, dimension(3,3,res(1),res(2),res(3)) :: x_scal
PetscScalar, dimension(3,3,res(1),res(2),res(3)):: f_scal
PetscInt :: iter, nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr)
call SNESGetIterationNumber(snes,iter,ierr)
use numerics, only: &
itmax, &
itmin
use math, only: &
math_rotate_backward33, &
math_transpose33, &
math_mul3333xx33
use mesh, only: &
res, &
wgt
use DAMASK_spectral_Utilities, only: &
field_real, &
Utilities_FFTforward, &
Utilities_FFTbackward, &
Utilities_fourierConvolution, &
Utilities_constitutiveResponse, &
Utilities_divergenceRMS
use IO, only : IO_intOut
implicit none
real(pReal), dimension(3,3) :: F_aim_lab_lastIter, F_aim_lab
DMDALocalInfo, dimension(*) :: myIn
PetscScalar, dimension(3,3,res(1),res(2),res(3)) :: x_scal
PetscScalar, dimension(3,3,res(1),res(2),res(3)):: f_scal
PetscInt :: iter, nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr)
CHKERRQ(ierr)
call SNESGetIterationNumber(snes,iter,ierr)
CHKERRQ(ierr)
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ 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,params%rotation_BC)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(coordinates,F_lastInc,x_scal,temperature,params%timeinc, &
f_scal,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .false.
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), &
' @ 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,params%rotation_BC)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(coordinates,F_lastInc,x_scal,temperature,params%timeinc, &
f_scal,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .false.
!--------------------------------------------------------------------------------------------------
! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc
err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc
F_aim_lab = math_rotate_backward33(F_aim,params%rotation_BC)
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%P_BC))) ! S = 0.0 for no bc
err_stress = maxval(abs(mask_stress * (P_av - params%P_BC))) ! mask = 0.0 for no bc
F_aim_lab = math_rotate_backward33(F_aim,params%rotation_BC)
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
field_real = 0.0_pReal
field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(f_scal,[res(1),res(2),res(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
err_div = Utilities_divergenceRMS()
call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab)
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2])
write(6,'(/,a)') '=========================================================================='
end subroutine BasicPETSc_formResidual
field_real = 0.0_pReal
field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(f_scal,[res(1),res(2),res(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_FFTforward()
err_div = Utilities_divergenceRMS()
call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab)
call Utilities_FFTbackward()
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(x_scal),order=[3,4,5,1,2])
write(6,'(/,a)') '=========================================================================='
end subroutine BasicPETSc_formResidual
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
err_div_tol, &
err_stress_tolrel, &
err_stress_tolabs
use math, only: &
math_mul33x33, &
math_eigenvalues33, &
math_transpose33
implicit none
subroutine BasicPETSc_converged(snes_local,it,xnorm,snorm,fnorm,reason,dummy,ierr)
SNES :: snes_local
PetscInt :: it
PetscReal :: xnorm, snorm, fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
logical :: Converged
real(pReal) :: pAvgDivL2
pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,math_transpose33(P_av)))))
Converged = (it > itmin .and. &
all([ err_div/pAvgDivL2/err_div_tol, &
err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal))
if (Converged) then
reason = 1
elseif (it > itmax) then
reason = -1
else
reason = 0
endif
write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/pAvgDivL2/err_div_tol,&
' (',err_div/pAvgDivL2,' N/m³)'
write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs), &
' (',err_stress,' Pa)'
use numerics, only: &
itmax, &
itmin, &
err_div_tol, &
err_stress_tolrel, &
err_stress_tolabs
use math, only: &
math_mul33x33, &
math_eigenvalues33, &
math_transpose33
implicit none
SNES :: snes_local
PetscInt :: it
PetscReal :: xnorm, snorm, fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
logical :: Converged
real(pReal) :: pAvgDivL2
pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,math_transpose33(P_av)))))
Converged = (it > itmin .and. &
all([ err_div/pAvgDivL2/err_div_tol, &
err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal))
if (Converged) then
reason = 1
elseif (it > itmax) then
reason = -1
else
reason = 0
endif
write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/pAvgDivL2/err_div_tol,&
' (',err_div/pAvgDivL2,' N/m³)'
write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs), &
' (',err_stress,' Pa)'
end subroutine BasicPETSc_converged
end subroutine BasicPETSc_converged
!--------------------------------------------------------------------------------------------------
!> @brief destroy routine
!--------------------------------------------------------------------------------------------------
subroutine BasicPETSc_destroy()
subroutine BasicPETSc_destroy()
use DAMASK_spectral_Utilities, only: &
Utilities_destroy
implicit none
PetscErrorCode :: ierr
use DAMASK_spectral_Utilities, only: &
Utilities_destroy
implicit none
PetscErrorCode :: ierr
call VecDestroy(solution_vec,ierr)
call SNESDestroy(snes,ierr)
call DMDestroy(da,ierr)
call PetscFinalize(ierr)
call Utilities_destroy()
call VecDestroy(solution_vec,ierr)
CHKERRQ(ierr)
call SNESDestroy(snes,ierr)
CHKERRQ(ierr)
call DMDestroy(da,ierr)
CHKERRQ(ierr)
call PetscFinalize(ierr)
CHKERRQ(ierr)
call Utilities_destroy()
CHKERRQ(ierr)
end subroutine BasicPETSc_destroy
end subroutine BasicPETSc_destroy
end module DAMASK_spectral_SolverBasicPETSc
end module DAMASK_spectral_SolverBasicPETSc

83
code/DAMASK_spectral_utilities.f90
View File

@ -17,26 +17,29 @@ module DAMASK_spectral_utilities
!--------------------------------------------------------------------------------------------------
! variables storing information for spectral method and FFTW
real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress of deformation) of field_fourier
type(C_PTR), private :: plan_forward, plan_backward !< plans for FFTW
real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !< real representation (some stress or deformation) of field_fourier
complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_fourier !< field on which the Fourier transform operates
real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method
real(pReal), private, dimension(:,:,:,:), allocatable :: xi !< wave vector field for divergence and for gamma operator
complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_fourier !< field on which the Fourier transform operates
real(pReal), private, dimension(3,3,3,3) :: C_ref !< reference stiffness
!--------------------------------------------------------------------------------------------------
! debug fftw
type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back !< plans for FFTW in case of debugging the Fourier transform
complex(pReal),private, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for debug of FFTW
complex(pReal),private, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field complex representation for debug of FFTW
complex(pReal),private, dimension(:,:,:), pointer :: scalarField_real, & !< scalar field real representation for debug of FFTW
scalarField_fourier !< scalar field complex representation for debug of FFTW
!--------------------------------------------------------------------------------------------------
! debug divergence
type(C_PTR), private :: plan_divergence !< plan for FFTW in case of debugging divergence calculation
real(pReal), private, dimension(:,:,:,:), pointer :: divergence_real !< scalar field real representation for debugging divergence calculation
complex(pReal),private, dimension(:,:,:,:), pointer :: divergence_fourier !< scalar field real representation for debugging divergence calculation
real(pReal), private, dimension(:,:,:,:), allocatable :: divergence_post !< data of divergence calculation using function from core modules (serves as a reference)
!--------------------------------------------------------------------------------------------------
! plans for FFTW
type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back !< plans for FFTW in case of debugging the Fourier transform
type(C_PTR), private :: plan_forward, plan_backward !< plans for FFTW
type(C_PTR), private :: plan_divergence !< plan for FFTW in case of debugging divergence calculation
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical, public :: &
@ -80,8 +83,8 @@ module DAMASK_spectral_utilities
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, sets debug flags, create plans for fftw
!> @details Sets the debug levels for general, divergence, restart and fftw from the biwise coding
!> @brief allocates all neccessary fields, sets debug flags, create plans for FFTW
!> @details Sets the debug levels for general, divergence, restart and FFTW from the biwise coding
!> provided by the debug module to logicals.
!> Allocates all fields used by FFTW and create the corresponding plans depending on the debug
!> level chosen.
@ -132,8 +135,8 @@ subroutine utilities_init()
!--------------------------------------------------------------------------------------------------
! allocation
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)
allocate (xi(3,res1_red,res(2),res(3)),source = 0.0_pReal) ! frequencies, only half the size for first dimension
tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate aligned 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
@ -148,7 +151,7 @@ subroutine utilities_init()
call fftw_set_timelimit(fftw_timelimit) ! set timelimit for plan creation
!--------------------------------------------------------------------------------------------------
! creating plans
! creating plans for the convolution
plan_forward = fftw_plan_many_dft_r2c(3, [res(3),res(2) ,res(1)], 9,& ! dimensions, logical length in each dimension in reversed order, no. of transforms
field_real, [res(3),res(2) ,res(1)+2_pInt],& ! input data, physical length in each dimension in reversed order
1, res(3)*res(2)*(res(1)+2_pInt),& ! striding, product of physical length in the 3 dimensions
@ -162,7 +165,7 @@ subroutine utilities_init()
1, res(3)*res(2)*(res(1)+2_pInt), fftw_planner_flag) ! striding, product of physical length in the 3 dimensions, planner precision
!--------------------------------------------------------------------------------------------------
! depending on (debug) options, allocate more memory and create additional plans
! depending on debug options, allocate more memory and create additional plans
if (debugDivergence) then
divergence = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T))
call c_f_pointer(divergence, divergence_real, [ res(1)+2_pInt,res(2),res(3),3])
@ -267,7 +270,7 @@ end subroutine utilities_updateGamma
!--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed
!> Does an unweighted FFT transform from real to complex.
!> @detailed Does an unweighted FFT transform from real to complex.
!> In case of debugging the FFT, also one component of the tensor (specified by row and column)
!> is independetly transformed complex to complex and compared to the whole tensor transform
!--------------------------------------------------------------------------------------------------
@ -326,7 +329,7 @@ end subroutine utilities_FFTforward
!--------------------------------------------------------------------------------------------------
!> @brief backward FFT of data in field_fourier to field_real
!> Does an inverse FFT transform from complex to real
!> @detailed Does an inverse FFT transform from complex to real
!> In case of debugging the FFT, also one component of the tensor (specified by row and column)
!> is independetly transformed complex to complex and compared to the whole tensor transform
!> results is weighted by number of points stored in wgt
@ -343,7 +346,7 @@ subroutine utilities_FFTbackward(row,column)
integer(pInt) :: i, j, k, m, n
!--------------------------------------------------------------------------------------------------
! unpack FFT data for conj complex symmetric part
! unpack FFT data for conj complex symmetric part. This data is not transformed when using c2r
if (debugFFTW) then
scalarField_fourier = field_fourier(1:res1_red,1:res(2),1:res(3),row,column)
do i = 0_pInt, res(1)/2_pInt-2_pInt
@ -377,6 +380,29 @@ subroutine utilities_FFTbackward(row,column)
field_real = field_real * wgt ! normalize the result by number of elements
!--------------------------------------------------------------------------------------------------
! 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
end subroutine utilities_FFTbackward
@ -397,7 +423,7 @@ subroutine utilities_fourierConvolution(fieldAim)
real(pReal), intent(in), dimension(3,3) :: fieldAim !< desired average value of the field after convolution
real(pReal), dimension(3,3) :: xiDyad, temp33_Real
real(pReal) :: filter !< weighting of current component
complex(pReal), dimension(3,3) :: temp33_complex
complex(pReal), dimension(3,3) :: temp33_complex
integer(pInt) :: &
i, j, k, &
l, m, n, o
@ -621,7 +647,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
CPFEM_general
implicit none
real(pReal), intent(inout), dimension(res(1),res(2),res(3)) :: temperature !< temperature field
real(pReal), intent(inout) :: temperature !< temperature (no field)
real(pReal), intent(in), dimension(res(1),res(2),res(3),3) :: coordinates !< coordinates field
real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: &
F_lastInc, & !< target deformation gradient
@ -655,7 +681,20 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
calcMode = 2_pInt
collectMode = 5_pInt
endif
!--------------------------------------------------------------------------------------------------
! 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
if (DebugGeneral) write(6,*) 'collect mode: ', collectMode,' calc mode: ', calcMode
flush(6)
@ -664,7 +703,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
ielem = ielem + 1_pInt
call CPFEM_general(collectMode,& ! collect cycle
coordinates(i,j,k,1:3), F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), &
temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF)
temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF)
collectMode = 3_pInt
enddo; enddo; enddo
@ -676,7 +715,7 @@ subroutine utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
ielem = ielem + 1_pInt
call CPFEM_general(calcMode,& ! first element in first iteration retains CPFEM_mode 1,
coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort)
temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF)
temperature,timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF)
calcMode = 2_pInt
C = C + dPdF
enddo; enddo; enddo
@ -711,7 +750,7 @@ pure function utilities_calculateRate(delta_aim,timeinc,timeinc_old,guess,field_
field !< data of current step
real(pReal), dimension(3,3,res(1),res(2),res(3)) :: utilities_calculateRate
if (guess) then
if(guess) then
utilities_calculateRate = (field-field_lastInc) / timeinc_old
else
utilities_calculateRate = spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3))/timeinc