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

changed order of most arrays to fortran-fast, whole new solver still work in progress

This commit is contained in:
Martin Diehl 2012-08-09 13:04:56 +00:00
parent fe4d4d9525
commit b324a98014
4 changed files with 150 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

163
code/DAMASK_spectral_SolverAL.f90
View File

@ -31,6 +31,7 @@ module DAMASK_spectral_SolverAL
#include <finclude/petscpc.h> #include <finclude/petscpc.h>
#include <finclude/petscsnes.h> #include <finclude/petscsnes.h>
#include <finclude/petscvec.h90> #include <finclude/petscvec.h90>
#include <finclude/petscdmda.h90>
#include <finclude/petscsnes.h90> #include <finclude/petscsnes.h90>
#endif #endif
@ -48,12 +49,13 @@ module DAMASK_spectral_SolverAL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM, private :: da
SNES, private :: snes SNES, private :: snes
Vec, private :: solution_vec Vec, private :: solution_vec
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! common pointwise data ! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, F_lastInc, F_lambda, F_lambda_lastInc, P real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc, P
real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates
real(pReal), private, dimension(:,:,:), allocatable :: temperature real(pReal), private, dimension(:,:,:), allocatable :: temperature
@ -72,8 +74,7 @@ module DAMASK_spectral_SolverAL
real(pReal), private :: err_stress, err_f, err_p real(pReal), private :: err_stress, err_f, err_p
logical, private :: ForwardData logical, private :: ForwardData
real(pReal), private, dimension(3,3) :: & real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal
mask_stress = 0.0_pReal
contains contains
@ -103,7 +104,8 @@ module DAMASK_spectral_SolverAL
use mesh, only: & use mesh, only: &
res, & res, &
geomdim geomdim, &
mesh_NcpElems
use math, only: & use math, only: &
math_invSym3333 math_invSym3333
@ -115,7 +117,7 @@ module DAMASK_spectral_SolverAL
PetscErrorCode :: ierr_psc PetscErrorCode :: ierr_psc
PetscObject :: dummy PetscObject :: dummy
PetscMPIInt :: rank PetscMPIInt :: rank
DM :: da PetscScalar, pointer :: xx_psc(:,:,:,:)
call Utilities_init() call Utilities_init()
@ -125,39 +127,61 @@ module DAMASK_spectral_SolverAL
#include "compilation_info.f90" #include "compilation_info.f90"
write(6,'(a)') '' write(6,'(a)') ''
allocate (F ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lambda ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (P (3,3, res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lambda_lastInc(res(1),res(2),res(3),3,3), source = 0.0_pReal)
allocate (P ( res(1), res(2),res(3),3,3), source = 0.0_pReal)
allocate (coordinates( res(1), res(2),res(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) allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! PETSc Init
call PetscInitialize(PETSC_NULL_CHARACTER,ierr_psc)
call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr_psc)
call SNESCreate(PETSC_COMM_WORLD,snes,ierr_psc)
call DMDACreate3d(PETSC_COMM_WORLD, &
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr_psc)
call DMCreateGlobalVector(da,solution_vec,ierr_psc)
call DMDASetLocalFunction(da,AL_formResidual,ierr_psc)
call SNESSetDM(snes,da,ierr_psc)
call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr_psc)
call PetscOptionsInsertString(petsc_options,ierr_psc)
call SNESSetFromOptions(snes,ierr_psc)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc)
if (restartInc == 1_pInt) then ! no deformation (no restart) if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F_lambda_lastInc = F_lastInc
xx_psc(0:8,:,:,:) = reshape(F_lastInc,[9,res(1),res(2),res(3)])
xx_psc(9:17,:,:,:) = xx_psc(0:8,:,:,:)
call flush(6)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
F(i,j,k,1:3,1:3) = math_I3
F_lastInc(i,j,k,1:3,1:3) = math_I3
F_lambda(i,j,k,1:3,1:3) = math_I3
F_lambda_lastInc(i,j,k,1:3,1:3) = math_I3
coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*res,pReal) - geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo enddo; enddo; enddo
print*, 'init done12'
elseif (restartInc > 1_pInt) then ! using old values from file elseif (restartInc > 1_pInt) then ! using old values from file
if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',&
restartInc - 1_pInt,' from file' restartInc - 1_pInt,' from file'
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',&
trim(getSolverJobName()),size(F)) trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F read (777,rec=1) xx_psc(0:8,:,:,:)
close (777) close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',&
trim(getSolverJobName()),size(F_lastInc)) trim(getSolverJobName()),size(F_lastInc))
read (777,rec=1) F_lastInc read (777,rec=1) F_lastInc
close (777) close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',& call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',&
trim(getSolverJobName()),size(F_lambda)) trim(getSolverJobName()),size(F_lambda_lastInc))
read (777,rec=1) F read (777,rec=1) xx_psc(9:17,:,:,:)
close (777) close (777)
call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',& call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',&
trim(getSolverJobName()),size(F_lambda_lastInc)) trim(getSolverJobName()),size(F_lambda_lastInc))
@ -173,9 +197,12 @@ module DAMASK_spectral_SolverAL
coordinates = 0.0 ! change it later!!! coordinates = 0.0 ! change it later!!!
endif endif
call Utilities_constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& call Utilities_constitutiveResponse(coordinates,reshape(xx_psc(0:8,:,:,:),shape(F_lastInc)),&
P,C,P_av,.false.,math_I3) reshape(xx_psc(0:8,:,:,:),shape(F_lastInc)),&
temperature,0.0_pReal,P,C,P_av,.false.,math_I3)
print*, 'const response'
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc)
print*, 'restored'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reference stiffness ! reference stiffness
if (restartInc == 1_pInt) then if (restartInc == 1_pInt) then
@ -192,25 +219,6 @@ module DAMASK_spectral_SolverAL
C_scale = C C_scale = C
S_scale = math_invSym3333(C) S_scale = math_invSym3333(C)
!--------------------------------------------------------------------------------------------------
! PETSc Init
call PetscInitialize(PETSC_NULL_CHARACTER,ierr_psc)
call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr_psc)
call SNESCreate(PETSC_COMM_WORLD,snes,ierr_psc)
call DMDACreate3d(PETSC_COMM_WORLD, &
DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, &
DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, &
18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr_psc)
call DMCreateGlobalVector(da,solution_vec,ierr_psc)
call DMDASetLocalFunction(da,AL_formResidual,ierr_psc)
call SNESSetDM(snes,da,ierr_psc)
call SNESSetConvergenceTest(snes,AL_converged,dummy,PETSC_NULL_FUNCTION,ierr_psc)
call PetscOptionsInsertString(petsc_options,ierr_psc)
call SNESSetFromOptions(snes,ierr_psc)
call DMDestroy(da,ierr_psc)
end subroutine AL_init end subroutine AL_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -255,7 +263,7 @@ module DAMASK_spectral_SolverAL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! !
PetscScalar, pointer :: xx_psc(:) PetscScalar, pointer :: xx_psc(:,:,:,:)
PetscErrorCode ierr_psc PetscErrorCode ierr_psc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -287,8 +295,12 @@ module DAMASK_spectral_SolverAL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update local deformation gradient and coordinates ! update local deformation gradient and coordinates
deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC)
call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc,F) call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc)
call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lambda_lastInc,F_lambda) call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc, &
xx_psc(1:9,1:res(1),1:res(2),1:res(3)))
!call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lambda_lastInc,&
! reshape(xx_psc(10:18,1:res(1),1:res(2),1:res(3)),shape(F_lambda_lastInc)))
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc)
call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -302,14 +314,6 @@ module DAMASK_spectral_SolverAL
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = timeinc params%timeinc = timeinc
ctr = 1_pInt
call VecGetArrayF90(solution_vec,xx_psc,ierr_psc)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
xx_psc(ctr:ctr+8_pInt) = reshape(F(i,j,k,1:3,1:3),[9])
xx_psc(ctr+9_pInt:ctr+17_pInt) = reshape(F_lambda(i,j,k,1:3,1:3),[9])
ctr = ctr + 18_pInt
enddo; enddo; enddo
call VecRestoreArrayF90(solution_vec,xx_psc,ierr_psc)
call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr_psc) call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr_psc)
call SNESGetConvergedReason(snes,reason,ierr_psc) call SNESGetConvergedReason(snes,reason,ierr_psc)
if (reason > 0 ) AL_solution%converged = .true. if (reason > 0 ) AL_solution%converged = .true.
@ -340,7 +344,7 @@ module DAMASK_spectral_SolverAL
implicit none implicit none
integer(pInt) :: i,j,k integer(pInt) :: i,j,k,l, ctr
real(pReal), dimension (3,3) :: temp33_real real(pReal), dimension (3,3) :: temp33_real
DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE) DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE)
@ -361,16 +365,10 @@ module DAMASK_spectral_SolverAL
write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',&
math_transpose33(F_aim) math_transpose33(F_aim)
do k=in(DMDA_LOCAL_INFO_ZS)+1,in(DMDA_LOCAL_INFO_ZS)+in(DMDA_LOCAL_INFO_ZM)
do j=in(DMDA_LOCAL_INFO_YS)+1,in(DMDA_LOCAL_INFO_YS)+in(DMDA_LOCAL_INFO_YM)
do i=in(DMDA_LOCAL_INFO_XS)+1,in(DMDA_LOCAL_INFO_XS)+in(DMDA_LOCAL_INFO_XM)
F(i,j,k,1:3,1:3) = reshape(x_scal(1:9,i,j,k),[3,3])
F_lambda(i,j,k,1:3,1:3) = reshape(x_scal(10:18,i,j,k),[3,3])
enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc,& call Utilities_constitutiveResponse(coordinates,F_lastInc,reshape(x_scal(1:9,1:res(1),1:res(2),1:res(3)),shape(F_lastInc)),&
temperature,params%timeinc,&
P,C,P_av,ForwardData,params%rotation_BC) P,C,P_av,ForwardData,params%rotation_BC)
ForwardData = .False. ForwardData = .False.
@ -383,37 +381,38 @@ module DAMASK_spectral_SolverAL
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing Fourier transform ! doing Fourier transform
field_real = 0.0_pReal field_real = 0.0_pReal
do j = 1_pInt, 3_pInt; do i = 1_pInt, 3_pInt
ctr = 1_pInt
do k = 1_pInt, 3_pInt; do l = 1_pInt, 3_pInt
field_real(1:res(1),1:res(2),1:res(3),i,j) = field_real(1:res(1),1:res(2),1:res(3),i,j) + &
C_scale(i,j,l,k)*(x_scal(9+ctr,1:res(1),1:res(2),1:res(3)) - &
x_scal(ctr,1:res(1),1:res(2),1:res(3)))
! P_temp(i,j,1:res(1),1:res(2),1:res(3)) = P_temp(i,j,1:res(1),1:res(2),1:res(3)) + &
! S_scale(i,j,l,k)*P(l,k,1:res(1),1:res(2),1:res(3))
ctr = ctr + 1_pInt
enddo; enddo
enddo; enddo
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_scale,F_lambda(i,j,k,1:3,1:3)-F(i,j,k,1:3,1:3)) P(1:3,1:3,i,j,k) = math_mul3333xx33(S_scale,P(1:3,1:3,i,j,k)) + math_I3
enddo; enddo; enddo enddo; enddo; enddo
call Utilities_forwardFFT() call Utilities_forwardFFT()
call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC)) call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC))
call Utilities_backwardFFT() call Utilities_backwardFFT()
err_f = 0.0_pReal f_scal(1:9,1:res(1),1:res(2),1:res(3)) = reshape(P,[9,res(1),res(2),res(3)]) - x_scal(10:18,1:res(1),1:res(2),1:res(3)) +&
err_p = 0.0_pReal x_scal(1:9,1:res(1),1:res(2),1:res(3)) - &
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),[9,res(1),res(2),res(3)],order=[2,3,4,1])
temp33_real = field_real(i,j,k,1:3,1:3) - F(i,j,k,1:3,1:3) f_scal(10:18,1:res(1),1:res(2),1:res(3)) = x_scal(1:9,1:res(1),1:res(2),1:res(3)) - &
err_f = err_f + sum(temp33_real*temp33_real) reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),&
[9,res(1),res(2),res(3)],order=[2,3,4,1])
temp33_real = F_lambda(i,j,k,1:3,1:3) - & err_f = wgt*sqrt(sum(f_scal(10:18,1:res(1),1:res(2),1:res(3))**2.0_pReal))
(math_mul3333xx33(S_scale,P(i,j,k,1:3,1:3)) + math_I3) err_p = wgt*sqrt(sum((f_scal( 1:9 ,1:res(1),1:res(2),1:res(3)) &
err_p = err_p + sum(temp33_real*temp33_real) -f_scal(10:18,1:res(1),1:res(2),1:res(3)))**2.0_pReal))
enddo; enddo; enddo
err_f = wgt*sqrt(err_f)
err_p = wgt*sqrt(err_p)
do k=in(DMDA_LOCAL_INFO_ZS)+1,in(DMDA_LOCAL_INFO_ZS)+in(DMDA_LOCAL_INFO_ZM)
do j=in(DMDA_LOCAL_INFO_YS)+1,in(DMDA_LOCAL_INFO_YS)+in(DMDA_LOCAL_INFO_YM)
do i=in(DMDA_LOCAL_INFO_XS)+1,in(DMDA_LOCAL_INFO_XS)+in(DMDA_LOCAL_INFO_XM)
temp33_real = math_mul3333xx33(S_scale,P(i,j,k,1:3,1:3)) + math_I3 - F_lambda(i,j,k,1:3,1:3) &
+ F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3)
f_scal(1:9,i,j,k) = reshape(temp33_real,[9])
f_scal(10:18,i,j,k) = reshape(F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3),[9])
enddo; enddo; enddo
return
end subroutine AL_formResidual end subroutine AL_formResidual
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

21
code/DAMASK_spectral_SolverBasic.f90
View File

@ -79,18 +79,18 @@ subroutine basic_init()
write(6,'(a)') '' write(6,'(a)') ''
allocate (F ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) allocate (F_lastInc ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal)
allocate (P ( res(1), res(2),res(3),3,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 (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal)
allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
if (restartInc == 1_pInt) then ! no deformation (no restart) if (restartInc == 1_pInt) then ! no deformation (no restart)
F = spread(spread(spread(math_I3,3,res(1)),4,res(2)),5,res(3)) ! initialize to identity
F_lastInc = F
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) 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) & coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) &
- geomdim/real(2_pInt*res,pReal) - geomdim/real(2_pInt*res,pReal)
enddo; enddo; enddo enddo; enddo; enddo
@ -114,7 +114,7 @@ subroutine basic_init()
coordinates = 0.0 ! change it later!!! coordinates = 0.0 ! change it later!!!
endif endif
!no rotation bc call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates)
call Utilities_constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& call Utilities_constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,&
P,C,temp33_Real,.false.,math_I3) P,C,temp33_Real,.false.,math_I3)
@ -259,17 +259,14 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme ! updated deformation gradient using fix point algorithm of basic scheme
field_real = 0.0_pReal field_real = 0.0_pReal
field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = P field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(P,[res(1),res(2),res(3),3,3],&
order=[4,5,1,2,3]) ! field real has a different order
call Utilities_forwardFFT() call Utilities_forwardFFT()
err_div = Utilities_divergenceRMS() err_div = Utilities_divergenceRMS()
call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab) call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab)
call Utilities_backwardFFT() call Utilities_backwardFFT()
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) F = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization
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
basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av) basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av)
if (basic_solution%converged .and. iter > itmin) exit if (basic_solution%converged .and. iter > itmin) exit
enddo convergenceLoop enddo convergenceLoop

98
code/DAMASK_spectral_Utilities.f90
View File

@ -71,6 +71,11 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, sets debug flags, create plans for fftw !> @brief allocates all neccessary fields, sets debug flags, create plans for fftw
!> @details Sets the debug levels for general, divergence, restart and fftw from the biwise coding
!> provided by the debug module to logicals.
!> Allocates all fields used by FFTW and create the corresponding plans depending on the debug
!> level chosen.
!> Initializes FFTW.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine Utilities_init() subroutine Utilities_init()
@ -185,15 +190,19 @@ subroutine Utilities_init()
enddo; enddo; enddo enddo; enddo; enddo
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(3,3,3,3,1,1,1), source = 0.0_pReal)
else ! precalculation of gamma_hat field else ! precalculation of gamma_hat field
allocate (gamma_hat(res1_red ,res(2),res(3),3,3,3,3), source =0.0_pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 allocate (gamma_hat(3,3,3,3,res1_red ,res(2),res(3)), source =0.0_pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
endif endif
end subroutine Utilities_init end subroutine Utilities_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief updates references stiffness and potentially precalculated gamma operator !> @brief updates references stiffness and potentially precalculated gamma operator
!> @details Sets the current reference stiffness to the stiffness given as an argument.
!> If the gamma operator is precalculated, it is calculated with this stiffness.
!> In case of a on-the-fly calculation, only the reference stiffness is updated.
!> The gamma operator is filtered depening on the filter selected in numerics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine Utilities_updateGamma(C) subroutine Utilities_updateGamma(C)
@ -218,15 +227,18 @@ subroutine Utilities_updateGamma(C)
temp33_Real = math_inv33(temp33_Real) temp33_Real = math_inv33(temp33_Real)
filter = Utilities_getFilter(xi(1:3,i,j,k)) filter = Utilities_getFilter(xi(1:3,i,j,k))
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)& forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)&
gamma_hat(i,j,k, l,m,n,o) = filter*temp33_Real(l,n)*xiDyad(m,o) gamma_hat(l,m,n,o, i,j,k) = filter*temp33_Real(l,n)*xiDyad(m,o)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 gamma_hat(1:3,1:3,1:3,1:3, 1,1,1) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1
endif endif
end subroutine Utilities_updateGamma end subroutine Utilities_updateGamma
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed !> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed
!> Does an unweighted FFT transform from real to complex.
!> In case of debugging the FFT, also one component of the tensor (specified by row and column)
!> is independetly transformed complex to complex and compared to the whole tensor transform
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine Utilities_forwardFFT(row,column) subroutine Utilities_forwardFFT(row,column)
use mesh, only : & use mesh, only : &
@ -276,6 +288,14 @@ subroutine Utilities_forwardFFT(row,column)
= cmplx(0.0_pReal,0.0_pReal,pReal) = cmplx(0.0_pReal,0.0_pReal,pReal)
end subroutine Utilities_forwardFFT end subroutine Utilities_forwardFFT
!--------------------------------------------------------------------------------------------------
!> @brief backward FFT of data in field_fourier to field_real
!> Does an inverse FFT transform from complex to real
!> In case of debugging the FFT, also one component of the tensor (specified by row and column)
!> is independetly transformed complex to complex and compared to the whole tensor transform
!> results is weighted by number of points stored in wgt
!--------------------------------------------------------------------------------------------------
subroutine Utilities_backwardFFT(row,column) subroutine Utilities_backwardFFT(row,column)
implicit none implicit none
@ -286,9 +306,7 @@ subroutine Utilities_backwardFFT(row,column)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 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 scalarField_fourier = field_fourier(1:res1_red,1:res(2),1:res(3),row,column)
scalarField_fourier(i,j,k) = field_fourier(i,j,k,row,column)
enddo; enddo; enddo
do i = 0_pInt, res(1)/2_pInt-2_pInt ! unpack fft data for conj complex symmetric part 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)
@ -303,8 +321,7 @@ subroutine Utilities_backwardFFT(row,column)
enddo; enddo enddo; enddo
endif endif
!--------------------------------------------------------------------------------------------------
! doing the inverse FT
call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -353,16 +370,16 @@ subroutine Utilities_fourierConvolution(fieldAim)
temp33_Real = math_inv33(temp33_Real) temp33_Real = math_inv33(temp33_Real)
filter = Utilities_getFilter(xi(1:3,i,j,k)) filter = Utilities_getFilter(xi(1:3,i,j,k))
forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)& forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)&
gamma_hat(1,1,1, l,m,n,o) = filter*temp33_Real(l,n)*xiDyad(m,o) gamma_hat(l,m,n,o, 1,1,1) = filter*temp33_Real(l,n)*xiDyad(m,o)
forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) &
temp33_Complex(l,m) = sum(gamma_hat(1,1,1, l,m, 1:3,1:3) * field_fourier(i,j,k,1:3,1:3)) temp33_Complex(l,m) = sum(gamma_hat(l,m,1:3,1:3, 1,1,1) * field_fourier(i,j,k,1:3,1:3))
field_fourier(i,j,k,1:3,1:3) = temp33_Complex field_fourier(i,j,k,1:3,1:3) = temp33_Complex
endif endif
enddo; enddo; enddo enddo; enddo; enddo
else ! use precalculated gamma-operator else ! use precalculated gamma-operator
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
forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt) & forall( m = 1_pInt:3_pInt, n = 1_pInt:3_pInt) &
temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) * field_fourier(i,j,k,1:3,1:3)) temp33_Complex(m,n) = sum(gamma_hat(m,n,1:3,1:3, i,j,k) * field_fourier(i,j,k,1:3,1:3))
field_fourier(i,j,k, 1:3,1:3) = temp33_Complex field_fourier(i,j,k, 1:3,1:3) = temp33_Complex
enddo; enddo; enddo enddo; enddo; enddo
endif endif
@ -422,18 +439,9 @@ real(pReal) function Utilities_divergenceRMS()
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 = sqrt(wgt*sum(divergence_real**2.0_pReal)) ! RMS in real space
err_post_div_RMS = 0.0_pReal err_post_div_RMS = sqrt(wgt*sum(divergence_post**2.0_pReal)) ! RMS in real space
err_real_div_max = 0.0_pReal err_real_div_max = sqrt(maxval(sum(divergence_real**2.0_pReal,dim=4))) ! max in real space
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
err_real_div_RMS = err_real_div_RMS + sum(divergence_real(i,j,k,1:3)**2.0_pReal) ! avg of squared L_2 norm of div(stress) in real space
err_post_div_RMS = err_post_div_RMS + sum(divergence_post(i,j,k,1:3)**2.0_pReal) ! avg of squared L_2 norm of div(stress) in real space
err_real_div_max = max(err_real_div_max,sum(divergence_real(i,j,k,1:3)**2.0_pReal)) ! max of squared L_2 norm of div(stress) in real space
enddo; enddo; enddo
err_real_div_RMS = sqrt(wgt*err_real_div_RMS) ! RMS in real space
err_post_div_RMS = sqrt(wgt*err_post_div_RMS) ! RMS in real space
err_real_div_max = sqrt( err_real_div_max) ! max in real space
err_div_max = sqrt( err_div_max) ! max in Fourier space err_div_max = sqrt( err_div_max) ! max in Fourier space
write(6,'(1x,a,es11.4)') 'error divergence FT RMS = ',err_div_RMS write(6,'(1x,a,es11.4)') 'error divergence FT RMS = ',err_div_RMS
@ -516,7 +524,7 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
real(pReal), dimension(res(1),res(2),res(3)) :: temperature real(pReal), dimension(res(1),res(2),res(3)) :: temperature
real(pReal), dimension(res(1),res(2),res(3),3) :: coordinates real(pReal), dimension(res(1),res(2),res(3),3) :: coordinates
real(pReal), dimension(res(1),res(2),res(3),3,3) :: F,F_lastInc, P real(pReal), dimension(3,3,res(1),res(2),res(3)) :: F,F_lastInc, P
real(pReal) :: timeinc real(pReal) :: timeinc
logical :: ForwardData logical :: ForwardData
integer(pInt) :: i, j, k, ielem integer(pInt) :: i, j, k, ielem
@ -524,7 +532,7 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
real(pReal), dimension(3,3,3,3) :: dPdF, C real(pReal), dimension(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, P_av, rotation_BC real(pReal), dimension(3,3) :: P_av, rotation_BC
write(6,'(a)') '' write(6,'(a)') ''
write(6,'(a)') '... evaluating constitutive response .................' write(6,'(a)') '... evaluating constitutive response .................'
@ -539,28 +547,27 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) 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(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(i,j,k,1:3,1:3),dPdF) temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),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
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(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(i,j,k,1:3,1:3),dPdF) temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),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)
enddo; enddo; enddo enddo; enddo; enddo
call debug_info() call debug_info()
P_av = math_rotate_forward33(sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt,rotation_BC) !average of P rotated
restartWrite = .false. restartWrite = .false.
P_av_lab = P_av_lab * wgt
P_av = math_rotate_forward33(P_av_lab,rotation_BC)
write (6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') ' Piola-Kirchhoff stress / MPa =',& write (6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') ' Piola-Kirchhoff stress / MPa =',&
math_transpose33(P_av)/1.e6_pReal math_transpose33(P_av)/1.e6_pReal
@ -574,17 +581,16 @@ subroutine Utilities_forwardField(delta_aim,timeinc,timeinc_old,guessmode,field_
real(pReal), intent(in), dimension(3,3) :: delta_aim real(pReal), intent(in), dimension(3,3) :: delta_aim
real(pReal), intent(in) :: timeinc, timeinc_old, guessmode real(pReal), intent(in) :: timeinc, timeinc_old, guessmode
real(pReal), intent(inout), dimension(res(1),res(2),res(3),3,3) :: field_lastInc,field real(pReal), intent(inout), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,field
integer(pInt) :: i,j,k
real(pReal), dimension(3,3) :: temp33_real if (guessmode == 1.0_pReal) then
field = field + (field-field_lastInc) * timeinc/timeinc_old
field_lastInc = (field + field_lastInc * timeinc/timeinc_old) /(1.0_pReal + timeinc/timeinc_old)
else
field_lastInc = field
field = field + spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3))
endif
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
temp33_Real = Field(i,j,k,1:3,1:3)
Field(i,j,k,1:3,1:3) = Field(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (field(i,j,k,1:3,1:3) - Field_lastInc(i,j,k,1:3,1:3))*timeinc/timeinc_old& ! guessing...
+ (1.0_pReal-guessmode) * delta_aim ! if not guessing, use prescribed average deformation where applicable
Field_lastInc(i,j,k,1:3,1:3) = temp33_Real
enddo; enddo; enddo
end subroutine Utilities_forwardField end subroutine Utilities_forwardField
real(pReal) function Utilities_getFilter(k) real(pReal) function Utilities_getFilter(k)
@ -614,9 +620,7 @@ subroutine Utilities_destroy()
implicit none implicit none
if (debugDivergence) then if (debugDivergence) call fftw_destroy_plan(plan_divergence)
call fftw_destroy_plan(plan_divergence)
endif
if (debugFFTW) then if (debugFFTW) then
call fftw_destroy_plan(plan_scalarField_forth) call fftw_destroy_plan(plan_scalarField_forth)

2
code/DAMASK_spectral_interface.f90
View File

@ -17,7 +17,7 @@
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>. ! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
! !
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!* $Id$ ! $Id$
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, 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 !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH