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

clearer naming, debug options for spectral do not work for MPI

This commit is contained in:
Martin Diehl 2015-09-11 08:52:03 +00:00
parent 751d1d7582
commit e88cedc6ae
8 changed files with 546 additions and 699 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
code/DAMASK_spectral_driver.f90
View File

@ -43,8 +43,8 @@ program DAMASK_spectral_Driver
debug_levelBasic
use math ! need to include the whole module for FFTW
use mesh, only: &
gridGlobal, &
geomSizeGlobal
grid, &
geomSize
use CPFEM, only: &
CPFEM_initAll
use FEsolving, only: &
@ -390,8 +390,8 @@ program DAMASK_spectral_Driver
write(resUnit) 'load:', trim(loadCaseFile) ! ... and write header
write(resUnit) 'workingdir:', trim(getSolverWorkingDirectoryName())
write(resUnit) 'geometry:', trim(geometryFile)
write(resUnit) 'grid:', gridGlobal
write(resUnit) 'size:', geomSizeGlobal
write(resUnit) 'grid:', grid
write(resUnit) 'size:', geomSize
write(resUnit) 'materialpoint_sizeResults:', materialpoint_sizeResults
write(resUnit) 'loadcases:', size(loadCases)
write(resUnit) 'frequencies:', loadCases%outputfrequency ! one entry per LoadCase

100
code/DAMASK_spectral_solverAL.f90
View File

@ -12,7 +12,7 @@ module DAMASK_spectral_solverAL
pReal
use math, only: &
math_I3
use DAMASK_spectral_utilities, only: &
use DAMASK_spectral_Utilities, only: &
tSolutionState, &
tSolutionParams
@ -121,8 +121,8 @@ subroutine AL_init
Utilities_updateGamma, &
Utilities_updateIPcoords
use mesh, only: &
gridLocal, &
gridGlobal
grid, &
grid3
use math, only: &
math_invSym3333
@ -145,36 +145,35 @@ subroutine AL_init
#include "compilation_info.f90"
endif
allocate (P (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (P (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (Fdot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_lambda_lastInc(3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_lambdaDot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (F_lambda_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (F_lambdaDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! PETSc Init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
do proc = 1, worldsize
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid
grid(1),grid(2),grid(3), & ! global grid
1 , 1, worldsize, &
18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
gridLocal (1),gridLocal (2),localK, & ! local grid
grid(1),grid(2),localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
call DMCreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr)
call DMDASNESSetFunctionLocal(da,INSERT_VALUES,AL_formResidual,dummy,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)
@ -184,13 +183,7 @@ subroutine AL_init
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(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,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) ! initialize to identity
F = reshape(F_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F_lambda = F
F_lambda_lastInc = F_lastInc
elseif (restartInc > 1_pInt) then ! read in F to calculate coordinates and initialize CPFEM general
restart: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment ', restartInc - 1_pInt, ' from file'
@ -218,16 +211,22 @@ subroutine AL_init
call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
read (777,rec=1) f_aimDot
close (777)
endif
elseif (restartInc == 1_pInt) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
F_lambda = F
F_lambda_lastInc = F_lastInc
endif restart
call utilities_updateIPcoords(F)
call Utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(F_lastInc,F,&
0.0_pReal,P,C_volAvg,C_minMaxAvg,temp33_Real,.false.,math_I3)
nullify(F)
nullify(F_lambda)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc
if (restartInc > 1_pInt) then ! using old values from files
readRestart: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -241,7 +240,7 @@ subroutine AL_init
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg
close (777)
endif
endif readRestart
call Utilities_updateGamma(C_minMaxAvg,.True.)
C_scale = C_minMaxAvg
@ -339,7 +338,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
polarBeta, &
worldrank
use mesh, only: &
gridLocal
grid3, &
grid
use IO, only: &
IO_intOut
use math, only: &
@ -350,7 +350,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
math_mul33x33
use DAMASK_spectral_Utilities, only: &
wgt, &
tensorField_realMPI, &
tensorField_real, &
utilities_FFTtensorForward, &
utilities_fourierGammaConvolution, &
Utilities_inverseLaplace, &
@ -408,7 +408,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
if(totalIter <= PETScIter) then ! new iteration
newIteration: if(totalIter <= PETScIter) then
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
totalIter = totalIter + 1_pInt
@ -420,15 +420,15 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim)
endif
flush(6)
endif
endif newIteration
!--------------------------------------------------------------------------------------------------
!
tensorField_realMPI = 0.0_pReal
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
tensorField_realMPI(1:3,1:3,i,j,k) = &
tensorField_real = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
tensorField_real(1:3,1:3,i,j,k) = &
polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_lambda(1:3,1:3,i,j,k) - math_I3))
@ -443,7 +443,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_lambda = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
residual_F_lambda = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
@ -455,8 +455,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! calculate divergence
tensorField_realMPI = 0.0_pReal
tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F
call utilities_FFTtensorForward()
err_div = Utilities_divergenceRMS()
call utilities_FFTtensorBackward()
@ -464,7 +464,7 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0_pInt
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
e = e + 1_pInt
residual_F(1:3,1:3,i,j,k) = math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - &
@ -475,8 +475,8 @@ subroutine AL_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! calculating curl
tensorField_realMPI = 0.0_pReal
tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward()
err_curl = Utilities_curlRMS()
call utilities_FFTtensorBackward()
@ -570,7 +570,8 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
use numerics, only: &
worldrank
use mesh, only: &
gridLocal
grid3, &
grid
use DAMASK_spectral_Utilities, only: &
Utilities_calculateRate, &
Utilities_forwardField, &
@ -639,14 +640,14 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
write (777,rec=1) C_volAvgLastInc
close(777)
endif
endif
call utilities_updateIPcoords(F)
if (cutBack) then
F_aim = F_aim_lastInc
F_lambda = reshape(F_lambda_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F = reshape(F_lastInc, [9,gridLocal(1),gridLocal(2),gridLocal(3)])
F_lambda = reshape(F_lambda_lastInc,[9,grid(1),grid(2),grid3])
F = reshape(F_lastInc, [9,grid(1),grid(2),grid3])
C_volAvg = C_volAvgLastInc
else
ForwardData = .True.
@ -667,22 +668,25 @@ subroutine AL_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,rotation_
! update coordinates and rate and forward last inc
call utilities_updateIPcoords(F)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]))
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]))
F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,gridLocal(1), &
gridLocal(2),gridLocal(3)]))
F_lastInc = reshape(F, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
F_lambda_lastInc = reshape(F_lambda,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
timeinc_old,guess,F_lambda_lastInc,reshape(F_lambda,[3,3,grid(1), &
grid(2),grid3]))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
F_lambda_lastInc = reshape(F_lambda,[3,3,grid(1),grid(2),grid3])
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
! update local deformation gradient
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim
math_rotate_backward33(F_aim,rotation_BC)), &
[9,gridLocal(1),gridLocal(2),gridLocal(3)])
[9,grid(1),grid(2),grid3])
F_lambda = reshape(Utilities_forwardField(timeinc,F_lambda_lastInc,F_lambdadot), &
[9,gridLocal(1),gridLocal(2),gridLocal(3)]) ! does not have any average value as boundary condition
[9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition
if (.not. guess) then ! large strain forwarding
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
F_lambda33 = reshape(F_lambda(1:9,i,j,k),[3,3])
F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
math_mul3333xx33(C_scale,&

59
code/DAMASK_spectral_solverBasicPETSc.f90
View File

@ -110,8 +110,8 @@ subroutine basicPETSc_init
utilities_updateIPcoords, &
wgt
use mesh, only: &
gridLocal, &
gridGlobal
grid, &
grid3
use math, only: &
math_invSym3333
@ -133,27 +133,27 @@ subroutine basicPETSc_init
#include "compilation_info.f90"
endif mainProcess
allocate (P (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (P (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (Fdot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
do proc = 1, worldsize
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid
grid(1),grid(2),grid(3), & ! global grid
1, 1, worldsize, &
9, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
gridLocal (1),gridLocal (2),localK, & ! local grid
grid (1),grid (2),localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
@ -169,10 +169,7 @@ subroutine basicPETSc_init
! init fields
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! get the data out of PETSc to work with
if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) ! initialize to identity
F = reshape(F_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
elseif (restartInc > 1_pInt) then ! using old values from file
restart: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment ', restartInc - 1_pInt, ' from file'
@ -189,7 +186,10 @@ subroutine basicPETSc_init
close (777)
F_aim = reshape(sum(sum(sum(F,dim=4),dim=3),dim=2) * wgt, [3,3]) ! average of F
F_aim_lastInc = sum(sum(sum(F_lastInc,dim=5),dim=4),dim=3) * wgt ! average of F_lastInc
endif
elseif (restartInc == 1_pInt) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restart
call Utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(F_lastInc, F, &
@ -202,7 +202,7 @@ subroutine basicPETSc_init
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! write data back to PETSc
if (restartInc > 1_pInt) then ! using old values from files
restartRead: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -216,7 +216,7 @@ subroutine basicPETSc_init
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg
close (777)
endif
endif restartRead
call Utilities_updateGamma(C_minmaxAvg,.True.)
@ -302,7 +302,8 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: &
worldrank
use mesh, only: &
gridLocal
grid, &
grid3
use math, only: &
math_rotate_backward33, &
math_transpose33, &
@ -312,9 +313,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
debug_spectral, &
debug_spectralRotation
use DAMASK_spectral_Utilities, only: &
wgt, &
tensorField_realMPI, &
tensorField_fourierMPI, &
tensorField_real, &
utilities_FFTtensorForward, &
utilities_FFTtensorBackward, &
utilities_fourierGammaConvolution, &
@ -345,7 +344,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
if (totalIter <= PETScIter) then ! new iteration
newIteration: if (totalIter <= PETScIter) then
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
totalIter = totalIter + 1_pInt
@ -357,9 +356,9 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim)
endif
flush(6)
endif
endif newIteration
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
@ -376,8 +375,8 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
tensorField_realMPI = 0.0_pReal
tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = f_scal
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = f_scal
call utilities_FFTtensorForward()
err_div = Utilities_divergenceRMS()
call utilities_fourierGammaConvolution(math_rotate_backward33(F_aim_lastIter-F_aim,params%rotation_BC))
@ -385,7 +384,7 @@ subroutine BasicPETSC_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
f_scal = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
end subroutine BasicPETSc_formResidual
@ -455,7 +454,8 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
math_mul33x33 ,&
math_rotate_backward33
use mesh, only: &
gridLocal
grid, &
grid3
use DAMASK_spectral_Utilities, only: &
Utilities_calculateRate, &
Utilities_forwardField, &
@ -513,9 +513,10 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
endif
call utilities_updateIPcoords(F)
if (cutBack) then
F_aim = F_aim_lastInc
F = reshape(F_lastInc, [9,gridLocal(1),gridLocal(2),gridLocal(3)])
F = reshape(F_lastInc, [9,grid(1),grid(2),grid3])
C_volAvg = C_volAvgLastInc
else
ForwardData = .True.
@ -536,15 +537,15 @@ subroutine BasicPETSc_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC,r
! update coordinates and rate and forward last inc
call utilities_updateIPcoords(F)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc,reshape(F,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]))
F_lastInc = reshape(F, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
timeinc_old,guess,F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
! update local deformation gradient
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim
math_rotate_backward33(F_aim,rotation_BC)),[9,gridLocal(1),gridLocal(2),gridLocal(3)])
math_rotate_backward33(F_aim,rotation_BC)),[9,grid(1),grid(2),grid3])
call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
end subroutine BasicPETSc_forward

104
code/DAMASK_spectral_solverPolarisation.f90
View File

@ -12,7 +12,7 @@ module DAMASK_spectral_solverPolarisation
pReal
use math, only: &
math_I3
use DAMASK_spectral_utilities, only: &
use DAMASK_spectral_Utilities, only: &
tSolutionState, &
tSolutionParams
@ -105,7 +105,7 @@ subroutine Polarisation_init
IO_intOut, &
IO_read_realFile, &
IO_timeStamp
use debug, only : &
use debug, only: &
debug_level, &
debug_spectral, &
debug_spectralRestart
@ -121,8 +121,8 @@ subroutine Polarisation_init
Utilities_updateGamma, &
Utilities_updateIPcoords
use mesh, only: &
gridLocal, &
gridGlobal
grid, &
grid3
use math, only: &
math_invSym3333
@ -145,29 +145,29 @@ subroutine Polarisation_init
#include "compilation_info.f90"
endif
allocate (P (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (P (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! allocate global fields
allocate (F_lastInc (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (Fdot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_tau_lastInc(3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_tauDot (3,3,gridLocal(1),gridLocal(2),gridLocal(3)),source = 0.0_pReal)
allocate (F_lastInc (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (Fdot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (F_tau_lastInc(3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
allocate (F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!--------------------------------------------------------------------------------------------------
! PETSc Init
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
do proc = 1, worldsize
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid
grid(1),grid(2),grid(3), & ! global grid
1 , 1, worldsize, &
18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
gridLocal (1),gridLocal (2),localK, & ! local grid
grid (1),grid (2),localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr)
@ -183,12 +183,7 @@ subroutine Polarisation_init
call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! places pointer xx_psc on PETSc data
F => xx_psc(0:8,:,:,:)
F_tau => xx_psc(9:17,:,:,:)
if (restartInc == 1_pInt) then ! no deformation (no restart)
F_lastInc = spread(spread(spread(math_I3,3,gridLocal(1)),4,gridLocal(2)),5,gridLocal(3)) ! initialize to identity
F = reshape(F_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F_tau = 2.0_pReal* F
F_tau_lastInc = 2.0_pReal*F_lastInc
elseif (restartInc > 1_pInt) then
restart: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading values of increment', restartInc - 1_pInt, 'from file'
@ -216,7 +211,12 @@ subroutine Polarisation_init
call IO_read_realFile(777,'F_aimDot',trim(getSolverJobName()),size(f_aimDot))
read (777,rec=1) f_aimDot
close (777)
endif
elseif (restartInc == 1_pInt) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
F_tau = 2.0_pReal* F
F_tau_lastInc = 2.0_pReal*F_lastInc
endif restart
call Utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(F_lastInc,F,&
@ -225,7 +225,7 @@ subroutine Polarisation_init
nullify(F_tau)
call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr); CHKERRQ(ierr) ! write data back to PETSc
if (restartInc > 1_pInt) then ! using old values from files
readRestart: if (restartInc > 1_pInt) then
if (iand(debug_level(debug_spectral),debug_spectralRestart)/= 0 .and. worldrank == 0_pInt) &
write(6,'(/,a,'//IO_intOut(restartInc-1_pInt)//',a)') &
'reading more values of increment', restartInc - 1_pInt, 'from file'
@ -239,7 +239,7 @@ subroutine Polarisation_init
call IO_read_realFile(777,'C_ref',trim(getSolverJobName()),size(C_minMaxAvg))
read (777,rec=1) C_minMaxAvg
close (777)
endif
endif readRestart
call Utilities_updateGamma(C_minMaxAvg,.True.)
C_scale = C_minMaxAvg
@ -335,10 +335,11 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
polarAlpha, &
polarBeta, &
worldrank
use mesh, only: &
grid3, &
grid
use IO, only: &
IO_intOut
use mesh, only: &
gridLocal
use math, only: &
math_rotate_backward33, &
math_transpose33, &
@ -347,7 +348,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
math_mul33x33
use DAMASK_spectral_Utilities, only: &
wgt, &
tensorField_realMPI, &
tensorField_real, &
utilities_FFTtensorForward, &
utilities_fourierGammaConvolution, &
Utilities_inverseLaplace, &
@ -405,7 +406,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if(nfuncs== 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
if (totalIter <= PETScIter) then ! new iteration
newIteration: if(totalIter <= PETScIter) then
!--------------------------------------------------------------------------------------------------
! report begin of new iteration
totalIter = totalIter + 1_pInt
@ -417,15 +418,15 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
math_transpose33(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' deformation gradient aim =', &
math_transpose33(F_aim)
endif
flush(6)
endif
endif newIteration
!--------------------------------------------------------------------------------------------------
!
tensorField_realMPI = 0.0_pReal
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
tensorField_realMPI(1:3,1:3,i,j,k) = &
tensorField_real = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
tensorField_real(1:3,1:3,i,j,k) = &
polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
@ -439,7 +440,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
residual_F_tau = polarBeta*F - tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
@ -451,8 +452,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! calculate divergence
tensorField_realMPI = 0.0_pReal
tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = residual_F
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F
call utilities_FFTtensorForward()
err_div = Utilities_divergenceRMS()
call utilities_FFTtensorBackward()
@ -460,7 +461,7 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0_pInt
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
e = e + 1_pInt
residual_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
@ -471,8 +472,8 @@ subroutine Polarisation_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! calculating curl
tensorField_realMPI = 0.0_pReal
tensorField_realMPI(1:3,1:3,1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = F
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward()
err_curl = Utilities_curlRMS()
call utilities_FFTtensorBackward()
@ -491,8 +492,8 @@ subroutine Polarisation_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,
err_div_tolAbs, &
err_curl_tolRel, &
err_curl_tolAbs, &
err_stress_tolabs, &
err_stress_tolrel, &
err_stress_tolAbs, &
err_stress_tolRel, &
worldrank
use math, only: &
math_mul3333xx33
@ -563,20 +564,21 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
math_mul3333xx33, &
math_transpose33, &
math_rotate_backward33
use numerics, only: &
worldrank
use mesh, only: &
grid3, &
grid
use DAMASK_spectral_Utilities, only: &
Utilities_calculateRate, &
Utilities_forwardField, &
Utilities_updateIPcoords, &
tBoundaryCondition, &
cutBack
use mesh, only: &
gridLocal
use IO, only: &
IO_write_JobRealFile
use FEsolving, only: &
restartWrite
use numerics, only: &
worldrank
implicit none
real(pReal), intent(in) :: &
@ -636,11 +638,10 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
endif
call utilities_updateIPcoords(F)
if (cutBack) then
F_aim = F_aim_lastInc
F_tau= reshape(F_tau_lastInc,[9,gridLocal(1),gridLocal(2),gridLocal(3)])
F = reshape(F_lastInc, [9,gridLocal(1),gridLocal(2),gridLocal(3)])
F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid3])
F = reshape(F_lastInc, [9,grid(1),grid(2),grid3])
C_volAvg = C_volAvgLastInc
else
ForwardData = .True.
@ -662,24 +663,25 @@ subroutine Polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,F_BC,P_BC
call utilities_updateIPcoords(F)
Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), &
timeinc_old,guess,F_lastInc, &
reshape(F,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]))
reshape(F,[3,3,grid(1),grid(2),grid3]))
F_tauDot = Utilities_calculateRate(math_rotate_backward33(2.0_pReal*f_aimDot,rotation_BC), &
timeinc_old,guess,F_tau_lastInc, &
reshape(F_tau,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)]))
F_lastInc = reshape(F, [3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
F_tau_lastInc = reshape(F_tau,[3,3,gridLocal(1),gridLocal(2),gridLocal(3)])
reshape(F_tau,[3,3,grid(1),grid(2),grid3]))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
endif
F_aim = F_aim + f_aimDot * timeinc
!--------------------------------------------------------------------------------------------------
! update local deformation gradient
F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! ensure that it matches rotated F_aim
math_rotate_backward33(F_aim,rotation_BC)), &
[9,gridLocal(1),gridLocal(2),gridLocal(3)])
[9,grid(1),grid(2),grid3])
F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & ! does not have any average value as boundary condition
[9,gridLocal(1),gridLocal(2),gridLocal(3)])
[9,grid(1),grid(2),grid3])
if (.not. guess) then ! large strain forwarding
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt, gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt, grid(1)
F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, &
math_mul3333xx33(C_scale,&

665
code/DAMASK_spectral_utilities.f90
View File

File diff suppressed because it is too large Load Diff

112
code/mesh.f90
View File

@ -30,15 +30,17 @@ module mesh
mesh_Nelems !< total number of elements in mesh
#ifdef Spectral
integer(pInt), dimension(3), public, protected :: &
grid !< (global) grid
integer(pInt), public, protected :: &
mesh_NcpElemsGlobal, & !< total number of CP elements in global mesh
gridGlobal(3), &
gridLocal (3), &
gridOffset
grid3, & !< (local) grid in 3rd direction
grid3Offset !< (local) grid offset in 3rd direction
real(pReal), dimension(3), public, protected :: &
geomSize
real(pReal), public, protected :: &
geomSizeGlobal(3), &
geomSizeLocal (3), &
geomSizeOffset
size3, & !< (local) size in 3rd direction
size3offset !< (local) size offset in 3rd direction
#endif
integer(pInt), dimension(:,:), allocatable, public, protected :: &
@ -565,31 +567,21 @@ subroutine mesh_init(ip,el)
#ifdef Spectral
#ifdef PETSc
call fftw_mpi_init()
#endif
call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file...
if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6)
grid = mesh_spectral_getGrid(fileUnit)
geomSize = mesh_spectral_getSize(fileUnit)
gridGlobal = mesh_spectral_getGrid(fileUnit)
gridMPI = int(gridGlobal,C_INTPTR_T)
#ifdef PETSc
gridMPI = int(grid,C_INTPTR_T)
alloc_local = fftw_mpi_local_size_3d(gridMPI(3), gridMPI(2), gridMPI(1)/2 +1, &
MPI_COMM_WORLD, local_K, local_K_offset)
gridLocal = [gridGlobal(1:2),int(local_K,pInt)]
gridOffset = int(local_K_offset,pInt)
grid3 = int(local_K,pInt)
grid3Offset = int(local_K_offset,pInt)
geomSizeGlobal = mesh_spectral_getSize(fileUnit)
geomSizeLocal = [geomSizeGlobal(1:2), &
geomSizeGlobal(3)*real(gridLocal(3),pReal)/real(gridGlobal(3),pReal)]
geomSizeOffset = geomSizeGlobal(3)*real(gridOffset,pReal) /real(gridGlobal(3),pReal)
#else
call IO_open_file(FILEUNIT,geometryFile) ! parse info from geometry file...
if (myDebug) write(6,'(a)') ' Opened geometry file'; flush(6)
gridGlobal = mesh_spectral_getGrid(fileUnit)
gridLocal = gridGlobal
gridOffset = 0_pInt
geomSizeGlobal = mesh_spectral_getSize(fileUnit)
geomSizeLocal = geomSizeGlobal
geomSizeOffset = 0.0_pReal
size3 = geomSize(3)*real(grid3,pReal) /real(grid(3),pReal)
size3Offset = geomSize(3)*real(grid3Offset,pReal)/real(grid(3),pReal)
#endif
if (myDebug) write(6,'(a)') ' Grid partitioned'; flush(6)
call mesh_spectral_count()
@ -1256,11 +1248,11 @@ subroutine mesh_spectral_count()
implicit none
mesh_Nelems = product(gridLocal)
mesh_Nelems = product(grid(1:2))*grid3
mesh_NcpElems= mesh_Nelems
mesh_Nnodes = product(gridLocal + 1_pInt)
mesh_Nnodes = product(grid(1:2) + 1_pInt)*(grid3 + 1_pInt)
mesh_NcpElemsGlobal = product(gridGlobal)
mesh_NcpElemsGlobal = product(grid)
end subroutine mesh_spectral_count
@ -1319,15 +1311,15 @@ subroutine mesh_spectral_build_nodes()
forall (n = 0_pInt:mesh_Nnodes-1_pInt)
mesh_node0(1,n+1_pInt) = mesh_unitlength * &
geomSizeLocal(1)*real(mod(n,(gridLocal(1)+1_pInt) ),pReal) &
/ real(gridLocal(1),pReal)
geomSize(1)*real(mod(n,(grid(1)+1_pInt) ),pReal) &
/ real(grid(1),pReal)
mesh_node0(2,n+1_pInt) = mesh_unitlength * &
geomSizeLocal(2)*real(mod(n/(gridLocal(1)+1_pInt),(gridLocal(2)+1_pInt)),pReal) &
/ real(gridLocal(2),pReal)
geomSize(2)*real(mod(n/(grid(1)+1_pInt),(grid(2)+1_pInt)),pReal) &
/ real(grid(2),pReal)
mesh_node0(3,n+1_pInt) = mesh_unitlength * &
geomSizeLocal(3)*real(mod(n/(gridLocal(1)+1_pInt)/(gridLocal(2)+1_pInt),(gridLocal(3)+1_pInt)),pReal) &
/ real(gridLocal(3),pReal) + &
geomSizeOffset
size3*real(mod(n/(grid(1)+1_pInt)/(grid(2)+1_pInt),(grid3+1_pInt)),pReal) &
/ real(grid3,pReal) + &
size3offset
end forall
mesh_node = mesh_node0
@ -1422,7 +1414,7 @@ subroutine mesh_spectral_build_elements(fileUnit)
enddo
elemType = FE_mapElemtype('C3D8R')
elemOffset = gridLocal(1)*gridLocal(2)*gridOffset
elemOffset = product(grid(1:2))*grid3Offset
e = 0_pInt
do while (e < mesh_NcpElems) ! fill expected number of elements, stop at end of data (or blank line!)
e = e+1_pInt ! valid element entry
@ -1430,15 +1422,15 @@ subroutine mesh_spectral_build_elements(fileUnit)
mesh_element( 2,e) = elemType ! elem type
mesh_element( 3,e) = homog ! homogenization
mesh_element( 4,e) = mesh_microGlobal(e+elemOffset) ! microstructure
mesh_element( 5,e) = e + (e-1_pInt)/gridLocal(1) + &
((e-1_pInt)/(gridLocal(1)*gridLocal(2)))*(gridLocal(1)+1_pInt) ! base node
mesh_element( 5,e) = e + (e-1_pInt)/grid(1) + &
((e-1_pInt)/(grid(1)*grid(2)))*(grid(1)+1_pInt) ! base node
mesh_element( 6,e) = mesh_element(5,e) + 1_pInt
mesh_element( 7,e) = mesh_element(5,e) + gridLocal(1) + 2_pInt
mesh_element( 8,e) = mesh_element(5,e) + gridLocal(1) + 1_pInt
mesh_element( 9,e) = mesh_element(5,e) +(gridLocal(1) + 1_pInt) * (gridLocal(2) + 1_pInt) ! second floor base node
mesh_element( 7,e) = mesh_element(5,e) + grid(1) + 2_pInt
mesh_element( 8,e) = mesh_element(5,e) + grid(1) + 1_pInt
mesh_element( 9,e) = mesh_element(5,e) +(grid(1) + 1_pInt) * (grid(2) + 1_pInt) ! second floor base node
mesh_element(10,e) = mesh_element(9,e) + 1_pInt
mesh_element(11,e) = mesh_element(9,e) + gridLocal(1) + 2_pInt
mesh_element(12,e) = mesh_element(9,e) + gridLocal(1) + 1_pInt
mesh_element(11,e) = mesh_element(9,e) + grid(1) + 2_pInt
mesh_element(12,e) = mesh_element(9,e) + grid(1) + 1_pInt
mesh_maxValStateVar(1) = max(mesh_maxValStateVar(1),mesh_element(3,e)) ! needed for statistics
mesh_maxValStateVar(2) = max(mesh_maxValStateVar(2),mesh_element(4,e))
enddo
@ -1465,32 +1457,32 @@ subroutine mesh_spectral_build_ipNeighborhood(fileUnit)
allocate(mesh_ipNeighborhood(3,mesh_maxNipNeighbors,mesh_maxNips,mesh_NcpElems),source=0_pInt)
e = 0_pInt
do z = 0_pInt,gridLocal(3)-1_pInt
do y = 0_pInt,gridLocal(2)-1_pInt
do x = 0_pInt,gridLocal(1)-1_pInt
do z = 0_pInt,grid3-1_pInt
do y = 0_pInt,grid(2)-1_pInt
do x = 0_pInt,grid(1)-1_pInt
e = e + 1_pInt
mesh_ipNeighborhood(1,1,1,e) = z * gridLocal(1) * gridLocal(2) &
+ y * gridLocal(1) &
+ modulo(x+1_pInt,gridLocal(1)) &
mesh_ipNeighborhood(1,1,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x+1_pInt,grid(1)) &
+ 1_pInt
mesh_ipNeighborhood(1,2,1,e) = z * gridLocal(1) * gridLocal(2) &
+ y * gridLocal(1) &
+ modulo(x-1_pInt,gridLocal(1)) &
mesh_ipNeighborhood(1,2,1,e) = z * grid(1) * grid(2) &
+ y * grid(1) &
+ modulo(x-1_pInt,grid(1)) &
+ 1_pInt
mesh_ipNeighborhood(1,3,1,e) = z * gridLocal(1) * gridLocal(2) &
+ modulo(y+1_pInt,gridLocal(2)) * gridLocal(1) &
mesh_ipNeighborhood(1,3,1,e) = z * grid(1) * grid(2) &
+ modulo(y+1_pInt,grid(2)) * grid(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,4,1,e) = z * gridLocal(1) * gridLocal(2) &
+ modulo(y-1_pInt,gridLocal(2)) * gridLocal(1) &
mesh_ipNeighborhood(1,4,1,e) = z * grid(1) * grid(2) &
+ modulo(y-1_pInt,grid(2)) * grid(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,gridLocal(3)) * gridLocal(1) * gridLocal(2) &
+ y * gridLocal(1) &
mesh_ipNeighborhood(1,5,1,e) = modulo(z+1_pInt,grid3) * grid(1) * grid(2) &
+ y * grid(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,gridLocal(3)) * gridLocal(1) * gridLocal(2) &
+ y * gridLocal(1) &
mesh_ipNeighborhood(1,6,1,e) = modulo(z-1_pInt,grid3) * grid(1) * grid(2) &
+ y * grid(1) &
+ x &
+ 1_pInt
mesh_ipNeighborhood(2,1:6,1,e) = 1_pInt

58
code/spectral_damage.f90
View File

@ -85,8 +85,8 @@ subroutine spectral_damage_init()
use DAMASK_spectral_Utilities, only: &
wgt
use mesh, only: &
gridLocal, &
gridGlobal
grid, &
grid3
use damage_nonlocal, only: &
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility
@ -112,17 +112,17 @@ subroutine spectral_damage_init()
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
do proc = 1, worldsize
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & !< cut off stencil at boundary
DMDA_STENCIL_BOX, & !< Moore (26) neighborhood around central point
gridGlobal(1),gridGlobal(2),gridGlobal(3), & !< global grid
grid(1),grid(2),grid(3), & !< global grid
1, 1, worldsize, &
1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap)
gridLocal (1),gridLocal (2),localK, & !< local grid
grid(1),grid(2),localK, & !< local grid
damage_grid,ierr) !< handle, error
CHKERRQ(ierr)
call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da
@ -150,16 +150,16 @@ subroutine spectral_damage_init()
yend = ystart + yend - 1
zend = zstart + zend - 1
call VecSet(solution,1.0,ierr); CHKERRQ(ierr)
allocate(damage_current(gridLocal(1),gridLocal(2),gridLocal(3)), source=1.0_pReal)
allocate(damage_lastInc(gridLocal(1),gridLocal(2),gridLocal(3)), source=1.0_pReal)
allocate(damage_stagInc(gridLocal(1),gridLocal(2),gridLocal(3)), source=1.0_pReal)
allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)
!--------------------------------------------------------------------------------------------------
! damage reference diffusion update
cell = 0_pInt
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell)
mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell)
@ -184,7 +184,8 @@ type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old
Utilities_maskedCompliance, &
Utilities_updateGamma
use mesh, only: &
gridLocal
grid, &
grid3
use damage_nonlocal, only: &
damage_nonlocal_putNonLocalDamage
@ -234,7 +235,7 @@ type(tSolutionState) function spectral_damage_solution(guess,timeinc,timeinc_old
!--------------------------------------------------------------------------------------------------
! updating damage state
cell = 0_pInt !< material point = 0
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt !< material point increase
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
enddo; enddo; enddo
@ -260,12 +261,13 @@ subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr)
use numerics, only: &
residualStiffness
use mesh, only: &
gridLocal
grid, &
grid3
use math, only: &
math_mul33x3
use DAMASK_spectral_Utilities, only: &
scalarField_realMPI, &
vectorField_realMPI, &
scalarField_real, &
vectorField_real, &
utilities_FFTvectorForward, &
utilities_FFTvectorBackward, &
utilities_FFTscalarForward, &
@ -295,26 +297,26 @@ subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr)
damage_current = x_scal
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_realMPI = 0.0_pReal
scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = damage_current
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = damage_current
call utilities_FFTscalarForward()
call utilities_fourierScalarGradient() !< calculate gradient of damage field
call utilities_FFTvectorBackward()
cell = 0_pInt
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
vectorField_realMPI(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
vectorField_realMPI(1:3,i,j,k))
vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward()
call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward()
cell = 0_pInt
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
call damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, damage_current(i,j,k), 1, cell)
mobility = damage_nonlocal_getMobility(1,cell)
scalarField_realMPI(i,j,k) = params%timeinc*scalarField_realMPI(i,j,k) + &
scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + &
params%timeinc*phiDot + &
mobility*damage_lastInc(i,j,k) - &
mobility*damage_current(i,j,k) + &
@ -326,13 +328,12 @@ subroutine spectral_damage_formResidual(in,x_scal,f_scal,dummy,ierr)
call utilities_FFTscalarForward()
call utilities_fourierGreenConvolution(D_ref, mobility_ref, params%timeinc)
call utilities_FFTscalarBackward()
where(scalarField_realMPI > damage_lastInc) scalarField_realMPI = damage_lastInc
where(scalarField_realMPI < residualStiffness) scalarField_realMPI = residualStiffness
where(scalarField_real > damage_lastInc) scalarField_real = damage_lastInc
where(scalarField_real < residualStiffness) scalarField_real = residualStiffness
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) - &
damage_current
f_scal = scalarField_real(1:grid(1),1:grid(2),1:grid3) - damage_current
end subroutine spectral_damage_formResidual
@ -341,7 +342,8 @@ end subroutine spectral_damage_formResidual
!--------------------------------------------------------------------------------------------------
subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime)
use mesh, only: &
gridLocal
grid, &
grid3
use DAMASK_spectral_Utilities, only: &
cutBack, &
wgt
@ -371,7 +373,7 @@ subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime)
call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = damage_current
call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr)
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
enddo; enddo; enddo
@ -382,7 +384,7 @@ subroutine spectral_damage_forward(guess,timeinc,timeinc_old,loadCaseTime)
cell = 0_pInt
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell)
mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell)

55
code/spectral_thermal.f90
View File

@ -85,8 +85,8 @@ subroutine spectral_thermal_init
use DAMASK_spectral_Utilities, only: &
wgt
use mesh, only: &
gridLocal, &
gridGlobal
grid, &
grid3
use thermal_conduction, only: &
thermal_conduction_getConductivity33, &
thermal_conduction_getMassDensity, &
@ -116,17 +116,17 @@ subroutine spectral_thermal_init
! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = gridLocal(3)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3
do proc = 1, worldsize
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
gridGlobal(1),gridGlobal(2),gridGlobal(3), & ! global grid
grid(1),grid(2),grid(3), & ! global grid
1, 1, worldsize, &
1, 0, & ! #dof (temperature field), ghost boundary width (domain overlap)
gridLocal (1),gridLocal (2),localK, & ! local grid
grid (1),grid(2),localK, & ! local grid
thermal_grid,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da
@ -142,11 +142,11 @@ subroutine spectral_thermal_init
xend = xstart + xend - 1
yend = ystart + yend - 1
zend = zstart + zend - 1
allocate(temperature_current(gridLocal(1),gridLocal(2),gridLocal(3)), source=0.0_pReal)
allocate(temperature_lastInc(gridLocal(1),gridLocal(2),gridLocal(3)), source=0.0_pReal)
allocate(temperature_stagInc(gridLocal(1),gridLocal(2),gridLocal(3)), source=0.0_pReal)
allocate(temperature_current(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(temperature_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(temperature_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
cell = 0_pInt
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
temperature_current(i,j,k) = temperature(mappingHomogenization(2,1,cell))% &
p(thermalMapping(mappingHomogenization(2,1,cell))%p(1,cell))
@ -160,7 +160,7 @@ subroutine spectral_thermal_init
cell = 0_pInt
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
D_ref = D_ref + thermal_conduction_getConductivity33(1,cell)
mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* &
@ -186,7 +186,8 @@ type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_ol
Utilities_maskedCompliance, &
Utilities_updateGamma
use mesh, only: &
gridLocal
grid, &
grid3
use thermal_conduction, only: &
thermal_conduction_putTemperatureAndItsRate
@ -236,7 +237,7 @@ type(tSolutionState) function spectral_thermal_solution(guess,timeinc,timeinc_ol
!--------------------------------------------------------------------------------------------------
! updating thermal state
cell = 0_pInt !< material point = 0
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt !< material point increase
call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), &
(temperature_current(i,j,k)-temperature_lastInc(i,j,k))/params%timeinc, &
@ -262,12 +263,13 @@ end function spectral_thermal_solution
!--------------------------------------------------------------------------------------------------
subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr)
use mesh, only: &
gridLocal
grid, &
grid3
use math, only: &
math_mul33x3
use DAMASK_spectral_Utilities, only: &
scalarField_realMPI, &
vectorField_realMPI, &
scalarField_real, &
vectorField_real, &
utilities_FFTvectorForward, &
utilities_FFTvectorBackward, &
utilities_FFTscalarForward, &
@ -298,25 +300,25 @@ subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr)
temperature_current = x_scal
!--------------------------------------------------------------------------------------------------
! evaluate polarization field
scalarField_realMPI = 0.0_pReal
scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3)) = temperature_current
scalarField_real = 0.0_pReal
scalarField_real(1:grid(1),1:grid(2),1:grid3) = temperature_current
call utilities_FFTscalarForward()
call utilities_fourierScalarGradient() !< calculate gradient of damage field
call utilities_FFTvectorBackward()
cell = 0_pInt
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
vectorField_realMPI(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, &
vectorField_realMPI(1:3,i,j,k))
vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k))
enddo; enddo; enddo
call utilities_FFTvectorForward()
call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward()
cell = 0_pInt
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell)
scalarField_realMPI(i,j,k) = params%timeinc*scalarField_realMPI(i,j,k) + &
scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + &
params%timeinc*Tdot + &
thermal_conduction_getMassDensity (1,cell)* &
thermal_conduction_getSpecificHeat(1,cell)*(temperature_lastInc(i,j,k) - &
@ -332,7 +334,7 @@ subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr)
!--------------------------------------------------------------------------------------------------
! constructing residual
f_scal = temperature_current - scalarField_realMPI(1:gridLocal(1),1:gridLocal(2),1:gridLocal(3))
f_scal = temperature_current - scalarField_real(1:grid(1),1:grid(2),1:grid3)
end subroutine spectral_thermal_formResidual
@ -341,7 +343,8 @@ end subroutine spectral_thermal_formResidual
!--------------------------------------------------------------------------------------------------
subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime)
use mesh, only: &
gridLocal
grid, &
grid3
use DAMASK_spectral_Utilities, only: &
cutBack, &
wgt
@ -373,7 +376,7 @@ subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime)
call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with
x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current
call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr)
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt !< material point increase
call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), &
(temperature_current(i,j,k) - &
@ -387,7 +390,7 @@ subroutine spectral_thermal_forward(guess,timeinc,timeinc_old,loadCaseTime)
cell = 0_pInt
D_ref = 0.0_pReal
mobility_ref = 0.0_pReal
do k = 1_pInt, gridLocal(3); do j = 1_pInt, gridLocal(2); do i = 1_pInt,gridLocal(1)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1)
cell = cell + 1_pInt
D_ref = D_ref + thermal_conduction_getConductivity33(1,cell)
mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* &