Allreduce better suited than Bcast

This commit is contained in:
Martin Diehl 2019-03-09 09:54:33 +01:00
parent a5dad0ca5c
commit 6f3771f6c4
4 changed files with 296 additions and 325 deletions

View File

@ -11,14 +11,9 @@ module spectral_damage
use prec, only: & use prec, only: &
pInt, & pInt, &
pReal pReal
use math, only: &
math_I3
use spectral_utilities, only: & use spectral_utilities, only: &
tSolutionState, & tSolutionState, &
tSolutionParams tSolutionParams
use numerics, only: &
worldrank, &
worldsize
implicit none implicit none
private private
@ -42,7 +37,7 @@ module spectral_damage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc. ! reference diffusion tensor, mobility etc.
integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment integer(pInt), private :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3), private :: D_ref real(pReal), dimension(3,3), private :: D_ref
real(pReal), private :: mobility_ref real(pReal), private :: mobility_ref
@ -57,93 +52,94 @@ contains
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info !> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine spectral_damage_init() subroutine spectral_damage_init()
use IO, only: & use IO, only: &
IO_intOut IO_intOut
use spectral_utilities, only: & use spectral_utilities, only: &
wgt wgt
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
use damage_nonlocal, only: & use damage_nonlocal, only: &
damage_nonlocal_getDiffusion33, & damage_nonlocal_getDiffusion33, &
damage_nonlocal_getMobility damage_nonlocal_getMobility
use numerics, only: &
implicit none worldrank, &
PetscInt, dimension(:), allocatable :: localK worldsize
integer(pInt) :: proc
integer(pInt) :: i, j, k, cell implicit none
DM :: damage_grid PetscInt, dimension(worldsize) :: localK
Vec :: uBound, lBound integer :: i, j, k, cell
PetscErrorCode :: ierr DM :: damage_grid
character(len=100) :: snes_type Vec :: uBound, lBound
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>' character(len=100) :: snes_type
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press, '
write(6,'(a,/)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018 ' write(6,'(/,a)') ' <<<+- spectral_damage init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press, '
write(6,'(a,/)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018 '
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,damage_snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(damage_snes,'damage_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 localK = 0
do proc = 1, worldsize localK(worldrank+1) = grid3
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
enddo call DMDACreate3D(PETSC_COMM_WORLD, &
call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & !< cut off stencil at boundary
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & !< cut off stencil at boundary DMDA_STENCIL_BOX, & !< Moore (26) neighborhood around central point
DMDA_STENCIL_BOX, & !< Moore (26) neighborhood around central point grid(1),grid(2),grid(3), & !< global grid
grid(1),grid(2),grid(3), & !< global grid 1, 1, worldsize, &
1, 1, worldsize, & 1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap)
1, 0, & !< #dof (damage phase field), ghost boundary width (domain overlap) [grid(1)],[grid(2)],localK, & !< local grid
[grid(1)],[grid(2)],localK, & !< local grid damage_grid,ierr) !< handle, error
damage_grid,ierr) !< handle, error CHKERRQ(ierr)
CHKERRQ(ierr) call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da
call SNESSetDM(damage_snes,damage_grid,ierr); CHKERRQ(ierr) !< connect snes to da call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr)
call DMsetFromOptions(damage_grid,ierr); CHKERRQ(ierr) call DMsetUp(damage_grid,ierr); CHKERRQ(ierr)
call DMsetUp(damage_grid,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(damage_grid,solution,ierr); CHKERRQ(ierr) !< global solution vector (grid x 1, i.e. every def grad tensor)
call DMCreateGlobalVector(damage_grid,solution,ierr); CHKERRQ(ierr) !< global solution vector (grid x 1, i.e. every def grad tensor) call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,spectral_damage_formResidual,&
call DMDASNESSetFunctionLocal(damage_grid,INSERT_VALUES,spectral_damage_formResidual,& PETSC_NULL_SNES,ierr) !< residual vector of same shape as solution vector
PETSC_NULL_SNES,ierr) !< residual vector of same shape as solution vector CHKERRQ(ierr)
CHKERRQ(ierr) call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional CLI arguments
call SNESSetFromOptions(damage_snes,ierr); CHKERRQ(ierr) !< pull it all together with additional CLI arguments call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr)
call SNESGetType(damage_snes,snes_type,ierr); CHKERRQ(ierr) if (trim(snes_type) == 'vinewtonrsls' .or. &
if (trim(snes_type) == 'vinewtonrsls' .or. & trim(snes_type) == 'vinewtonssls') then
trim(snes_type) == 'vinewtonssls') then call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr)
call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr) call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr)
call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr) call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc.
call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc. call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) endif
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr) call DMDAGetCorners(damage_grid,xstart,ystart,zstart,xend,yend,zend,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
xend = xstart + xend - 1 xend = xstart + xend - 1
yend = ystart + yend - 1 yend = ystart + yend - 1
zend = zstart + zend - 1 zend = zstart + zend - 1
call VecSet(solution,1.0_pReal,ierr); CHKERRQ(ierr) call VecSet(solution,1.0_pReal,ierr); CHKERRQ(ierr)
allocate(damage_current(grid(1),grid(2),grid3), 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_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_stagInc(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 ! damage reference diffusion update
cell = 0_pInt cell = 0_pInt
D_ref = 0.0_pReal D_ref = 0.0_pReal
mobility_ref = 0.0_pReal mobility_ref = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1_pInt cell = cell + 1
D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell) D_ref = D_ref + damage_nonlocal_getDiffusion33(1,cell)
mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell) mobility_ref = mobility_ref + damage_nonlocal_getMobility(1,cell)
enddo; enddo; enddo enddo; enddo; enddo
D_ref = D_ref*wgt D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine spectral_damage_init end subroutine spectral_damage_init
@ -151,74 +147,69 @@ end subroutine spectral_damage_init
!> @brief solution for the spectral damage scheme with internal iterations !> @brief solution for the spectral damage scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadCaseTime) type(tSolutionState) function spectral_damage_solution(timeinc,timeinc_old,loadCaseTime)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
err_damage_tolAbs, & err_damage_tolAbs, &
err_damage_tolRel err_damage_tolRel
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
use damage_nonlocal, only: & use damage_nonlocal, only: &
damage_nonlocal_putNonLocalDamage damage_nonlocal_putNonLocalDamage
implicit none
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
loadCaseTime !< remaining time of current load case
implicit none integer :: i, j, k, cell
PetscInt ::position
!-------------------------------------------------------------------------------------------------- PetscReal :: minDamage, maxDamage, stagNorm, solnNorm
! input data for solution PetscErrorCode :: ierr
real(pReal), intent(in) :: & SNESConvergedReason :: reason
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
loadCaseTime !< remaining time of current load case
integer(pInt) :: i, j, k, cell
PetscInt ::position
PetscReal :: minDamage, maxDamage, stagNorm, solnNorm
!--------------------------------------------------------------------------------------------------
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
spectral_damage_solution%converged =.false. spectral_damage_solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old params%timeincOld = timeinc_old
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) call SNESSolve(damage_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)
if (reason < 1) then if (reason < 1) then
spectral_damage_solution%converged = .false. spectral_damage_solution%converged = .false.
spectral_damage_solution%iterationsNeeded = itmax spectral_damage_solution%iterationsNeeded = itmax
else else
spectral_damage_solution%converged = .true. spectral_damage_solution%converged = .true.
spectral_damage_solution%iterationsNeeded = totalIter spectral_damage_solution%iterationsNeeded = totalIter
endif endif
stagNorm = maxval(abs(damage_current - damage_stagInc)) stagNorm = maxval(abs(damage_current - damage_stagInc))
solnNorm = maxval(abs(damage_current)) solnNorm = maxval(abs(damage_current))
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
damage_stagInc = damage_current damage_stagInc = damage_current
spectral_damage_solution%stagConverged = stagNorm < err_damage_tolAbs & spectral_damage_solution%stagConverged = stagNorm < err_damage_tolAbs &
.or. stagNorm < err_damage_tolRel*solnNorm .or. stagNorm < err_damage_tolRel*solnNorm
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updating damage state ! updating damage state
cell = 0_pInt !< material point = 0 cell = 0
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1_pInt !< material point increase cell = cell + 1
call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell) call damage_nonlocal_putNonLocalDamage(damage_current(i,j,k),1,cell)
enddo; enddo; enddo enddo; enddo; enddo
call VecMin(solution,position,minDamage,ierr); CHKERRQ(ierr) call VecMin(solution,position,minDamage,ierr); CHKERRQ(ierr)
call VecMax(solution,position,maxDamage,ierr); CHKERRQ(ierr) call VecMax(solution,position,maxDamage,ierr); CHKERRQ(ierr)
if (spectral_damage_solution%converged) & if (spectral_damage_solution%converged) &
write(6,'(/,a)') ' ... nonlocal damage converged .....................................' write(6,'(/,a)') ' ... nonlocal damage converged .....................................'
write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',& write(6,'(/,a,f8.6,2x,f8.6,2x,f8.6,/)',advance='no') ' Minimum|Maximum|Delta Damage = ',&
minDamage, maxDamage, stagNorm minDamage, maxDamage, stagNorm
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)
end function spectral_damage_solution end function spectral_damage_solution

View File

@ -108,8 +108,8 @@ subroutine basic_init
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: F PetscScalar, pointer, dimension(:,:,:,:) :: F
PetscInt, dimension(:), allocatable :: localK PetscInt, dimension(worldsize) :: localK
integer :: proc, fileUnit integer :: fileUnit
character(len=1024) :: rankStr character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>'
@ -125,10 +125,9 @@ subroutine basic_init
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 localK = 0
do proc = 1, worldsize !ToDo: there are smarter options in MPI localK(worldrank+1) = grid3
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, & call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point

View File

@ -118,8 +118,8 @@ subroutine Polarisation_init
FandF_tau, & ! overall pointer to solution data FandF_tau, & ! overall pointer to solution data
F, & ! specific (sub)pointer F, & ! specific (sub)pointer
F_tau ! specific (sub)pointer F_tau ! specific (sub)pointer
PetscInt, dimension(:), allocatable :: localK PetscInt, dimension(worldsize) :: localK
integer :: proc, fileUnit integer :: fileUnit
character(len=1024) :: rankStr character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>' write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverPolarisation init -+>>>'
@ -137,10 +137,9 @@ subroutine Polarisation_init
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 localK = 0
do proc = 1, worldsize !ToDo: there are smarter options in MPI localK(worldrank+1) = grid3
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
enddo
call DMDACreate3d(PETSC_COMM_WORLD, & call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point

View File

@ -9,16 +9,10 @@ module spectral_thermal
use PETScdmda use PETScdmda
use PETScsnes use PETScsnes
use prec, only: & use prec, only: &
pInt, &
pReal pReal
use math, only: &
math_I3
use spectral_utilities, only: & use spectral_utilities, only: &
tSolutionState, & tSolutionState, &
tSolutionParams tSolutionParams
use numerics, only: &
worldrank, &
worldsize
implicit none implicit none
private private
@ -42,7 +36,7 @@ module spectral_thermal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reference diffusion tensor, mobility etc. ! reference diffusion tensor, mobility etc.
integer(pInt), private :: totalIter = 0_pInt !< total iteration in current increment integer, private :: totalIter = 0 !< total iteration in current increment
real(pReal), dimension(3,3), private :: D_ref real(pReal), dimension(3,3), private :: D_ref
real(pReal), private :: mobility_ref real(pReal), private :: mobility_ref
@ -57,104 +51,96 @@ contains
!> @brief allocates all neccessary fields and fills them with data, potentially from restart info !> @brief allocates all neccessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine spectral_thermal_init subroutine spectral_thermal_init
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800 use spectral_utilities, only: &
use, intrinsic :: iso_fortran_env, only: & wgt
compiler_version, & use mesh, only: &
compiler_options grid, &
#endif grid3
use IO, only: & use thermal_conduction, only: &
IO_timeStamp thermal_conduction_getConductivity33, &
use spectral_utilities, only: & thermal_conduction_getMassDensity, &
wgt thermal_conduction_getSpecificHeat
use mesh, only: & use material, only: &
grid, & mappingHomogenization, &
grid3 temperature, &
use thermal_conduction, only: & thermalMapping
thermal_conduction_getConductivity33, & use numerics, only: &
thermal_conduction_getMassDensity, & worldrank, &
thermal_conduction_getSpecificHeat worldsize
use material, only: &
mappingHomogenization, &
temperature, &
thermalMapping
implicit none
integer(pInt), dimension(:), allocatable :: localK
integer(pInt) :: proc
integer(pInt) :: i, j, k, cell
DM :: thermal_grid
PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>' implicit none
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press,' integer, dimension(worldsize) :: localK
write(6,'(/,a)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018' integer :: i, j, k, cell
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() DM :: thermal_grid
#include "compilation_info.f90" PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- spectral_thermal init -+>>>'
write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, volume in press,'
write(6,'(/,a)') ' chapter Spectral Solvers for Crystal Plasticity and Multi-Physics Simulations. Springer, 2018'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,thermal_snes,ierr); CHKERRQ(ierr)
call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(thermal_snes,'thermal_',ierr);CHKERRQ(ierr)
allocate(localK(worldsize), source = 0); localK(worldrank+1) = grid3 localK = 0
do proc = 1, worldsize localK(worldrank+1) = grid3
call MPI_Bcast(localK(proc),1,MPI_INTEGER,proc-1,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
enddo call DMDACreate3D(PETSC_COMM_WORLD, &
call DMDACreate3D(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point grid(1),grid(2),grid(3), & ! global grid
grid(1),grid(2),grid(3), & ! global grid 1, 1, worldsize, &
1, 1, worldsize, & 1, 0, & !< #dof (thermal phase field), ghost boundary width (domain overlap)
1, 0, & !< #dof (thermal phase field), ghost boundary width (domain overlap) [grid(1)],[grid(2)],localK, & !< local grid
[grid(1)],[grid(2)],localK, & !< local grid thermal_grid,ierr) !< handle, error
thermal_grid,ierr) !< handle, error CHKERRQ(ierr)
CHKERRQ(ierr) call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da
call SNESSetDM(thermal_snes,thermal_grid,ierr); CHKERRQ(ierr) ! connect snes to da call DMsetFromOptions(thermal_grid,ierr); CHKERRQ(ierr)
call DMsetFromOptions(thermal_grid,ierr); CHKERRQ(ierr) call DMsetUp(thermal_grid,ierr); CHKERRQ(ierr)
call DMsetUp(thermal_grid,ierr); CHKERRQ(ierr) call DMCreateGlobalVector(thermal_grid,solution ,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor)
call DMCreateGlobalVector(thermal_grid,solution ,ierr); CHKERRQ(ierr) ! global solution vector (grid x 1, i.e. every def grad tensor) call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,spectral_thermal_formResidual,&
call DMDASNESSetFunctionLocal(thermal_grid,INSERT_VALUES,spectral_thermal_formResidual,& PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector CHKERRQ(ierr)
CHKERRQ(ierr) call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
call SNESSetFromOptions(thermal_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,ierr) call DMDAGetCorners(thermal_grid,xstart,ystart,zstart,xend,yend,zend,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
xend = xstart + xend - 1 xend = xstart + xend - 1
yend = ystart + yend - 1 yend = ystart + yend - 1
zend = zstart + zend - 1 zend = zstart + zend - 1
allocate(temperature_current(grid(1),grid(2),grid3), 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_lastInc(grid(1),grid(2),grid3), source=0.0_pReal)
allocate(temperature_stagInc(grid(1),grid(2),grid3), source=0.0_pReal) allocate(temperature_stagInc(grid(1),grid(2),grid3), source=0.0_pReal)
cell = 0_pInt cell = 0
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1_pInt cell = cell + 1
temperature_current(i,j,k) = temperature(mappingHomogenization(2,1,cell))% & temperature_current(i,j,k) = temperature(mappingHomogenization(2,1,cell))% &
p(thermalMapping(mappingHomogenization(2,1,cell))%p(1,cell)) p(thermalMapping(mappingHomogenization(2,1,cell))%p(1,cell))
temperature_lastInc(i,j,k) = temperature_current(i,j,k) temperature_lastInc(i,j,k) = temperature_current(i,j,k)
temperature_stagInc(i,j,k) = temperature_current(i,j,k) temperature_stagInc(i,j,k) = temperature_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
call DMDAVecGetArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with call DMDAVecGetArrayF90(thermal_grid,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 x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current
call DMDAVecRestoreArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(thermal_grid,solution,x_scal,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! thermal reference diffusion update ! thermal reference diffusion update
cell = 0_pInt cell = 0
D_ref = 0.0_pReal D_ref = 0.0_pReal
mobility_ref = 0.0_pReal mobility_ref = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1_pInt cell = cell + 1
D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) D_ref = D_ref + thermal_conduction_getConductivity33(1,cell)
mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* & mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* &
thermal_conduction_getSpecificHeat(1,cell) thermal_conduction_getSpecificHeat(1,cell)
enddo; enddo; enddo enddo; enddo; enddo
D_ref = D_ref*wgt D_ref = D_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,D_ref,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
mobility_ref = mobility_ref*wgt mobility_ref = mobility_ref*wgt
call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,mobility_ref,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine spectral_thermal_init end subroutine spectral_thermal_init
@ -162,76 +148,72 @@ end subroutine spectral_thermal_init
!> @brief solution for the spectral thermal scheme with internal iterations !> @brief solution for the spectral thermal scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,loadCaseTime) type(tSolutionState) function spectral_thermal_solution(timeinc,timeinc_old,loadCaseTime)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
err_thermal_tolAbs, & err_thermal_tolAbs, &
err_thermal_tolRel err_thermal_tolRel
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
use thermal_conduction, only: & use thermal_conduction, only: &
thermal_conduction_putTemperatureAndItsRate thermal_conduction_putTemperatureAndItsRate
implicit none
implicit none real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
loadCaseTime !< remaining time of current load case
integer :: i, j, k, cell
PetscInt :: position
PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm
!-------------------------------------------------------------------------------------------------- PetscErrorCode :: ierr
! input data for solution SNESConvergedReason :: reason
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old, & !< increment in time of last increment
loadCaseTime !< remaining time of current load case
integer(pInt) :: i, j, k, cell
PetscInt :: position
PetscReal :: minTemperature, maxTemperature, stagNorm, solnNorm
!-------------------------------------------------------------------------------------------------- spectral_thermal_solution%converged =.false.
! PETSc Data
PetscErrorCode :: ierr
SNESConvergedReason :: reason
spectral_thermal_solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old params%timeincOld = timeinc_old
call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr)
if (reason < 1) then if (reason < 1) then
spectral_thermal_solution%converged = .false. spectral_thermal_solution%converged = .false.
spectral_thermal_solution%iterationsNeeded = itmax spectral_thermal_solution%iterationsNeeded = itmax
else else
spectral_thermal_solution%converged = .true. spectral_thermal_solution%converged = .true.
spectral_thermal_solution%iterationsNeeded = totalIter spectral_thermal_solution%iterationsNeeded = totalIter
endif endif
stagNorm = maxval(abs(temperature_current - temperature_stagInc)) stagNorm = maxval(abs(temperature_current - temperature_stagInc))
solnNorm = maxval(abs(temperature_current)) solnNorm = maxval(abs(temperature_current))
call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,stagNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,solnNorm,1,MPI_DOUBLE,MPI_MAX,PETSC_COMM_WORLD,ierr)
temperature_stagInc = temperature_current temperature_stagInc = temperature_current
spectral_thermal_solution%stagConverged = stagNorm < err_thermal_tolAbs & spectral_thermal_solution%stagConverged = stagNorm < err_thermal_tolAbs &
.or. stagNorm < err_thermal_tolRel*solnNorm .or. stagNorm < err_thermal_tolRel*solnNorm
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updating thermal state ! updating thermal state
cell = 0_pInt !< material point = 0 cell = 0
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1_pInt !< material point increase cell = cell + 1
call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), &
(temperature_current(i,j,k)-temperature_lastInc(i,j,k))/params%timeinc, & (temperature_current(i,j,k)-temperature_lastInc(i,j,k))/params%timeinc, &
1,cell) 1,cell)
enddo; enddo; enddo enddo; enddo; enddo
call VecMin(solution,position,minTemperature,ierr); CHKERRQ(ierr) call VecMin(solution,position,minTemperature,ierr); CHKERRQ(ierr)
call VecMax(solution,position,maxTemperature,ierr); CHKERRQ(ierr) call VecMax(solution,position,maxTemperature,ierr); CHKERRQ(ierr)
if (spectral_thermal_solution%converged) & if (spectral_thermal_solution%converged) &
write(6,'(/,a)') ' ... thermal conduction converged ..................................' write(6,'(/,a)') ' ... thermal conduction converged ..................................'
write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',& write(6,'(/,a,f8.4,2x,f8.4,2x,f8.4,/)',advance='no') ' Minimum|Maximum|Delta Temperature / K = ',&
minTemperature, maxTemperature, stagNorm minTemperature, maxTemperature, stagNorm
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)
end function spectral_thermal_solution end function spectral_thermal_solution
@ -272,7 +254,7 @@ subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr)
f_scal f_scal
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(pInt) :: i, j, k, cell integer :: i, j, k, cell
real(pReal) :: Tdot, dTdot_dT real(pReal) :: Tdot, dTdot_dT
temperature_current = x_scal temperature_current = x_scal
@ -283,18 +265,18 @@ subroutine spectral_thermal_formResidual(in,x_scal,f_scal,dummy,ierr)
call utilities_FFTscalarForward() call utilities_FFTscalarForward()
call utilities_fourierScalarGradient() !< calculate gradient of damage field call utilities_fourierScalarGradient() !< calculate gradient of damage field
call utilities_FFTvectorBackward() call utilities_FFTvectorBackward()
cell = 0_pInt cell = 0
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1_pInt cell = cell + 1
vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward() call utilities_FFTvectorForward()
call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field call utilities_fourierVectorDivergence() !< calculate damage divergence in fourier field
call utilities_FFTscalarBackward() call utilities_FFTscalarBackward()
cell = 0_pInt cell = 0
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1_pInt cell = cell + 1
call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell) call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, temperature_current(i,j,k), 1, cell)
scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + & scalarField_real(i,j,k) = params%timeinc*scalarField_real(i,j,k) + &
params%timeinc*Tdot + & params%timeinc*Tdot + &
@ -333,10 +315,10 @@ subroutine spectral_thermal_forward()
thermal_conduction_getSpecificHeat thermal_conduction_getSpecificHeat
implicit none implicit none
integer(pInt) :: i, j, k, cell integer :: i, j, k, cell
DM :: dm_local DM :: dm_local
PetscScalar, dimension(:,:,:), pointer :: x_scal PetscScalar, dimension(:,:,:), pointer :: x_scal
PetscErrorCode :: ierr PetscErrorCode :: ierr
if (cutBack) then if (cutBack) then
temperature_current = temperature_lastInc temperature_current = temperature_lastInc
@ -344,13 +326,13 @@ subroutine spectral_thermal_forward()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reverting thermal field state ! reverting thermal field state
cell = 0_pInt !< material point = 0 cell = 0
call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr) call SNESGetDM(thermal_snes,dm_local,ierr); CHKERRQ(ierr)
call DMDAVecGetArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) !< get the data out of PETSc to work with 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 x_scal(xstart:xend,ystart:yend,zstart:zend) = temperature_current
call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(dm_local,solution,x_scal,ierr); CHKERRQ(ierr)
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1_pInt !< material point increase cell = cell + 1
call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), & call thermal_conduction_putTemperatureAndItsRate(temperature_current(i,j,k), &
(temperature_current(i,j,k) - & (temperature_current(i,j,k) - &
temperature_lastInc(i,j,k))/params%timeinc, & temperature_lastInc(i,j,k))/params%timeinc, &
@ -360,11 +342,11 @@ subroutine spectral_thermal_forward()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update rate and forward last inc ! update rate and forward last inc
temperature_lastInc = temperature_current temperature_lastInc = temperature_current
cell = 0_pInt cell = 0
D_ref = 0.0_pReal D_ref = 0.0_pReal
mobility_ref = 0.0_pReal mobility_ref = 0.0_pReal
do k = 1_pInt, grid3; do j = 1_pInt, grid(2); do i = 1_pInt,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1_pInt cell = cell + 1
D_ref = D_ref + thermal_conduction_getConductivity33(1,cell) D_ref = D_ref + thermal_conduction_getConductivity33(1,cell)
mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* & mobility_ref = mobility_ref + thermal_conduction_getMassDensity(1,cell)* &
thermal_conduction_getSpecificHeat(1,cell) thermal_conduction_getSpecificHeat(1,cell)