2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
2019-03-12 10:36:59 +05:30
!> @brief Grid solver for mechanics: Spectral basic
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-12 10:36:59 +05:30
module grid_mech_spectral_basic
2018-05-17 15:34:21 +05:30
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
use PETScdmda
use PETScsnes
2019-03-12 10:36:59 +05:30
use prec , only : &
2012-08-14 22:28:23 +05:30
pInt , &
pReal
use math , only : &
math_I3
2015-12-16 02:15:54 +05:30
use spectral_utilities , only : &
2015-06-03 23:00:31 +05:30
tSolutionState , &
tSolutionParams
2012-10-02 20:56:56 +05:30
2012-08-14 22:28:23 +05:30
implicit none
2012-11-09 03:03:58 +05:30
private
2014-09-11 18:58:15 +05:30
2012-08-14 22:28:23 +05:30
character ( len = * ) , parameter , public :: &
2019-03-12 10:36:59 +05:30
GRID_MECH_SPECTRAL_BASIC_LABEL = 'spectral_basic'
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
! derived types
2012-10-02 20:56:56 +05:30
type ( tSolutionParams ) , private :: params
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
! PETSc data
2012-11-12 19:44:39 +05:30
DM , private :: da
SNES , private :: snes
Vec , private :: solution_vec
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
! common pointwise data
2015-03-25 21:36:19 +05:30
real ( pReal ) , private , dimension ( : , : , : , : , : ) , allocatable :: F_lastInc , Fdot
2012-11-12 19:44:39 +05:30
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
2012-11-12 19:44:39 +05:30
real ( pReal ) , private , dimension ( 3 , 3 ) :: &
2018-02-16 20:06:18 +05:30
F_aimDot = 0.0_pReal , & !< assumed rate of average deformation gradient
F_aim = math_I3 , & !< current prescribed deformation gradient
F_aim_lastInc = math_I3 , & !< previous average deformation gradient
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
character ( len = 1024 ) , private :: incInfo !< time and increment information
2012-11-12 19:44:39 +05:30
real ( pReal ) , private , dimension ( 3 , 3 , 3 , 3 ) :: &
2013-08-08 14:43:29 +05:30
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
2018-02-16 20:06:18 +05:30
C_minMaxAvgLastInc = 0.0_pReal , & !< previous (min+max)/2 stiffness
S = 0.0_pReal !< current compliance (filled up with zeros)
real ( pReal ) , private :: &
err_BC , & !< deviation from stress BC
err_div !< RMS of div of P
2013-08-08 14:43:29 +05:30
integer ( pInt ) , private :: &
totalIter = 0_pInt !< total iteration in current increment
2018-02-16 20:06:18 +05:30
2013-01-10 03:49:32 +05:30
public :: &
2019-03-12 10:36:59 +05:30
grid_mech_spectral_basic_init , &
grid_mech_spectral_basic_solution , &
grid_mech_spectral_basic_forward
2013-01-10 03:49:32 +05:30
contains
2012-11-09 03:03:58 +05:30
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
2018-02-16 20:06:18 +05:30
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-12 10:36:59 +05:30
subroutine grid_mech_spectral_basic_init
2012-10-02 20:56:56 +05:30
use IO , only : &
2013-08-08 14:43:29 +05:30
IO_intOut , &
2018-11-17 21:20:19 +05:30
IO_error , &
2019-03-09 03:46:56 +05:30
IO_open_jobFile_binary
2013-05-08 21:22:29 +05:30
use debug , only : &
2015-09-11 14:22:03 +05:30
debug_level , &
debug_spectral , &
debug_spectralRestart
2012-10-02 20:56:56 +05:30
use FEsolving , only : &
restartInc
2015-03-25 21:36:19 +05:30
use numerics , only : &
worldrank , &
2019-03-12 10:36:59 +05:30
worldsize , &
petsc_options
2018-02-16 20:06:18 +05:30
use homogenization , only : &
materialpoint_F0
2012-10-02 20:56:56 +05:30
use DAMASK_interface , only : &
getSolverJobName
2015-12-16 02:15:54 +05:30
use spectral_utilities , only : &
2019-03-12 10:36:59 +05:30
utilities_constitutiveResponse , &
utilities_updateGamma , &
utilities_updateIPcoords , &
2015-03-25 21:36:19 +05:30
wgt
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid , &
grid3
2012-10-02 20:56:56 +05:30
use math , only : &
math_invSym3333
2013-08-08 14:43:29 +05:30
2012-10-02 20:56:56 +05:30
implicit none
2015-09-12 23:56:25 +05:30
real ( pReal ) , dimension ( 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ) :: P
2012-12-14 20:48:04 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: &
temp33_Real = 0.0_pReal
2016-06-27 21:20:43 +05:30
PetscErrorCode :: ierr
PetscScalar , pointer , dimension ( : , : , : , : ) :: F
2019-03-09 14:24:33 +05:30
PetscInt , dimension ( worldsize ) :: localK
integer :: fileUnit
2015-03-25 21:36:19 +05:30
character ( len = 1024 ) :: rankStr
2019-03-12 10:36:59 +05:30
write ( 6 , '(/,a)' ) ' <<<+- grid_mech_spectral_basic init -+>>>'
2019-03-09 15:32:12 +05:30
write ( 6 , '(/,a)' ) ' Eisenlohr et al., International Journal of Plasticity 46:37– 53, 2013'
write ( 6 , '(a)' ) ' https://doi.org/10.1016/j.ijplas.2012.09.012'
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'
2013-12-18 14:39:32 +05:30
2019-03-12 10:36:59 +05:30
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
call PETScOptionsInsertString ( PETSC_NULL_OPTIONS , '-mech_snes_type ngmres' , ierr )
CHKERRQ ( ierr )
call PETScOptionsInsertString ( PETSC_NULL_OPTIONS , trim ( petsc_options ) , ierr )
CHKERRQ ( ierr )
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! allocate global fields
2018-02-16 20:06:18 +05:30
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 )
2012-08-14 22:28:23 +05:30
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
2013-01-10 03:49:32 +05:30
call SNESCreate ( PETSC_COMM_WORLD , snes , ierr ) ; CHKERRQ ( ierr )
2015-06-03 23:00:31 +05:30
call SNESSetOptionsPrefix ( snes , 'mech_' , ierr ) ; CHKERRQ ( ierr )
2019-03-09 14:24:33 +05:30
localK = 0
localK ( worldrank + 1 ) = grid3
call MPI_Allreduce ( MPI_IN_PLACE , localK , worldsize , MPI_INTEGER , MPI_SUM , PETSC_COMM_WORLD , ierr )
2012-11-06 21:30:51 +05:30
call DMDACreate3d ( PETSC_COMM_WORLD , &
2014-11-19 11:05:10 +05:30
DM_BOUNDARY_NONE , DM_BOUNDARY_NONE , DM_BOUNDARY_NONE , & ! cut off stencil at boundary
DMDA_STENCIL_BOX , & ! Moore (26) neighborhood around central point
2015-12-14 23:42:09 +05:30
grid ( 1 ) , grid ( 2 ) , grid ( 3 ) , & ! global grid
2016-06-27 21:20:43 +05:30
1 , 1 , worldsize , &
2015-03-25 21:36:19 +05:30
9 , 0 , & ! #dof (F tensor), ghost boundary width (domain overlap)
2018-05-17 15:34:21 +05:30
[ grid ( 1 ) ] , [ grid ( 2 ) ] , localK , & ! local grid
2014-11-19 11:05:10 +05:30
da , ierr ) ! handle, error
2015-03-25 21:36:19 +05:30
CHKERRQ ( ierr )
2018-02-16 20:06:18 +05:30
call SNESSetDM ( snes , da , ierr ) ; CHKERRQ ( ierr ) ! connect snes to da
2018-05-17 15:34:21 +05:30
call DMsetFromOptions ( da , ierr ) ; CHKERRQ ( ierr )
call DMsetUp ( da , ierr ) ; CHKERRQ ( ierr )
call DMcreateGlobalVector ( da , solution_vec , ierr ) ; CHKERRQ ( ierr ) ! global solution vector (grid x 9, i.e. every def grad tensor)
2019-03-12 10:36:59 +05:30
call DMDASNESsetFunctionLocal ( da , INSERT_VALUES , grid_mech_spectral_basic_formResidual , PETSC_NULL_SNES , ierr ) ! residual vector of same shape as solution vector
2013-12-17 21:07:14 +05:30
CHKERRQ ( ierr )
2019-03-12 10:36:59 +05:30
call SNESsetConvergenceTest ( snes , grid_mech_spectral_basic_converged , PETSC_NULL_SNES , PETSC_NULL_FUNCTION , ierr ) ! specify custom convergence check function "_converged"
2013-11-19 21:36:53 +05:30
CHKERRQ ( ierr )
2018-05-17 15:34:21 +05:30
call SNESsetFromOptions ( snes , ierr ) ; CHKERRQ ( ierr ) ! pull it all together with additional CLI arguments
2012-08-14 22:28:23 +05:30
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! init fields
2013-12-18 14:39:32 +05:30
call DMDAVecGetArrayF90 ( da , solution_vec , F , ierr ) ; CHKERRQ ( ierr ) ! get the data out of PETSc to work with
2018-02-17 03:11:07 +05:30
restart : if ( restartInc > 0_pInt ) then
2018-02-16 20:06:18 +05:30
if ( iand ( debug_level ( debug_spectral ) , debug_spectralRestart ) / = 0 ) then
2018-02-17 03:11:07 +05:30
write ( 6 , '(/,a,' / / IO_intOut ( restartInc ) / / ',a)' ) &
'reading values of increment ' , restartInc , ' from file'
2018-02-16 20:06:18 +05:30
flush ( 6 )
endif
2019-03-09 03:46:56 +05:30
fileUnit = IO_open_jobFile_binary ( 'F_aimDot' )
read ( fileUnit ) F_aimDot ; close ( fileUnit )
2015-03-25 21:36:19 +05:30
write ( rankStr , '(a1,i0)' ) '_' , worldrank
2019-03-09 03:46:56 +05:30
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 )
2014-03-25 21:14:16 +05:30
F_aim = reshape ( sum ( sum ( sum ( F , dim = 4 ) , dim = 3 ) , dim = 2 ) * wgt , [ 3 , 3 ] ) ! average of F
2018-11-17 21:20:19 +05:30
call MPI_Allreduce ( MPI_IN_PLACE , F_aim , 9 , MPI_DOUBLE , MPI_SUM , PETSC_COMM_WORLD , ierr )
if ( ierr / = 0_pInt ) call IO_error ( 894_pInt , ext_msg = 'F_aim' )
2014-03-25 21:14:16 +05:30
F_aim_lastInc = sum ( sum ( sum ( F_lastInc , dim = 5 ) , dim = 4 ) , dim = 3 ) * wgt ! average of F_lastInc
2018-11-17 21:20:19 +05:30
call MPI_Allreduce ( MPI_IN_PLACE , F_aim_lastInc , 9 , MPI_DOUBLE , MPI_SUM , PETSC_COMM_WORLD , ierr )
if ( ierr / = 0_pInt ) call IO_error ( 894_pInt , ext_msg = 'F_aim_lastInc' )
2018-02-17 03:11:07 +05:30
elseif ( restartInc == 0_pInt ) then restart
2015-09-11 14:22:03 +05:30
F_lastInc = spread ( spread ( spread ( math_I3 , 3 , grid ( 1 ) ) , 4 , grid ( 2 ) ) , 5 , grid3 ) ! initialize to identity
F = reshape ( F_lastInc , [ 9 , grid ( 1 ) , grid ( 2 ) , grid3 ] )
endif restart
2014-02-03 21:27:04 +05:30
2018-02-16 20:06:18 +05:30
materialpoint_F0 = reshape ( F_lastInc , [ 3 , 3 , 1 , product ( grid ( 1 : 2 ) ) * grid3 ] ) ! set starting condition for materialpoint_stressAndItsTangent
2015-09-24 23:08:49 +05:30
call Utilities_updateIPcoords ( reshape ( F , shape ( F_lastInc ) ) )
2018-02-16 20:06:18 +05:30
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
2014-03-25 21:14:16 +05:30
call DMDAVecRestoreArrayF90 ( da , solution_vec , F , ierr ) ; CHKERRQ ( ierr ) ! write data back to PETSc
2018-02-16 20:06:18 +05:30
! QUESTION: why not writing back right after reading (l.189)?
2012-12-14 01:50:04 +05:30
2019-03-12 10:36:59 +05:30
restartRead : if ( restartInc > 0_pInt ) then
if ( iand ( debug_level ( debug_spectral ) , debug_spectralRestart ) / = 0 ) &
2018-02-17 03:11:07 +05:30
write ( 6 , '(/,a,' / / IO_intOut ( restartInc ) / / ',a)' ) &
'reading more values of increment ' , restartInc , ' from file'
2014-02-03 21:27:04 +05:30
flush ( 6 )
2019-03-09 03:46:56 +05:30
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 )
fileUnit = IO_open_jobFile_binary ( 'C_ref' )
read ( fileUnit ) C_minMaxAvg ; close ( fileUnit )
2015-09-11 14:22:03 +05:30
endif restartRead
2018-02-16 20:06:18 +05:30
call Utilities_updateGamma ( C_minMaxAvg , . true . )
2012-11-12 19:44:39 +05:30
2019-03-12 10:36:59 +05:30
end subroutine grid_mech_spectral_basic_init
2018-02-16 20:06:18 +05:30
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
2018-05-17 19:57:36 +05:30
!> @brief solution for the Basic scheme with internal iterations
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-12 10:36:59 +05:30
function grid_mech_spectral_basic_solution ( incInfoIn , timeinc , timeinc_old , stress_BC , rotation_BC ) result ( solution )
2012-11-12 19:44:39 +05:30
use numerics , only : &
2015-12-15 01:34:59 +05:30
update_gamma
2015-12-16 02:15:54 +05:30
use spectral_utilities , only : &
2012-11-12 19:44:39 +05:30
tBoundaryCondition , &
2019-03-12 10:36:59 +05:30
utilities_maskedCompliance , &
utilities_updateGamma
2012-11-12 19:44:39 +05:30
use FEsolving , only : &
restartWrite , &
terminallyIll
2013-01-10 03:49:32 +05:30
2012-11-12 19:44:39 +05:30
implicit none
2013-08-08 14:43:29 +05:30
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
! input data for solution
2018-07-04 00:56:11 +05:30
character ( len = * ) , intent ( in ) :: &
2018-02-16 20:06:18 +05:30
incInfoIn
2018-07-04 00:56:11 +05:30
real ( pReal ) , intent ( in ) :: &
2019-03-12 10:36:59 +05:30
timeinc , & !< time increment of current solution
timeinc_old !< time increment of last successful increment
2018-07-04 00:56:11 +05:30
type ( tBoundaryCondition ) , intent ( in ) :: &
2016-10-26 02:24:32 +05:30
stress_BC
2012-11-12 19:44:39 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: rotation_BC
2019-03-12 10:36:59 +05:30
type ( tSolutionState ) :: &
solution
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
2013-08-08 14:43:29 +05:30
! PETSc Data
2019-03-12 10:36:59 +05:30
PetscErrorCode :: ierr
2012-11-12 19:44:39 +05:30
SNESConvergedReason :: reason
2016-06-27 21:20:43 +05:30
2012-11-12 19:44:39 +05:30
incInfo = incInfoIn
2012-12-14 20:48:04 +05:30
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
2016-10-26 02:24:32 +05:30
S = Utilities_maskedCompliance ( rotation_BC , stress_BC % maskLogical , C_volAvg )
2018-02-16 20:06:18 +05:30
if ( update_gamma ) call Utilities_updateGamma ( C_minMaxAvg , restartWrite )
2012-11-12 19:44:39 +05:30
2013-01-10 19:03:43 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-12 10:36:59 +05:30
! set module wide available data
2018-07-04 00:56:11 +05:30
params % stress_mask = stress_BC % maskFloat
2016-10-26 02:24:32 +05:30
params % stress_BC = stress_BC % values
2012-11-12 19:44:39 +05:30
params % rotation_BC = rotation_BC
2016-10-26 02:24:32 +05:30
params % timeinc = timeinc
params % timeincOld = timeinc_old
2012-11-12 19:44:39 +05:30
2015-12-02 04:06:19 +05:30
!--------------------------------------------------------------------------------------------------
! solve BVP
2018-05-17 15:34:21 +05:30
call SNESsolve ( snes , PETSC_NULL_VEC , solution_vec , ierr ) ; CHKERRQ ( ierr )
2015-12-02 04:06:19 +05:30
!--------------------------------------------------------------------------------------------------
! check convergence
2018-02-16 20:06:18 +05:30
call SNESGetConvergedReason ( snes , reason , ierr ) ; CHKERRQ ( ierr )
2019-03-12 10:36:59 +05:30
solution % converged = reason > 0
solution % iterationsNeeded = totalIter
solution % termIll = terminallyIll
2012-11-12 19:44:39 +05:30
terminallyIll = . false .
2012-08-14 22:28:23 +05:30
2019-03-12 10:36:59 +05:30
end function grid_mech_spectral_basic_solution
2012-08-14 22:28:23 +05:30
2013-01-10 03:49:32 +05:30
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
2018-05-17 19:57:36 +05:30
!> @brief forms the basic residual vector
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-12 10:36:59 +05:30
subroutine grid_mech_spectral_basic_formResidual ( in , & ! DMDA info (needs to be named "in" for XRANGE, etc. macros to work)
2018-07-04 00:56:11 +05:30
F , & ! defgrad field on grid
residuum , & ! residuum field on grid
dummy , &
ierr )
2012-11-12 19:44:39 +05:30
use numerics , only : &
itmax , &
itmin
2015-03-25 21:36:19 +05:30
use mesh , only : &
2015-09-11 14:22:03 +05:30
grid , &
grid3
2012-11-12 19:44:39 +05:30
use math , only : &
math_rotate_backward33 , &
math_mul3333xx33
2013-05-08 21:22:29 +05:30
use debug , only : &
debug_level , &
debug_spectral , &
debug_spectralRotation
2015-12-16 02:15:54 +05:30
use spectral_utilities , only : &
2015-09-11 14:22:03 +05:30
tensorField_real , &
2015-06-03 23:00:31 +05:30
utilities_FFTtensorForward , &
utilities_fourierGammaConvolution , &
2016-06-27 21:20:43 +05:30
utilities_FFTtensorBackward , &
2019-03-12 10:36:59 +05:30
utilities_constitutiveResponse , &
utilities_divergenceRMS
2013-05-08 21:22:29 +05:30
use IO , only : &
IO_intOut
2015-03-25 21:36:19 +05:30
use FEsolving , only : &
terminallyIll
2012-11-12 19:44:39 +05:30
implicit none
2018-02-16 20:06:18 +05:30
DMDALocalInfo , dimension ( DMDA_LOCAL_INFO_SIZE ) :: in
PetscScalar , &
2018-07-04 00:56:11 +05:30
dimension ( 3 , 3 , XG_RANGE , YG_RANGE , ZG_RANGE ) , intent ( in ) :: F
2018-02-16 20:06:18 +05:30
PetscScalar , &
2018-07-04 00:56:11 +05:30
dimension ( 3 , 3 , X_RANGE , Y_RANGE , Z_RANGE ) , intent ( out ) :: residuum
2019-03-12 10:36:59 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: &
deltaF_aim
2013-01-10 19:03:43 +05:30
PetscInt :: &
2013-08-08 14:43:29 +05:30
PETScIter , &
2013-01-10 19:03:43 +05:30
nfuncs
2012-11-12 19:44:39 +05:30
PetscObject :: dummy
PetscErrorCode :: ierr
2013-01-10 03:49:32 +05:30
call SNESGetNumberFunctionEvals ( snes , nfuncs , ierr ) ; CHKERRQ ( ierr )
2013-08-08 14:43:29 +05:30
call SNESGetIterationNumber ( snes , PETScIter , ierr ) ; CHKERRQ ( ierr )
2012-12-16 05:22:06 +05:30
2018-02-16 20:06:18 +05:30
if ( nfuncs == 0 . and . PETScIter == 0 ) totalIter = - 1_pInt ! new increment
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
2018-02-16 20:06:18 +05:30
! begin of new iteration
newIteration : if ( totalIter < = PETScIter ) then
2013-08-08 14:43:29 +05:30
totalIter = totalIter + 1_pInt
2018-02-16 20:06:18 +05:30
write ( 6 , '(1x,a,3(a,' / / IO_intOut ( itmax ) / / '))' ) &
trim ( incInfo ) , ' @ Iteration ' , itmin , '≤' , totalIter , '≤' , itmax
2016-10-26 02:24:32 +05:30
if ( iand ( debug_level ( debug_spectral ) , debug_spectralRotation ) / = 0 ) &
2018-02-16 20:06:18 +05:30
write ( 6 , '(/,a,/,3(3(f12.7,1x)/))' , advance = 'no' ) &
2018-02-25 17:18:58 +05:30
' deformation gradient aim (lab) =' , transpose ( math_rotate_backward33 ( F_aim , params % rotation_BC ) )
2018-02-16 20:06:18 +05:30
write ( 6 , '(/,a,/,3(3(f12.7,1x)/))' , advance = 'no' ) &
2018-02-25 17:18:58 +05:30
' deformation gradient aim =' , transpose ( F_aim )
2016-10-26 02:24:32 +05:30
flush ( 6 )
2015-09-11 14:22:03 +05:30
endif newIteration
2013-07-26 21:55:37 +05:30
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
2018-07-04 00:56:11 +05:30
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 )
2015-03-25 21:36:19 +05:30
call MPI_Allreduce ( MPI_IN_PLACE , terminallyIll , 1 , MPI_LOGICAL , MPI_LOR , PETSC_COMM_WORLD , ierr )
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
! stress BC handling
2018-02-16 20:06:18 +05:30
deltaF_aim = math_mul3333xx33 ( S , P_av - params % stress_BC )
F_aim = F_aim - deltaF_aim
2018-07-04 00:56:11 +05:30
err_BC = maxval ( abs ( params % stress_mask * ( P_av - params % stress_BC ) ) ) ! mask = 0.0 when no stress bc
2015-12-15 01:34:59 +05:30
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme
2015-09-11 14:22:03 +05:30
tensorField_real = 0.0_pReal
2018-07-04 00:56:11 +05:30
tensorField_real ( 1 : 3 , 1 : 3 , 1 : grid ( 1 ) , 1 : grid ( 2 ) , 1 : grid3 ) = residuum ! store fPK field for subsequent FFT forward transform
2018-02-16 20:06:18 +05:30
call utilities_FFTtensorForward ( ) ! FFT forward of global "tensorField_real"
2018-07-04 00:56:11 +05:30
err_div = Utilities_divergenceRMS ( ) ! divRMS of tensorField_fourier for later use
2018-02-16 20:06:18 +05:30
call utilities_fourierGammaConvolution ( math_rotate_backward33 ( deltaF_aim , params % rotation_BC ) ) ! convolution of Gamma and tensorField_fourier, with arg
call utilities_FFTtensorBackward ( ) ! FFT backward of global tensorField_fourier
2012-11-12 19:44:39 +05:30
!--------------------------------------------------------------------------------------------------
2013-08-08 14:43:29 +05:30
! constructing residual
2018-07-04 00:56:11 +05:30
residuum = tensorField_real ( 1 : 3 , 1 : 3 , 1 : grid ( 1 ) , 1 : grid ( 2 ) , 1 : grid3 ) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too
2013-12-18 14:39:32 +05:30
2019-03-12 10:36:59 +05:30
end subroutine grid_mech_spectral_basic_formResidual
2012-11-12 19:44:39 +05:30
2012-08-14 22:28:23 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
2019-03-12 10:36:59 +05:30
subroutine grid_mech_spectral_basic_converged ( snes_local , PETScIter , xnorm , snorm , fnorm , reason , dummy , ierr )
2012-11-12 19:44:39 +05:30
use numerics , only : &
itmax , &
itmin , &
2013-08-08 14:43:29 +05:30
err_div_tolRel , &
err_div_tolAbs , &
2013-08-07 22:50:05 +05:30
err_stress_tolRel , &
2016-10-26 02:24:32 +05:30
err_stress_tolAbs
2013-03-06 20:01:13 +05:30
use FEsolving , only : &
terminallyIll
2016-06-27 21:20:43 +05:30
2012-11-12 19:44:39 +05:30
implicit none
SNES :: snes_local
2013-08-08 14:43:29 +05:30
PetscInt :: PETScIter
2013-01-10 03:49:32 +05:30
PetscReal :: &
2018-07-04 00:56:11 +05:30
xnorm , & ! not used
snorm , & ! not used
fnorm ! not used
2012-11-12 19:44:39 +05:30
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
2013-01-10 19:03:43 +05:30
real ( pReal ) :: &
2013-08-07 22:50:05 +05:30
divTol , &
2018-02-16 20:06:18 +05:30
BCTol
2016-06-27 21:20:43 +05:30
2018-02-16 20:06:18 +05:30
divTol = max ( maxval ( abs ( P_av ) ) * err_div_tolRel , err_div_tolAbs )
BCTol = max ( maxval ( abs ( P_av ) ) * err_stress_tolRel , err_stress_tolAbs )
2016-06-27 21:20:43 +05:30
2013-10-23 20:52:12 +05:30
converged : if ( ( totalIter > = itmin . and . &
2013-12-18 14:39:32 +05:30
all ( [ err_div / divTol , &
2018-02-16 20:06:18 +05:30
err_BC / BCTol ] < 1.0_pReal ) ) &
2013-08-08 14:43:29 +05:30
. or . terminallyIll ) then
2012-11-12 19:44:39 +05:30
reason = 1
2013-08-07 22:50:05 +05:30
elseif ( totalIter > = itmax ) then converged
2012-11-12 19:44:39 +05:30
reason = - 1
2013-08-07 22:50:05 +05:30
else converged
2012-11-12 19:44:39 +05:30
reason = 0
2013-08-07 22:50:05 +05:30
endif converged
!--------------------------------------------------------------------------------------------------
! report
2016-10-26 02:24:32 +05:30
write ( 6 , '(1/,a)' ) ' ... reporting .............................................................'
write ( 6 , '(1/,a,f12.2,a,es8.2,a,es9.2,a)' ) ' error divergence = ' , &
2018-02-16 20:06:18 +05:30
err_div / divTol , ' (' , err_div , ' / m, tol = ' , divTol , ')'
write ( 6 , '(a,f12.2,a,es8.2,a,es9.2,a)' ) ' error stress BC = ' , &
err_BC / BCTol , ' (' , err_BC , ' Pa, tol = ' , BCTol , ')'
2016-10-26 02:24:32 +05:30
write ( 6 , '(/,a)' ) ' ==========================================================================='
2019-03-12 10:36:59 +05:30
flush ( 6 )
2013-01-10 03:49:32 +05:30
2019-03-12 10:36:59 +05:30
end subroutine grid_mech_spectral_basic_converged
2012-08-14 22:28:23 +05:30
2013-12-18 15:05:05 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
2018-02-16 20:06:18 +05:30
!> @details find new boundary conditions and best F estimate for end of current timestep
!> possibly writing restart information, triggering of state increment in DAMASK, and updating of IPcoordinates
2013-12-18 15:05:05 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-12 10:36:59 +05:30
subroutine grid_mech_spectral_basic_forward ( guess , timeinc , timeinc_old , loadCaseTime , deformation_BC , stress_BC , rotation_BC )
2018-02-16 20:06:18 +05:30
use math , only : &
math_mul33x33 , &
math_rotate_backward33
use numerics , only : &
2019-03-12 10:36:59 +05:30
worldrank
2018-02-16 20:06:18 +05:30
use homogenization , only : &
materialpoint_F0
use mesh , only : &
grid , &
grid3
use CPFEM2 , only : &
CPFEM_age
use spectral_utilities , only : &
2019-03-12 10:36:59 +05:30
utilities_calculateRate , &
utilities_forwardField , &
utilities_updateIPcoords , &
2018-02-16 20:06:18 +05:30
tBoundaryCondition , &
cutBack
use IO , only : &
2019-03-09 03:46:56 +05:30
IO_open_jobFile_binary
2018-02-16 20:06:18 +05:30
use FEsolving , only : &
restartWrite
implicit none
2018-07-04 00:56:11 +05:30
logical , intent ( in ) :: &
2018-02-16 20:06:18 +05:30
guess
2018-07-04 00:56:11 +05:30
real ( pReal ) , intent ( in ) :: &
2018-02-16 20:06:18 +05:30
timeinc_old , &
timeinc , &
loadCaseTime !< remaining time of current load case
type ( tBoundaryCondition ) , intent ( in ) :: &
stress_BC , &
deformation_BC
2018-07-04 00:56:11 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
2018-02-16 20:06:18 +05:30
rotation_BC
2019-03-12 10:36:59 +05:30
PetscErrorCode :: ierr
2018-02-16 20:06:18 +05:30
PetscScalar , dimension ( : , : , : , : ) , pointer :: F
2019-03-09 03:46:56 +05:30
integer :: fileUnit
2018-02-16 20:06:18 +05:30
character ( len = 32 ) :: rankStr
call DMDAVecGetArrayF90 ( da , solution_vec , F , ierr ) ; CHKERRQ ( ierr )
if ( cutBack ) then
C_volAvg = C_volAvgLastInc ! QUESTION: where is this required?
C_minMaxAvg = C_minMaxAvgLastInc ! QUESTION: where is this required?
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 )
2019-03-09 03:46:56 +05:30
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_aimDot' , 'w' )
write ( fileUnit ) F_aimDot ; close ( fileUnit )
2018-02-16 20:06:18 +05:30
endif
write ( rankStr , '(a1,i0)' ) '_' , worldrank
2019-03-09 03:46:56 +05:30
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 )
2018-02-16 20:06:18 +05:30
endif
call CPFEM_age ( ) ! age state and kinematics
call utilities_updateIPcoords ( F )
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
2018-02-25 17:18:58 +05:30
F_aimDot = merge ( stress_BC % maskFloat * ( F_aim - F_aim_lastInc ) / timeinc_old , 0.0_pReal , guess )
2018-02-16 20:06:18 +05:30
F_aim_lastInc = F_aim
2018-02-25 17:18:58 +05:30
2018-02-16 20:06:18 +05:30
!--------------------------------------------------------------------------------------------------
! calculate rate for aim
if ( deformation_BC % myType == 'l' ) then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC % maskFloat * math_mul33x33 ( deformation_BC % values , F_aim_lastInc )
elseif ( deformation_BC % myType == 'fdot' ) then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC % maskFloat * deformation_BC % values
elseif ( deformation_BC % myType == 'f' ) then ! aim at end of load case is prescribed
F_aimDot = &
F_aimDot + deformation_BC % maskFloat * ( deformation_BC % values - F_aim_lastInc ) / loadCaseTime
endif
Fdot = Utilities_calculateRate ( guess , &
F_lastInc , reshape ( F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ] ) , timeinc_old , &
math_rotate_backward33 ( F_aimDot , rotation_BC ) )
F_lastInc = reshape ( F , [ 3 , 3 , grid ( 1 ) , grid ( 2 ) , grid3 ] ) ! winding F forward
materialpoint_F0 = reshape ( F_lastInc , [ 3 , 3 , 1 , product ( grid ( 1 : 2 ) ) * grid3 ] ) ! set starting condition for materialpoint_stressAndItsTangent
endif
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
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
math_rotate_backward33 ( F_aim , rotation_BC ) ) , [ 9 , grid ( 1 ) , grid ( 2 ) , grid3 ] )
call DMDAVecRestoreArrayF90 ( da , solution_vec , F , ierr ) ; CHKERRQ ( ierr )
2019-03-12 10:36:59 +05:30
end subroutine grid_mech_spectral_basic_forward
2013-10-23 20:52:12 +05:30
2019-03-12 10:36:59 +05:30
end module grid_mech_spectral_basic