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

unified style

This commit is contained in:
Martin Diehl 2019-04-15 15:53:46 +02:00
parent 18f9deef1a
commit 4793f964f8
2 changed files with 430 additions and 433 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
src/grid_mech_spectral_basic.f90
View File

@ -27,8 +27,7 @@ module grid_mech_spectral_basic
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
type, private :: tNumerics type, private :: tNumerics
logical :: & logical :: update_gamma !< update gamma operator with current stiffness
update_gamma !< update gamma operator with current stiffness
end type tNumerics end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
@ -42,8 +41,8 @@ module grid_mech_spectral_basic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! common pointwise data ! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: & real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, & F_lastInc, & !< field of previous compatible deformation gradients
Fdot Fdot !< field of assumed rate of compatible deformation gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
@ -101,8 +100,7 @@ subroutine grid_mech_spectral_basic_init
use spectral_utilities, only: & use spectral_utilities, only: &
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
utilities_updateGamma, & utilities_updateGamma, &
utilities_updateIPcoords, & utilities_updateIPcoords
wgt
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
@ -179,6 +177,7 @@ subroutine grid_mech_spectral_basic_init
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
call HDF5_read(fileHandle,F_aim, 'F_aim') call HDF5_read(fileHandle,F_aim, 'F_aim')
call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc') call HDF5_read(fileHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(fileHandle,F_aimDot, 'F_aimDot') call HDF5_read(fileHandle,F_aimDot, 'F_aimDot')
@ -226,7 +225,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
terminallyIll terminallyIll
implicit none implicit none
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
@ -248,8 +246,8 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! 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) call Utilities_updateGamma(C_minMaxAvg,restartWrite) if (num%update_gamma) call utilities_updateGamma(C_minMaxAvg,restartWrite)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide available data ! set module wide available data
@ -329,7 +327,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
else else
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! restart information for spectral solver ! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set? if (restartWrite) then
write(6,'(/,a)') ' writing converged results for restart';flush(6) write(6,'(/,a)') ' writing converged results for restart';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
@ -344,6 +342,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
call HDF5_write(fileHandle,C_volAvg, 'C_volAvg') call HDF5_write(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_write(fileHandle,C_minMaxAvg, 'C_minMaxAvg') call HDF5_write(fileHandle,C_minMaxAvg, 'C_minMaxAvg')
call HDF5_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
endif endif
@ -390,7 +389,7 @@ end subroutine grid_mech_spectral_basic_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convergence check !> @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: & use numerics, only: &
itmax, & itmax, &
itmin, & itmin, &
@ -403,11 +402,11 @@ subroutine converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
implicit none implicit none
SNES :: snes_local SNES :: snes_local
PetscInt :: PETScIter PetscInt, intent(in) :: PETScIter
PetscReal :: & PetscReal, intent(in) :: &
xnorm, & ! not used devNull1, &
snorm, & ! not used devNull2, &
fnorm ! not used devNull3
SNESConvergedReason :: reason SNESConvergedReason :: reason
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
@ -506,7 +505,7 @@ subroutine formResidual(in, F, &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! 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, & P_av,C_volAvg,C_minMaxAvg, &
F,params%timeinc,params%rotation_BC) F,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)

830
src/grid_mech_spectral_polarisation.f90
View File

@ -7,80 +7,79 @@
module grid_mech_spectral_polarisation module grid_mech_spectral_polarisation
#include <petsc/finclude/petscsnes.h> #include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h> #include <petsc/finclude/petscdmda.h>
use DAMASK_interface use DAMASK_interface
use HDF5_utilities use HDF5_utilities
use PETScdmda use PETScdmda
use PETScsnes use PETScsnes
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_I3 math_I3
use spectral_utilities, only: & use spectral_utilities, only: &
tSolutionState, & tSolutionState, &
tSolutionParams tSolutionParams
implicit none implicit none
private private
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type(tSolutionParams), private :: params type(tSolutionParams), private :: params
type, private :: tNumerics type, private :: tNumerics
logical :: & logical :: update_gamma !< update gamma operator with current stiffness
update_gamma !< update gamma operator with current stiffness
end type tNumerics end type tNumerics
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM, private :: da DM, private :: da
SNES, private :: snes SNES, private :: snes
Vec, private :: solution_vec Vec, private :: solution_vec
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! common pointwise data ! common pointwise data
real(pReal), private, dimension(:,:,:,:,:), allocatable :: & real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, & !< field of previous compatible deformation gradients F_lastInc, & !< field of previous compatible deformation gradients
F_tau_lastInc, & !< field of previous incompatible deformation gradient F_tau_lastInc, & !< field of previous incompatible deformation gradient
Fdot, & !< field of assumed rate of compatible deformation gradient Fdot, & !< field of assumed rate of compatible deformation gradient
F_tauDot !< field of assumed rate of incopatible deformation gradient F_tauDot !< field of assumed rate of incopatible deformation gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc. ! stress, stiffness and compliance average etc.
real(pReal), private, dimension(3,3) :: & real(pReal), private, dimension(3,3) :: &
F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient
F_aim = math_I3, & !< current prescribed deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient
F_aim_lastInc = math_I3, & !< previous average deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient
F_av = 0.0_pReal, & !< average incompatible def grad field F_av = 0.0_pReal, & !< average incompatible def grad field
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character(len=1024), private :: incInfo !< time and increment information
real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness
S = 0.0_pReal, & !< current compliance (filled up with zeros)
C_scale = 0.0_pReal, &
S_scale = 0.0_pReal
real(pReal), private :: &
err_BC, & !< deviation from stress BC
err_curl, & !< RMS of curl of F
err_div !< RMS of div of P
integer, private :: & character(len=1024), private :: incInfo !< time and increment information
totalIter = 0 !< total iteration in current increment real(pReal), private, dimension(3,3,3,3) :: &
C_volAvg = 0.0_pReal, & !< current volume average stiffness
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness
C_minMaxAvgLastInc = 0.0_pReal, & !< previous (min+max)/2 stiffness
S = 0.0_pReal, & !< current compliance (filled up with zeros)
C_scale = 0.0_pReal, &
S_scale = 0.0_pReal
public :: & real(pReal), private :: &
grid_mech_spectral_polarisation_init, & err_BC, & !< deviation from stress BC
grid_mech_spectral_polarisation_solution, & err_curl, & !< RMS of curl of F
grid_mech_spectral_polarisation_forward err_div !< RMS of div of P
private :: &
converged, & integer, private :: &
formResidual totalIter = 0 !< total iteration in current increment
public :: &
grid_mech_spectral_polarisation_init, &
grid_mech_spectral_polarisation_solution, &
grid_mech_spectral_polarisation_forward
private :: &
converged, &
formResidual
contains contains
@ -88,141 +87,142 @@ contains
!> @brief allocates all necessary fields and fills them with data, potentially from restart info !> @brief allocates all necessary fields and fills them with data, potentially from restart info
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_init subroutine grid_mech_spectral_polarisation_init
use IO, only: & use IO, only: &
IO_intOut, & IO_intOut, &
IO_error, & IO_error, &
IO_open_jobFile_binary IO_open_jobFile_binary
use FEsolving, only: & use FEsolving, only: &
restartInc restartInc
use config, only :& use config, only :&
config_numerics config_numerics
use numerics, only: & use numerics, only: &
worldrank, & worldrank, &
worldsize, & worldsize, &
petsc_options petsc_options
use homogenization, only: & use homogenization, only: &
materialpoint_F0 materialpoint_F0
use DAMASK_interface, only: & use DAMASK_interface, only: &
getSolverJobName getSolverJobName
use spectral_utilities, only: & use spectral_utilities, only: &
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
utilities_updateGamma, & utilities_updateGamma, &
utilities_updateIPcoords utilities_updateIPcoords
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
use math, only: & use math, only: &
math_invSym3333 math_invSym3333
implicit none implicit none
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal temp33_Real = 0.0_pReal
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
FandF_tau, & ! overall pointer to solution data
F, & ! specific (sub)pointer
F_tau ! specific (sub)pointer
PetscInt, dimension(worldsize) :: localK
integer(HID_T) :: fileHandle
integer :: fileUnit
character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>' PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' FandF_tau, & ! overall pointer to solution data
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006' F, & ! specific (sub)pointer
F_tau ! specific (sub)pointer
PetscInt, dimension(worldsize) :: localK
integer(HID_T) :: fileHandle
integer :: fileUnit
character(len=1024) :: rankStr
write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'
num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0 write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:31–45, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
num%update_gamma = config_numerics%getInt('update_gamma',defaultVal=0) > 0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc ! set default and user defined options for PETSc
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr) call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate global fields ! allocate global fields
allocate (F_lastInc (3,3,grid(1),grid(2),grid3),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(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_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) allocate(F_tauDot (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 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)
localK = 0 localK = 0
localK(worldrank+1) = grid3 localK(worldrank+1) = grid3
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
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, &
18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap) 18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & ! local grid [grid(1)],[grid(2)],localK, & ! local grid
da,ierr) ! handle, error da,ierr) ! handle, error
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da
call DMsetFromOptions(da,ierr); CHKERRQ(ierr) call DMsetFromOptions(da,ierr); CHKERRQ(ierr)
call DMsetUp(da,ierr); CHKERRQ(ierr) call DMsetUp(da,ierr); CHKERRQ(ierr)
call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor) call DMcreateGlobalVector(da,solution_vec,ierr); CHKERRQ(ierr) ! global solution vector (grid x 18, i.e. every def grad tensor)
call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector call DMDASNESsetFunctionLocal(da,INSERT_VALUES,formResidual,PETSC_NULL_SNES,ierr) ! residual vector of same shape as solution vector
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged" call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged"
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
F => FandF_tau( 0: 8,:,:,:) F => FandF_tau( 0: 8,:,:,:)
F_tau => FandF_tau( 9:17,:,:,:) F_tau => FandF_tau( 9:17,:,:,:)
restart: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file'
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5')
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
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
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file'
call HDF5_read(fileHandle,C_volAvg, 'C_volAvg')
call HDF5_read(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeFile(fileHandle)
restart: if (restartInc > 0) then fileUnit = IO_open_jobFile_binary('C_ref')
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading values of increment ', restartInc, ' from file' read(fileUnit) C_minMaxAvg; close(fileUnit)
endif restartRead
write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5') call Utilities_updateGamma(C_minMaxAvg,.true.)
C_scale = C_minMaxAvg
call HDF5_read(fileHandle,F_aim, 'F_aim') S_scale = math_invSym3333(C_minMaxAvg)
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
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
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateIPcoords(reshape(F,shape(F_lastInc)))
call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
0.0_pReal, & ! time increment
math_I3) ! no rotation of boundary condition
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
restartRead: if (restartInc > 0) then
write(6,'(/,a,'//IO_intOut(restartInc)//',a)') ' reading more values of increment ', restartInc, ' from file'
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.)
C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg)
end subroutine grid_mech_spectral_polarisation_init end subroutine grid_mech_spectral_polarisation_init
@ -231,66 +231,65 @@ end subroutine grid_mech_spectral_polarisation_init
!> @brief solution for the Polarisation scheme with internal iterations !> @brief solution for the Polarisation scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
use math, only: & use math, only: &
math_invSym3333 math_invSym3333
use spectral_utilities, only: & use spectral_utilities, only: &
tBoundaryCondition, & tBoundaryCondition, &
utilities_maskedCompliance, & utilities_maskedCompliance, &
utilities_updateGamma utilities_updateGamma
use FEsolving, only: & use FEsolving, only: &
restartWrite, & restartWrite, &
terminallyIll terminallyIll
implicit none implicit none
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc, & !< time increment of current solution timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: & type(tBoundaryCondition), intent(in) :: &
stress_BC stress_BC
real(pReal), dimension(3,3), intent(in) :: rotation_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC
type(tSolutionState) :: & type(tSolutionState) :: &
solution solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc Data ! PETSc Data
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
incInfo = incInfoIn incInfo = incInfoIn
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator) ! 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 if (num%update_gamma) then
call utilities_updateGamma(C_minMaxAvg,restartWrite) call utilities_updateGamma(C_minMaxAvg,restartWrite)
C_scale = C_minMaxAvg C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg) S_scale = math_invSym3333(C_minMaxAvg)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide available data ! set module wide available data
params%stress_mask = stress_BC%maskFloat params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old params%timeincOld = timeinc_old
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve BVP ! solve BVP
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! check convergence ! check convergence
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
solution%converged = reason > 0 solution%converged = reason > 0
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
end function grid_mech_spectral_polarisation_solution end function grid_mech_spectral_polarisation_solution
@ -324,13 +323,12 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
use FEsolving, only: & use FEsolving, only: &
restartWrite restartWrite
implicit none
logical, intent(in) :: & logical, intent(in) :: &
guess guess
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc_old, & timeinc_old, &
timeinc, & timeinc, &
loadCaseTime !< remaining time of current load case loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: & type(tBoundaryCondition), intent(in) :: &
stress_BC, & stress_BC, &
deformation_BC deformation_BC
@ -349,30 +347,29 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
F_tau => FandF_tau( 9:17,:,:,:) F_tau => FandF_tau( 9:17,:,:,:)
if (cutBack) then if (cutBack) then
C_volAvg = C_volAvgLastInc ! QUESTION: where is this required? C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required? C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required?
else else
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! restart information for spectral solver ! restart information for spectral solver
if (restartWrite) then ! QUESTION: where is this logical properly set? if (restartWrite) then
write(6,'(/,a)') ' writing converged results for restart';flush(6) write(6,'(/,a)') ' writing converged results for restart';flush(6)
write(rankStr,'(a1,i0)')'_',worldrank write(rankStr,'(a1,i0)')'_',worldrank
fileHandle = HDF5_openFile(trim(getSolverJobName())//trim(rankStr)//'.hdf5','w') 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_volAvg, 'C_volAvg')
call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc') call HDF5_write(fileHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_write(fileHandle,F_aim, 'F_aim')
call HDF5_write(fileHandle,F_aimDot, 'F_aimDot')
call HDF5_write(fileHandle,F_aim_lastInc, 'F_aim_lastInc')
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_closeFile(fileHandle) call HDF5_closeFile(fileHandle)
endif endif
call CPFEM_age ! age state and kinematics call CPFEM_age ! age state and kinematics
@ -386,13 +383,13 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = & F_aimDot = &
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = & F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values F_aimDot + deformation_BC%maskFloat * deformation_BC%values
elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed
F_aimDot = & F_aimDot = &
F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime
endif endif
@ -404,33 +401,33 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
F_tauDot = utilities_calculateRate(guess, & F_tauDot = utilities_calculateRate(guess, &
F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, & F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
math_rotate_backward33(F_aimDot,rotation_BC)) math_rotate_backward33(F_aimDot,rotation_BC))
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) ! winding F forward
F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau forward F_tau_lastInc = reshape(F_tau, [3,3,grid(1),grid(2),grid3]) ! winding F_tau forward
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc F_aim = F_aim_lastInc + F_aimDot * timeinc
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
math_rotate_backward33(F_aim,rotation_BC)),& math_rotate_backward33(F_aim,rotation_BC)),&
[9,grid(1),grid(2),grid3]) [9,grid(1),grid(2),grid3])
if (guess) then if (guess) then
F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), & F_tau = reshape(Utilities_forwardField(timeinc,F_tau_lastInc,F_taudot), &
[9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition [9,grid(1),grid(2),grid3]) ! does not have any average value as boundary condition
else else
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, & F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, &
math_mul3333xx33(C_scale,& math_mul3333xx33(C_scale,&
matmul(transpose(F_lambda33),& matmul(transpose(F_lambda33),&
F_lambda33)-math_I3))*0.5_pReal)& F_lambda33)-math_I3))*0.5_pReal)&
+ math_I3 + math_I3
F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
endif endif
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_forward end subroutine grid_mech_spectral_polarisation_forward
@ -438,61 +435,61 @@ end subroutine grid_mech_spectral_polarisation_forward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convergence check !> @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: & use numerics, only: &
itmax, & itmax, &
itmin, & itmin, &
err_div_tolRel, & err_div_tolRel, &
err_div_tolAbs, & err_div_tolAbs, &
err_curl_tolRel, & err_curl_tolRel, &
err_curl_tolAbs, & err_curl_tolAbs, &
err_stress_tolRel, & err_stress_tolRel, &
err_stress_tolAbs err_stress_tolAbs
use FEsolving, only: & use FEsolving, only: &
terminallyIll terminallyIll
implicit none implicit none
SNES :: snes_local SNES :: snes_local
PetscInt :: PETScIter PetscInt, intent(in) :: PETScIter
PetscReal :: & PetscReal, intent(in) :: &
xnorm, & ! not used devNull1, &
snorm, & ! not used devNull2, &
fnorm ! not used devNull3
SNESConvergedReason :: reason SNESConvergedReason :: reason
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
real(pReal) :: & real(pReal) :: &
curlTol, & curlTol, &
divTol, & divTol, &
BCTol BCTol
curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs) curlTol = max(maxval(abs(F_aim-math_I3))*err_curl_tolRel ,err_curl_tolAbs)
divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs) divTol = max(maxval(abs(P_av)) *err_div_tolRel ,err_div_tolAbs)
BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs) BCTol = max(maxval(abs(P_av)) *err_stress_tolRel,err_stress_tolAbs)
if ((totalIter >= itmin .and. & if ((totalIter >= itmin .and. &
all([ err_div /divTol, & all([ err_div /divTol, &
err_curl/curlTol, & err_curl/curlTol, &
err_BC /BCTol ] < 1.0_pReal)) & err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then .or. terminallyIll) then
reason = 1 reason = 1
elseif (totalIter >= itmax) then elseif (totalIter >= itmax) then
reason = -1 reason = -1
else else
reason = 0 reason = 0
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report ! report
write(6,'(1/,a)') ' ... reporting .............................................................' write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', & write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')' err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', & write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')' err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', & write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', &
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')' err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)
end subroutine converged end subroutine converged
@ -502,144 +499,145 @@ end subroutine converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine formResidual(in, FandF_tau, & subroutine formResidual(in, FandF_tau, &
residuum, dummy,ierr) residuum, dummy,ierr)
use numerics, only: & use numerics, only: &
itmax, & itmax, &
itmin, & itmin, &
polarAlpha, & polarAlpha, &
polarBeta polarBeta
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
use math, only: & use math, only: &
math_rotate_forward33, & math_rotate_forward33, &
math_rotate_backward33, & math_rotate_backward33, &
math_mul3333xx33, & math_mul3333xx33, &
math_invSym3333 math_invSym3333
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_spectral, & debug_spectral, &
debug_spectralRotation debug_spectralRotation
use spectral_utilities, only: & use spectral_utilities, only: &
wgt, & wgt, &
tensorField_real, & tensorField_real, &
utilities_FFTtensorForward, & utilities_FFTtensorForward, &
utilities_fourierGammaConvolution, & utilities_fourierGammaConvolution, &
utilities_FFTtensorBackward, & utilities_FFTtensorBackward, &
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
utilities_divergenceRMS, & utilities_divergenceRMS, &
utilities_curlRMS utilities_curlRMS
use IO, only: & use IO, only: &
IO_intOut IO_intOut
use homogenization, only: & use homogenization, only: &
materialpoint_dPdF materialpoint_dPdF
use FEsolving, only: & use FEsolving, only: &
terminallyIll terminallyIll
implicit none implicit none
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work) 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), & PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), &
target, intent(in) :: FandF_tau target, intent(in) :: FandF_tau
PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),& PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),&
target, intent(out) :: residuum !< residuum field target, intent(out) :: residuum !< residuum field
PetscScalar, pointer, dimension(:,:,:,:,:) :: & PetscScalar, pointer, dimension(:,:,:,:,:) :: &
F, & F, &
F_tau, & F_tau, &
residual_F, & residual_F, &
residual_F_tau residual_F_tau
PetscInt :: & PetscInt :: &
PETScIter, & PETScIter, &
nfuncs nfuncs
PetscObject :: dummy PetscObject :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer :: & integer :: &
i, j, k, e i, j, k, e
F => FandF_tau(1:3,1:3,1,& F => FandF_tau(1:3,1:3,1,&
XG_RANGE,YG_RANGE,ZG_RANGE) XG_RANGE,YG_RANGE,ZG_RANGE)
F_tau => FandF_tau(1:3,1:3,2,& F_tau => FandF_tau(1:3,1:3,2,&
XG_RANGE,YG_RANGE,ZG_RANGE) XG_RANGE,YG_RANGE,ZG_RANGE)
residual_F => residuum(1:3,1:3,1,& residual_F => residuum(1:3,1:3,1,&
X_RANGE, Y_RANGE, Z_RANGE) X_RANGE, Y_RANGE, Z_RANGE)
residual_F_tau => residuum(1:3,1:3,2,& residual_F_tau => residuum(1:3,1:3,2,&
X_RANGE, Y_RANGE, Z_RANGE) X_RANGE, Y_RANGE, Z_RANGE)
F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr) call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! begin of new iteration ! begin of new iteration
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') & write(6,'(1x,a,3(a,'//IO_intOut(itmax)//'))') &
trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) & if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC)) ' deformation gradient aim (lab) =', transpose(math_rotate_backward33(F_aim,params%rotation_BC))
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim) ' deformation gradient aim =', transpose(F_aim)
flush(6) flush(6)
endif newIteration endif newIteration
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! !
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
tensorField_real(1:3,1:3,i,j,k) = & tensorField_real(1:3,1:3,i,j,k) = &
polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
polarAlpha*matmul(F(1:3,1:3,i,j,k), & polarAlpha*matmul(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)) math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
enddo; enddo; enddo enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing convolution in Fourier space ! doing convolution in Fourier space
call utilities_FFTtensorForward call utilities_FFTtensorForward
call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC)) call utilities_fourierGammaConvolution(math_rotate_backward33(polarBeta*F_aim,params%rotation_BC))
call utilities_FFTtensorBackward call utilities_FFTtensorBackward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, & P_av,C_volAvg,C_minMaxAvg, &
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC) F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc
err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim &
-math_rotate_forward33(F_av,params%rotation_BC)) + & -math_rotate_forward33(F_av,params%rotation_BC)) + &
params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc
! calculate divergence ! calculate divergence
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise
call utilities_FFTtensorForward call utilities_FFTtensorForward
err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual
e = 0 e = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
e = e + 1 e = e + 1
residual_F(1:3,1:3,i,j,k) = & 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), & 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) - matmul(F(1:3,1:3,i,j,k), & residual_F(1:3,1:3,i,j,k) - matmul(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))) & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
+ residual_F_tau(1:3,1:3,i,j,k) + residual_F_tau(1:3,1:3,i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating curl ! calculating curl
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward call utilities_FFTtensorForward
err_curl = Utilities_curlRMS() err_curl = Utilities_curlRMS()
end subroutine formResidual end subroutine formResidual