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

Merge branch 'SolverRestart-in-HDF5' into 'development' 

Solver restart in hdf5

See merge request damask/DAMASK!74
This commit is contained in:
Karo 2019-04-30 19:18:34 +02:00
parent c3a7b7b068 5a0e408d59
commit 1c762860f8
6 changed files with 966 additions and 1031 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/CPFEM2.f90
View File

@ -257,7 +257,7 @@ subroutine CPFEM_age()
write(6,'(a)') '<< CPFEM >> writing restart variables of last converged step to hdf5 file'
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','a')
call HDF5_write(fileHandle,material_phase, 'recordedPhase')
call HDF5_write(fileHandle,crystallite_F0, 'convergedF')

4
src/grid/grid_damage_spectral.f90
View File

@ -64,7 +64,6 @@ subroutine grid_damage_spectral_init
worldsize, &
petsc_options
implicit none
PetscInt, dimension(worldsize) :: localK
integer :: i, j, k, cell
DM :: damage_grid
@ -164,7 +163,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old,loadCaseTime) result(
use damage_nonlocal, only: &
damage_nonlocal_putNonLocalDamage
implicit none
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
@ -236,7 +234,6 @@ subroutine grid_damage_spectral_forward
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility
implicit none
integer :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
@ -301,7 +298,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension( &

123
src/grid/grid_mech_FEM.f90
View File

@ -7,10 +7,11 @@
module grid_mech_FEM
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use DAMASK_interface
use HDF5_utilities
use PETScdmda
use PETScsnes
use prec, only: &
pInt, &
pReal
use math, only: &
math_I3
@ -59,8 +60,8 @@ module grid_mech_FEM
real(pReal), private :: &
err_BC !< deviation from stress BC
integer(pInt), private :: &
totalIter = 0_pInt !< total iteration in current increment
integer, private :: &
totalIter = 0 !< total iteration in current increment
public :: &
grid_mech_FEM_init, &
@ -98,7 +99,6 @@ subroutine grid_mech_FEM_init
use math, only: &
math_invSym3333
implicit none
real(pReal) :: HGCoeff = 0e-2_pReal
PetscInt, dimension(:), allocatable :: localK
real(pReal), dimension(3,3) :: &
@ -113,8 +113,8 @@ subroutine grid_mech_FEM_init
1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, &
1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8])
PetscErrorCode :: ierr
integer(pInt) :: rank
integer :: fileUnit
integer :: rank
integer(HID_T) :: fileHandle
character(len=1024) :: rankStr
real(pReal), dimension(3,3,3,3) :: devNull
PetscScalar, pointer, dimension(:,:,:,:) :: &
@ -165,8 +165,8 @@ subroutine grid_mech_FEM_init
CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr)
CHKERRQ(ierr)
call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr) ! specify custom convergence check function "_converged"
call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
@ -203,31 +203,24 @@ subroutine grid_mech_FEM_init
restart: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('F_aim')
read(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc')
read(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr))
read(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('u'//trim(rankStr))
read(fileUnit) u_current; close (fileUnit)
fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr))
read(fileUnit) u_lastInc; close (fileUnit)
call HDF5_read(fileHandle,F_aim, 'F_aim')
call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_read(fileHandle,F, 'F')
call HDF5_read(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_read(fileHandle,u_current, 'u')
call HDF5_read(fileHandle,u_lastInc, 'u_lastInc')
elseif (restartInc == 0) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restart
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateIPcoords(F)
call Utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
call utilities_updateIPcoords(F)
call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
@ -236,12 +229,11 @@ subroutine grid_mech_FEM_init
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr)
CHKERRQ(ierr)
restartRead: if (restartInc > 0_pInt) then
restartRead: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
read(fileUnit) C_volAvgLastInc; close(fileUnit)
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
endif restartRead
end subroutine grid_mech_FEM_init
@ -260,8 +252,6 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
character(len=*), intent(in) :: &
@ -274,7 +264,6 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
real(pReal), dimension(3,3), intent(in) :: rotation_BC
type(tSolutionState) :: &
solution
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
@ -284,7 +273,7 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
!--------------------------------------------------------------------------------------------------
! set module wide available data
params%stress_mask = stress_BC%maskFloat
@ -335,7 +324,6 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
use FEsolving, only: &
restartWrite
implicit none
logical, intent(in) :: &
guess
real(pReal), intent(in) :: &
@ -348,7 +336,7 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
real(pReal), dimension(3,3), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
integer :: fileUnit
integer(HID_T) :: fileHandle
character(len=32) :: rankStr
PetscScalar, pointer, dimension(:,:,:,:) :: &
u_current,u_lastInc
@ -361,35 +349,27 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim','w')
write(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w')
write(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
if (restartWrite) then
write(6,'(/,a)') ' writing converged results for restart';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('u'//trim(rankStr),'w')
write(fileUnit) u_current; close (fileUnit)
fileUnit = IO_open_jobFile_binary('u_lastInc'//trim(rankStr),'w')
write(fileUnit) u_lastInc; close (fileUnit)
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
call HDF5_write(fileHandle,F_aim, 'F_aim')
call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc')
call HDF5_write(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_write(fileHandle,F, 'F')
call HDF5_write(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_write(fileHandle,u_current, 'u')
call HDF5_write(fileHandle,u_lastInc, 'u_lastInc')
call HDF5_write(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
endif
call CPFEM_age() ! age state and kinematics
call CPFEM_age ! age state and kinematics
call utilities_updateIPcoords(F)
C_volAvgLastInc = C_volAvg
@ -441,7 +421,7 @@ end subroutine grid_mech_FEM_forward
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,ierr)
use mesh
use spectral_utilities
use numerics, only: &
@ -454,12 +434,11 @@ use spectral_utilities
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, & ! not used
snorm, & ! not used
PetscInt, intent(in) :: PETScIter
PetscReal, intent(in) :: &
devNull1, &
devNull2, &
fnorm
SNESConvergedReason :: reason
PetscObject :: dummy
@ -473,7 +452,6 @@ use spectral_utilities
divTol = max(maxval(abs(P_av))*err_div_tolRel ,err_div_tolAbs)
BCTol = max(maxval(abs(P_av))*err_stress_tolRel,err_stress_tolAbs)
if ((totalIter >= itmin .and. &
all([ err_div/divTol, &
err_BC /BCTol ] < 1.0_pReal)) &
@ -501,7 +479,8 @@ end subroutine converged
!--------------------------------------------------------------------------------------------------
!> @brief forms the residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(da_local,x_local,f_local,dummy,ierr)
subroutine formResidual(da_local,x_local, &
f_local,dummy,ierr)
use numerics, only: &
itmax, &
itmin
@ -525,7 +504,6 @@ subroutine formResidual(da_local,x_local,f_local,dummy,ierr)
use homogenization, only: &
materialpoint_dPdF
implicit none
DM :: da_local
Vec :: x_local, f_local
PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal
@ -544,12 +522,12 @@ subroutine formResidual(da_local,x_local,f_local,dummy,ierr)
call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1_pInt ! new increment
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1_pInt
totalIter = totalIter + 1
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter+1, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
@ -644,7 +622,6 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
use homogenization, only: &
materialpoint_dPdF
implicit none
DM :: da_local
Vec :: x_local, coordinates

103
src/grid/grid_mech_spectral_basic.f90
View File

@ -7,6 +7,8 @@
module grid_mech_spectral_basic
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use DAMASK_interface
use HDF5_utilities
use PETScdmda
use PETScsnes
use prec, only: &
@ -25,8 +27,7 @@ module grid_mech_spectral_basic
type(tSolutionParams), private :: params
type, private :: tNumerics
logical :: &
update_gamma !< update gamma operator with current stiffness
logical :: update_gamma !< update gamma operator with current stiffness
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
@ -40,8 +41,8 @@ module grid_mech_spectral_basic
!--------------------------------------------------------------------------------------------------
! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, &
Fdot
F_lastInc, & !< field of previous compatible deformation gradients
Fdot !< field of assumed rate of compatible deformation gradient
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
@ -99,15 +100,13 @@ subroutine grid_mech_spectral_basic_init
use spectral_utilities, only: &
utilities_constitutiveResponse, &
utilities_updateGamma, &
utilities_updateIPcoords, &
wgt
utilities_updateIPcoords
use mesh, only: &
grid, &
grid3
use math, only: &
math_invSym3333
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
@ -116,6 +115,7 @@ subroutine grid_mech_spectral_basic_init
PetscScalar, pointer, dimension(:,:,:,:) :: &
F ! pointer to solution data
PetscInt, dimension(worldsize) :: localK
integer(HID_T) :: fileHandle
integer :: fileUnit
character(len=1024) :: rankStr
@ -174,19 +174,14 @@ subroutine grid_mech_spectral_basic_init
restart: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('F_aim')
read(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc')
read(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr))
read(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
call HDF5_read(fileHandle,F_aim, 'F_aim')
call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_read(fileHandle,F, 'F')
call HDF5_read(fileHandle,F_lastInc, 'F_lastInc')
elseif (restartInc == 0) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
@ -203,15 +198,15 @@ subroutine grid_mech_spectral_basic_init
restartRead: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') 'reading more values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
read(fileUnit) C_volAvgLastInc; close(fileUnit)
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
fileUnit = IO_open_jobFile_binary('C_ref')
read(fileUnit) C_minMaxAvg; close(fileUnit)
endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.true.)
call utilities_updateGamma(C_minMaxAvg,.true.)
end subroutine grid_mech_spectral_basic_init
@ -228,8 +223,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
character(len=*), intent(in) :: &
@ -251,8 +244,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (num%update_gamma) call Utilities_updateGamma(C_minMaxAvg,restartWrite)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg,restartWrite)
!--------------------------------------------------------------------------------------------------
! set module wide available data
@ -306,7 +299,6 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
use FEsolving, only: &
restartWrite
implicit none
logical, intent(in) :: &
guess
real(pReal), intent(in) :: &
@ -321,7 +313,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: F
integer :: fileUnit
integer(HID_T) :: fileHandle
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr)
@ -332,28 +324,23 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim','w')
write(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w')
write(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
if (restartWrite) then
write(6,'(/,a)') ' writing converged results for restart';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
call HDF5_write(fileHandle,F_aim, 'F_aim')
call HDF5_write(fileHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_write(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_write(fileHandle,F, 'F')
call HDF5_write(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_write(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_write(fileHandle,C_minMaxAvg, 'C_minMaxAvg')
call HDF5_closeFile(fileHandle)
endif
call CPFEM_age ! age state and kinematics
@ -399,7 +386,7 @@ end subroutine grid_mech_spectral_basic_forward
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
@ -410,13 +397,12 @@ subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, & ! not used
snorm, & ! not used
fnorm ! not used
PetscInt, intent(in) :: PETScIter
PetscReal, intent(in) :: &
devNull1, &
devNull2, &
devNull3
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
@ -452,7 +438,7 @@ end subroutine converged
!--------------------------------------------------------------------------------------------------
!> @brief forms the basic residual vector
!> @brief forms the residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in, F, &
residuum, dummy, ierr)
@ -481,7 +467,6 @@ subroutine formResidual(in, F, &
use FEsolving, only: &
terminallyIll
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work)
PetscScalar, dimension(3,3,XG_RANGE,YG_RANGE,ZG_RANGE), &
intent(in) :: F !< deformation gradient field
@ -515,7 +500,7 @@ subroutine formResidual(in, F, &
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
call Utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory)
call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, &
F,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)

107
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -7,6 +7,8 @@
module grid_mech_spectral_polarisation
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use DAMASK_interface
use HDF5_utilities
use PETScdmda
use PETScsnes
use prec, only: &
@ -25,8 +27,7 @@ module grid_mech_spectral_polarisation
type(tSolutionParams), private :: params
type, private :: tNumerics
logical :: &
update_gamma !< update gamma operator with current stiffness
logical :: update_gamma !< update gamma operator with current stiffness
end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name?
@ -105,15 +106,13 @@ subroutine grid_mech_spectral_polarisation_init
use spectral_utilities, only: &
utilities_constitutiveResponse, &
utilities_updateGamma, &
utilities_updateIPcoords, &
wgt
utilities_updateIPcoords
use mesh, only: &
grid, &
grid3
use math, only: &
math_invSym3333
implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
@ -124,6 +123,7 @@ subroutine grid_mech_spectral_polarisation_init
F, & ! specific (sub)pointer
F_tau ! specific (sub)pointer
PetscInt, dimension(worldsize) :: localK
integer(HID_T) :: fileHandle
integer :: fileUnit
character(len=1024) :: rankStr
@ -183,23 +183,16 @@ subroutine grid_mech_spectral_polarisation_init
restart: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('F_aim')
read(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc')
read(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot')
read(fileUnit) F_aimDot; close(fileUnit)
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr))
read(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr))
read(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr))
read(fileUnit) F_tau; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr))
read(fileUnit) F_tau_lastInc; close (fileUnit)
call HDF5_read(fileHandle,F_aim, 'F_aim')
call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_read(fileHandle,F, 'F')
call HDF5_read(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_read(fileHandle,F_tau, 'F_tau')
call HDF5_read(fileHandle,F_tau_lastInc,'F_tau_lastInc')
elseif (restartInc == 0) then restart
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
@ -218,15 +211,15 @@ subroutine grid_mech_spectral_polarisation_init
restartRead: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file'
fileUnit = IO_open_jobFile_binary('C_volAvg')
read(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv')
read(fileUnit) C_volAvgLastInc; close(fileUnit)
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
fileUnit = IO_open_jobFile_binary('C_ref')
read(fileUnit) C_minMaxAvg; close(fileUnit)
endif restartRead
call Utilities_updateGamma(C_minMaxAvg,.true.)
call utilities_updateGamma(C_minMaxAvg,.true.)
C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg)
@ -247,8 +240,6 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
restartWrite, &
terminallyIll
implicit none
!--------------------------------------------------------------------------------------------------
! input data for solution
character(len=*), intent(in) :: &
@ -270,7 +261,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
S = Utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (num%update_gamma) then
call utilities_updateGamma(C_minMaxAvg,restartWrite)
C_scale = C_minMaxAvg
@ -330,7 +321,6 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
use FEsolving, only: &
restartWrite
implicit none
logical, intent(in) :: &
guess
real(pReal), intent(in) :: &
@ -347,7 +337,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
integer :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33
integer :: fileUnit
integer(HID_T) :: fileHandle
character(len=32) :: rankStr
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
@ -360,32 +350,24 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
else
!--------------------------------------------------------------------------------------------------
! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set?
write(6,'(/,a)') ' writing converged results for restart'
flush(6)
if (worldrank == 0) then
fileUnit = IO_open_jobFile_binary('C_volAvg','w')
write(fileUnit) C_volAvg; close(fileUnit)
fileUnit = IO_open_jobFile_binary('C_volAvgLastInv','w')
write(fileUnit) C_volAvgLastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim','w')
write(fileUnit) F_aim; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aim_lastInc','w')
write(fileUnit) F_aim_lastInc; close(fileUnit)
fileUnit = IO_open_jobFile_binary('F_aimDot','w')
write(fileUnit) F_aimDot; close(fileUnit)
endif
if (restartWrite) then
write(6,'(/,a)') ' writing converged results for restart';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank
fileUnit = IO_open_jobFile_binary('F'//trim(rankStr),'w')
write(fileUnit) F; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_lastInc'//trim(rankStr),'w')
write(fileUnit) F_lastInc; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau'//trim(rankStr),'w')
write(fileUnit) F_tau; close (fileUnit)
fileUnit = IO_open_jobFile_binary('F_tau_lastInc'//trim(rankStr),'w')
write(fileUnit) F_tau_lastInc; close (fileUnit)
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w')
call HDF5_write(fileHandle,F_aim, 'F_aim')
call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc')
call HDF5_write(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_write(fileHandle,F, 'F')
call HDF5_write(fileHandle,F_lastInc, 'F_lastInc')
call HDF5_write(fileHandle,F_tau, 'F_tau')
call HDF5_write(fileHandle,F_tau_lastInc, 'F_tau_lastInc')
call HDF5_write(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
endif
call CPFEM_age ! age state and kinematics
@ -451,7 +433,7 @@ end subroutine grid_mech_spectral_polarisation_forward
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr)
use numerics, only: &
itmax, &
itmin, &
@ -464,13 +446,12 @@ subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
use FEsolving, only: &
terminallyIll
implicit none
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: &
xnorm, & ! not used
snorm, & ! not used
fnorm ! not used
PetscInt, intent(in) :: PETScIter
PetscReal, intent(in) :: &
devNull1, &
devNull2, &
devNull3
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
@ -511,7 +492,7 @@ end subroutine converged
!--------------------------------------------------------------------------------------------------
!> @brief forms the polarisation residual vector
!> @brief forms the residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in, FandF_tau, &
residuum, dummy,ierr)
@ -548,7 +529,6 @@ subroutine formResidual(in, FandF_tau, &
use FEsolving, only: &
terminallyIll
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work)
PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), &
target, intent(in) :: FandF_tau
@ -635,6 +615,7 @@ subroutine formResidual(in, FandF_tau, &
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise
call utilities_FFTtensorForward
err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress
!--------------------------------------------------------------------------------------------------
! constructing residual
e = 0

4
src/grid/grid_thermal_spectral.f90
View File

@ -69,7 +69,6 @@ subroutine grid_thermal_spectral_init
worldsize, &
petsc_options
implicit none
PetscInt, dimension(worldsize) :: localK
integer :: i, j, k, cell
DM :: thermal_grid
@ -167,7 +166,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old,loadCaseTime) result
use thermal_conduction, only: &
thermal_conduction_putTemperatureAndItsRate
implicit none
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
@ -242,7 +240,6 @@ subroutine grid_thermal_spectral_forward
thermal_conduction_getMassDensity, &
thermal_conduction_getSpecificHeat
implicit none
integer :: i, j, k, cell
DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal
@ -311,7 +308,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
thermal_conduction_getMassDensity, &
thermal_conduction_getSpecificHeat
implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: &
in
PetscScalar, dimension( &