DAMASK_EICMD/src/grid/grid_mech_spectral_polarisa...

651 lines
34 KiB
Fortran
Raw Normal View History

2013-07-24 18:36:16 +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
!> @brief Grid solver for mechanics: Spectral Polarisation
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-12 23:26:28 +05:30
module grid_mech_spectral_polarisation
#include <petsc/finclude/petscsnes.h>
#include <petsc/finclude/petscdmda.h>
2019-04-15 19:23:46 +05:30
use PETScdmda
use PETScsnes
2019-06-10 00:50:38 +05:30
use prec
use DAMASK_interface
use HDF5_utilities
use math
2019-12-03 00:36:58 +05:30
use rotations
2019-06-10 00:50:38 +05:30
use spectral_utilities
use FEsolving
use config
use numerics
use homogenization
2020-03-20 11:28:55 +05:30
use discretization_grid
2019-06-10 00:50:38 +05:30
use debug
2020-06-16 22:17:19 +05:30
use YAML_types
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
implicit none
private
2020-02-25 23:23:12 +05:30
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
2016-06-27 21:20:43 +05:30
! derived types
type(tSolutionParams), private :: params
2020-02-25 23:23:12 +05:30
type, private :: tNumerics
2019-04-15 19:23:46 +05:30
logical :: update_gamma !< update gamma operator with current stiffness
end type tNumerics
2020-02-25 23:23:12 +05:30
type(tNumerics) :: num ! numerics parameters. Better name?
2020-02-25 23:23:12 +05:30
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! PETSc data
2019-04-15 19:23:46 +05:30
DM, private :: da
SNES, private :: snes
Vec, private :: solution_vec
2016-06-27 21:20:43 +05:30
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! common pointwise data
2019-04-15 19:23:46 +05:30
real(pReal), private, dimension(:,:,:,:,:), allocatable :: &
F_lastInc, & !< field of previous compatible deformation gradients
2020-02-25 23:23:12 +05:30
F_tau_lastInc, & !< field of previous incompatible deformation gradient
2019-04-15 19:23:46 +05:30
Fdot, & !< field of assumed rate of compatible deformation gradient
F_tauDot !< field of assumed rate of incopatible deformation gradient
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
2019-04-15 19:23:46 +05:30
real(pReal), private, dimension(3,3) :: &
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
F_av = 0.0_pReal, & !< average incompatible def grad field
P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress
2020-02-25 23:23:12 +05:30
2020-01-04 23:31:36 +05:30
character(len=pStringLen), private :: incInfo !< time and increment information
2019-04-15 19:23:46 +05:30
real(pReal), private, dimension(3,3,3,3) :: &
2020-02-25 23:23:12 +05:30
C_volAvg = 0.0_pReal, & !< current volume average stiffness
2019-04-15 19:23:46 +05:30
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
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
real(pReal), private :: &
err_BC, & !< deviation from stress BC
err_curl, & !< RMS of curl of F
err_div !< RMS of div of P
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
integer, private :: &
totalIter = 0 !< total iteration in current increment
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
public :: &
grid_mech_spectral_polarisation_init, &
grid_mech_spectral_polarisation_solution, &
2019-10-24 17:16:36 +05:30
grid_mech_spectral_polarisation_forward, &
grid_mech_spectral_polarisation_updateCoords, &
2019-10-24 17:16:36 +05:30
grid_mech_spectral_polarisation_restartWrite
2013-07-24 18:36:16 +05:30
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
2019-03-12 23:26:28 +05:30
subroutine grid_mech_spectral_polarisation_init
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P
real(pReal), dimension(3,3) :: &
temp33_Real = 0.0_pReal
2020-06-16 22:17:19 +05:30
class (tNode), pointer :: &
2020-06-18 02:30:03 +05:30
num_grid, &
num_generic
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
PetscErrorCode :: ierr
PetscScalar, pointer, dimension(:,:,:,:) :: &
FandF_tau, & ! overall pointer to solution data
F, & ! specific (sub)pointer
F_tau ! specific (sub)pointer
2020-02-25 23:23:12 +05:30
PetscInt, dimension(0:worldsize-1) :: localK
2019-12-07 19:54:45 +05:30
integer(HID_T) :: fileHandle, groupHandle
integer :: fileUnit
2020-06-18 02:30:03 +05:30
character(len=pStringLen) :: &
fileName, &
petsc_options
2020-02-25 23:23:12 +05:30
2019-12-07 19:54:45 +05:30
write(6,'(/,a)') ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(6)
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
write(6,'(/,a)') ' Shanthraj et al., International Journal of Plasticity 66:3145, 2015'
write(6,'(a)') ' https://doi.org/10.1016/j.ijplas.2014.02.006'
2020-02-25 23:23:12 +05:30
2020-06-18 02:30:03 +05:30
!-------------------------------------------------------------------------------------------------
! read numerical parameters
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
petsc_options = num_generic%get_asString('petsc_options',defaultVal='')
2020-06-16 22:17:19 +05:30
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
2020-06-17 18:51:51 +05:30
num%update_gamma = num_grid%get_asBool('update_gamma',defaultVal=.false.)
2013-07-24 18:36:16 +05:30
2019-03-12 23:26:28 +05:30
!--------------------------------------------------------------------------------------------------
! set default and user defined options for PETSc
2019-04-15 19:23:46 +05:30
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr)
CHKERRQ(ierr)
call PETScOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_options),ierr)
CHKERRQ(ierr)
2019-03-12 23:26:28 +05:30
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! allocate global fields
2019-04-15 19:23:46 +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)
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)
2020-02-25 23:23:12 +05:30
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
2016-06-27 21:20:43 +05:30
! initialize solver specific parts of PETSc
2019-04-15 19:23:46 +05:30
call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr)
2020-02-25 23:23:12 +05:30
call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr)
2020-01-12 05:19:03 +05:30
localK = 0
localK(worldrank) = grid3
2019-04-15 19:23:46 +05:30
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
call DMDACreate3d(PETSC_COMM_WORLD, &
DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary
DMDA_STENCIL_BOX, & ! Moore (26) neighborhood around central point
grid(1),grid(2),grid(3), & ! global grid
1 , 1, worldsize, &
18, 0, & ! #dof (F tensor), ghost boundary width (domain overlap)
[grid(1)],[grid(2)],localK, & ! local grid
da,ierr) ! handle, error
CHKERRQ(ierr)
call SNESSetDM(snes,da,ierr); CHKERRQ(ierr) ! connect snes to da
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 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
2020-02-25 23:23:12 +05:30
CHKERRQ(ierr)
2019-04-15 19:23:46 +05:30
call SNESsetConvergenceTest(snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "converged"
CHKERRQ(ierr)
call SNESsetFromOptions(snes,ierr); CHKERRQ(ierr) ! pull it all together with additional CLI arguments
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
2020-02-25 23:23:12 +05:30
! init fields
2019-04-15 19:23:46 +05:30
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! places pointer on PETSc data
2019-10-28 17:47:05 +05:30
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
2020-02-25 23:23:12 +05:30
restartRead: if (interface_restartInc > 0) then
2019-12-07 19:54:45 +05:30
write(6,'(/,a,i0,a)') ' reading restart data of increment ', interface_restartInc, ' from file'
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName)
2019-10-28 17:47:05 +05:30
groupHandle = HDF5_openGroup(fileHandle,'solver')
2020-02-25 23:23:12 +05:30
2019-10-28 17:47:05 +05:30
call HDF5_read(groupHandle,F_aim, 'F_aim')
call HDF5_read(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_read(groupHandle,F_aimDot, 'F_aimDot')
call HDF5_read(groupHandle,F, 'F')
call HDF5_read(groupHandle,F_lastInc, 'F_lastInc')
call HDF5_read(groupHandle,F_tau, 'F_tau')
call HDF5_read(groupHandle,F_tau_lastInc,'F_tau_lastInc')
2020-02-25 23:23:12 +05:30
elseif (interface_restartInc == 0) then restartRead
2019-04-15 19:23:46 +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])
F_tau = 2.0_pReal*F
F_tau_lastInc = 2.0_pReal*F_lastInc
endif restartRead
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call Utilities_updateCoords(reshape(F,shape(F_lastInc)))
2019-04-15 19:23:46 +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
2019-12-03 00:53:50 +05:30
0.0_pReal) ! time increment
2019-04-15 19:23:46 +05:30
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer
2020-02-25 23:23:12 +05:30
restartRead2: if (interface_restartInc > 0) then
2019-12-07 19:54:45 +05:30
write(6,'(/,a,i0,a)') ' reading more restart data of increment ', interface_restartInc, ' from file'
2019-10-28 17:47:05 +05:30
call HDF5_read(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeGroup(groupHandle)
2019-04-15 19:23:46 +05:30
call HDF5_closeFile(fileHandle)
call MPI_File_open(PETSC_COMM_WORLD, trim(getSolverJobName())//'.C_ref', &
MPI_MODE_RDONLY,MPI_INFO_NULL,fileUnit,ierr)
call MPI_File_read(fileUnit,C_minMaxAvg,81,MPI_DOUBLE,MPI_STATUS_IGNORE,ierr)
call MPI_File_close(fileUnit,ierr)
endif restartRead2
2020-02-25 23:23:12 +05:30
2019-10-24 02:20:01 +05:30
call utilities_updateGamma(C_minMaxAvg)
call utilities_saveReferenceStiffness
2019-04-15 19:23:46 +05:30
C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg)
2020-02-25 23:23:12 +05:30
2019-03-12 23:26:28 +05:30
end subroutine grid_mech_spectral_polarisation_init
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief solution for the Polarisation scheme with internal iterations
!--------------------------------------------------------------------------------------------------
2019-03-12 23:26:28 +05:30
function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution)
2016-06-27 21:20:43 +05:30
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! input data for solution
2019-04-15 19:23:46 +05:30
character(len=*), intent(in) :: &
incInfoIn
2019-12-03 00:36:58 +05:30
real(pReal), intent(in) :: &
2019-04-15 19:23:46 +05:30
timeinc, & !< time increment of current solution
timeinc_old !< time increment of last successful increment
type(tBoundaryCondition), intent(in) :: &
stress_BC
2019-12-03 00:36:58 +05:30
type(rotation), intent(in) :: &
rotation_BC
2019-04-15 19:23:46 +05:30
type(tSolutionState) :: &
solution
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! PETSc Data
2020-02-25 23:23:12 +05:30
PetscErrorCode :: ierr
2019-04-15 19:23:46 +05:30
SNESConvergedReason :: reason
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
incInfo = incInfoIn
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! update stiffness (and gamma operator)
2019-04-15 19:23:46 +05:30
S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg)
if (num%update_gamma) then
2019-10-24 02:20:01 +05:30
call utilities_updateGamma(C_minMaxAvg)
2019-04-15 19:23:46 +05:30
C_scale = C_minMaxAvg
S_scale = math_invSym3333(C_minMaxAvg)
2020-02-25 23:23:12 +05:30
endif
!--------------------------------------------------------------------------------------------------
2020-02-25 23:23:12 +05:30
! set module wide available data
2019-04-15 19:23:46 +05:30
params%stress_mask = stress_BC%maskFloat
params%stress_BC = stress_BC%values
params%rotation_BC = rotation_BC
params%timeinc = timeinc
params%timeincOld = timeinc_old
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
2020-02-25 23:23:12 +05:30
! solve BVP
2019-04-15 19:23:46 +05:30
call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
! check convergence
2019-04-15 19:23:46 +05:30
call SNESGetConvergedReason(snes,reason,ierr); CHKERRQ(ierr)
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
solution%converged = reason > 0
solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll
terminallyIll = .false.
2013-07-24 18:36:16 +05:30
2019-03-12 23:26:28 +05:30
end function grid_mech_spectral_polarisation_solution
2013-07-24 18:36:16 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!> @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
!--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,&
deformation_BC,stress_BC,rotation_BC)
logical, intent(in) :: &
cutBack, &
guess
real(pReal), intent(in) :: &
timeinc_old, &
timeinc, &
2019-04-15 19:23:46 +05:30
loadCaseTime !< remaining time of current load case
type(tBoundaryCondition), intent(in) :: &
stress_BC, &
deformation_BC
2019-12-03 00:36:58 +05:30
type(rotation), intent(in) :: &
rotation_BC
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
integer :: i, j, k
real(pReal), dimension(3,3) :: F_lambda33
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
if (cutBack) then
C_volAvg = C_volAvgLastInc
C_minMaxAvg = C_minMaxAvgLastInc
2020-02-25 23:23:12 +05:30
else
C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess)
F_aim_lastInc = F_aim
2019-10-24 17:16:36 +05:30
!-----------------------------------------------------------------------------------------------
! calculate rate for aim
2019-04-15 19:23:46 +05:30
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
2019-04-15 19:23:46 +05:30
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values
2019-04-15 19:23:46 +05:30
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
2019-10-25 11:25:23 +05:30
Fdot = utilities_calculateRate(guess, &
F_lastInc,reshape(F,[3,3,grid(1),grid(2),grid3]),timeinc_old, &
2020-02-25 23:23:12 +05:30
rotation_BC%rotate(F_aimDot,active=.true.))
2019-10-25 11:25:23 +05:30
F_tauDot = utilities_calculateRate(guess, &
F_tau_lastInc,reshape(F_tau,[3,3,grid(1),grid(2),grid3]), timeinc_old, &
2020-02-25 23:23:12 +05:30
rotation_BC%rotate(F_aimDot,active=.true.))
2019-10-30 03:45:02 +05:30
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
2020-02-25 23:23:12 +05:30
2019-10-30 03:45:02 +05:30
materialpoint_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3])
endif
!--------------------------------------------------------------------------------------------------
! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * timeinc
2019-04-15 19:23:46 +05:30
F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average
2020-02-25 23:23:12 +05:30
rotation_BC%rotate(F_aim,active=.true.)),&
[9,grid(1),grid(2),grid3])
2019-04-15 19:23:46 +05:30
if (guess) then
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
else
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 = math_mul3333xx33(S_scale,matmul(F_lambda33, &
math_mul3333xx33(C_scale,&
matmul(transpose(F_lambda33),&
2019-12-03 00:36:58 +05:30
F_lambda33)-math_I3))*0.5_pReal) &
+ math_I3
2019-04-15 19:23:46 +05:30
F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
enddo; enddo; enddo
endif
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
2019-03-12 23:26:28 +05:30
end subroutine grid_mech_spectral_polarisation_forward
!--------------------------------------------------------------------------------------------------
!> @brief Age
!--------------------------------------------------------------------------------------------------
2019-12-07 19:54:45 +05:30
subroutine grid_mech_spectral_polarisation_updateCoords
PetscErrorCode :: ierr
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau
2020-02-25 23:23:12 +05:30
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
call utilities_updateCoords(FandF_tau(0:8,:,:,:))
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_updateCoords
2019-10-24 17:16:36 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Write current solver and constitutive data for restart to file
!--------------------------------------------------------------------------------------------------
2019-12-07 19:54:45 +05:30
subroutine grid_mech_spectral_polarisation_restartWrite
2019-10-24 17:16:36 +05:30
PetscErrorCode :: ierr
2019-12-07 19:54:45 +05:30
integer(HID_T) :: fileHandle, groupHandle
2019-10-24 17:16:36 +05:30
PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau
2019-12-07 19:54:45 +05:30
character(len=pStringLen) :: fileName
2019-10-24 17:16:36 +05:30
call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
F => FandF_tau(0: 8,:,:,:)
F_tau => FandF_tau(9:17,:,:,:)
2019-12-07 19:54:45 +05:30
write(6,'(a)') ' writing solver data required for restart to file'; flush(6)
2019-10-24 17:16:36 +05:30
2019-12-07 19:54:45 +05:30
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
fileHandle = HDF5_openFile(fileName,'w')
2019-10-28 17:47:05 +05:30
groupHandle = HDF5_addGroup(fileHandle,'solver')
2020-02-25 23:23:12 +05:30
2019-10-28 17:47:05 +05:30
call HDF5_write(groupHandle,F_aim, 'F_aim')
call HDF5_write(groupHandle,F_aim_lastInc,'F_aim_lastInc')
call HDF5_write(groupHandle,F_aimDot, 'F_aimDot')
call HDF5_write(groupHandle,F, 'F')
call HDF5_write(groupHandle,F_lastInc, 'F_lastInc')
call HDF5_write(groupHandle,F_tau, 'F_tau')
call HDF5_write(groupHandle,F_tau_lastInc,'F_tau_lastInc')
call HDF5_write(groupHandle,C_volAvg, 'C_volAvg')
call HDF5_write(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
call HDF5_closeGroup(groupHandle)
2019-10-24 17:16:36 +05:30
call HDF5_closeFile(fileHandle)
2020-02-25 23:23:12 +05:30
2019-10-25 04:12:59 +05:30
if(num%update_gamma) call utilities_saveReferenceStiffness
2020-02-25 23:23:12 +05:30
2019-10-24 17:16:36 +05:30
call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr)
end subroutine grid_mech_spectral_polarisation_restartWrite
!--------------------------------------------------------------------------------------------------
!> @brief convergence check
!--------------------------------------------------------------------------------------------------
2019-04-15 19:23:46 +05:30
subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dummy,ierr)
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
SNES :: snes_local
PetscInt, intent(in) :: PETScIter
PetscReal, intent(in) :: &
devNull1, &
devNull2, &
2020-02-25 23:23:12 +05:30
devNull3
2019-04-15 19:23:46 +05:30
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
2020-06-18 00:16:03 +05:30
integer :: &
itmin, & !< minimum number of iterations
itmax !< maximum number of iterations
2019-04-15 19:23:46 +05:30
real(pReal) :: &
curlTol, &
divTol, &
2020-06-17 00:19:58 +05:30
BCTol, &
2020-06-18 00:16:03 +05:30
eps_div_atol, & !< absolute tolerance for equilibrium
eps_div_rtol, & !< relative tolerance for equilibrium
eps_curl_atol, & !< absolute tolerance for compatibility
eps_curl_rtol, & !< relative tolerance for compatibility
eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC
eps_stress_rtol !< relative tolerance for fullfillment of stress BC
2020-06-17 00:19:58 +05:30
class(tNode), pointer :: &
2020-06-18 00:16:03 +05:30
num_grid, &
num_generic
!-----------------------------------------------------------------------------------
! reading numerical parameters and do sanity check
2020-06-17 00:19:58 +05:30
num_grid => numerics_root%get('grid',defaultVal=emptyDict)
eps_div_atol = num_grid%get_asFloat('eps_div_atol',defaultVal=1.0e-4_pReal)
eps_div_rtol = num_grid%get_asFloat('eps_div_rtol',defaultVal=5.0e-4_pReal)
eps_curl_atol = num_grid%get_asFloat('eps_curl_atol',defaultVal=1.0e-10_pReal)
eps_curl_rtol = num_grid%get_asFloat('eps_curl_rtol',defaultVal=5.0e-4_pReal)
eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=0.01_pReal)
2020-06-18 00:16:03 +05:30
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmin = num_generic%get_asInt('itmin',defaultVal=1)
itmax = num_generic%get_asInt('itmax',defaultVal=250)
if (eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
if (eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
if (eps_curl_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_curl_atol')
if (eps_curl_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_curl_rtol')
if (eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol')
if (eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol')
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
!------------------------------------------------------------------------------------
2020-02-25 23:23:12 +05:30
2020-06-17 00:19:58 +05:30
curlTol = max(maxval(abs(F_aim-math_I3))*eps_curl_rtol ,eps_curl_atol)
divTol = max(maxval(abs(P_av)) *eps_div_rtol ,eps_div_atol)
BCTol = max(maxval(abs(P_av)) *eps_stress_rtol,eps_stress_atol)
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
if ((totalIter >= itmin .and. &
all([ err_div /divTol, &
err_curl/curlTol, &
err_BC /BCTol ] < 1.0_pReal)) &
.or. terminallyIll) then
reason = 1
elseif (totalIter >= itmax) then
reason = -1
else
reason = 0
endif
!--------------------------------------------------------------------------------------------------
! report
2019-04-15 19:23:46 +05:30
write(6,'(1/,a)') ' ... reporting .............................................................'
write(6,'(1/,a,f12.2,a,es8.2,a,es9.2,a)') ' error divergence = ', &
err_div/divTol, ' (',err_div, ' / m, tol = ',divTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error curl = ', &
err_curl/curlTol,' (',err_curl,' -, tol = ',curlTol,')'
write(6, '(a,f12.2,a,es8.2,a,es9.2,a)') ' error BC = ', &
2020-02-25 23:23:12 +05:30
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
2019-04-15 19:23:46 +05:30
write(6,'(/,a)') ' ==========================================================================='
2020-02-25 23:23:12 +05:30
flush(6)
end subroutine converged
!--------------------------------------------------------------------------------------------------
2019-04-15 19:46:14 +05:30
!> @brief forms the residual vector
!--------------------------------------------------------------------------------------------------
subroutine formResidual(in, FandF_tau, &
residuum, dummy,ierr)
2019-04-15 19:23:46 +05:30
DMDALocalInfo, dimension(DMDA_LOCAL_INFO_SIZE) :: in !< DMDA info (needs to be named "in" for macros like XRANGE to work)
PetscScalar, dimension(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE), &
target, intent(in) :: FandF_tau
PetscScalar, dimension(3,3,2,X_RANGE,Y_RANGE,Z_RANGE),&
target, intent(out) :: residuum !< residuum field
PetscScalar, pointer, dimension(:,:,:,:,:) :: &
F, &
F_tau, &
residual_F, &
residual_F_tau
PetscInt :: &
PETScIter, &
nfuncs
PetscObject :: dummy
PetscErrorCode :: ierr
2020-06-17 00:19:58 +05:30
class(tNode), pointer :: &
2020-06-18 00:16:03 +05:30
num_grid, &
num_generic
2020-06-17 00:19:58 +05:30
real(pReal) :: &
2020-06-17 20:49:21 +05:30
polarAlpha, & !< polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme
polarBeta !< polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme
2019-04-15 19:23:46 +05:30
integer :: &
2020-06-18 00:16:03 +05:30
i, j, k, e, &
itmin, itmax
2019-04-15 19:23:46 +05:30
2020-06-18 00:16:03 +05:30
!--------------------------------------------------------------------------------------------------
! read numerical paramteters and do sanity checks
2020-06-17 00:19:58 +05:30
num_grid => numerics_root%get('grid',defaultVal = emptyDict)
polarAlpha = num_grid%get_asFloat('polaralpha',defaultVal=1.0_pReal)
polarBeta = num_grid%get_asFloat('polarbeta', defaultVal=1.0_pReal)
2020-06-18 00:16:03 +05:30
num_generic => numerics_root%get('generic',defaultVal=emptyDict)
itmin = num_generic%get_asInt('itmin',defaultVal=1)
itmax = num_generic%get_asInt('itmax',defaultVal=250)
if (itmax <= 1) call IO_error(301,ext_msg='itmax')
if (itmin > itmax .or. itmin < 1) call IO_error(301,ext_msg='itmin')
2020-06-17 18:51:51 +05:30
if (polarAlpha <= 0.0_pReal .or. polarAlpha > 2.0_pReal) call IO_error(301,ext_msg='polarAlpha')
if (polarBeta < 0.0_pReal .or. polarBeta > 2.0_pReal) call IO_error(301,ext_msg='polarBeta')
2020-06-18 00:16:03 +05:30
!---------------------------------------------------------------------------------------------------
2020-06-17 18:51:51 +05:30
2019-10-28 17:47:05 +05:30
F => FandF_tau(1:3,1:3,1,&
2019-04-15 19:23:46 +05:30
XG_RANGE,YG_RANGE,ZG_RANGE)
2019-10-28 17:47:05 +05:30
F_tau => FandF_tau(1:3,1:3,2,&
2019-04-15 19:23:46 +05:30
XG_RANGE,YG_RANGE,ZG_RANGE)
2019-10-28 17:47:05 +05:30
residual_F => residuum(1:3,1:3,1,&
2019-04-15 19:23:46 +05:30
X_RANGE, Y_RANGE, Z_RANGE)
2019-10-28 17:47:05 +05:30
residual_F_tau => residuum(1:3,1:3,2,&
2019-04-15 19:23:46 +05:30
X_RANGE, Y_RANGE, Z_RANGE)
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)
2020-02-25 23:23:12 +05:30
2019-04-15 19:23:46 +05:30
call SNESGetNumberFunctionEvals(snes,nfuncs,ierr); CHKERRQ(ierr)
call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr)
2019-04-15 19:23:46 +05:30
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
!--------------------------------------------------------------------------------------------------
! begin of new iteration
2019-04-15 19:23:46 +05:30
newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1
2019-12-11 23:54:29 +05:30
write(6,'(1x,a,3(a,i0))') trim(incInfo), ' @ Iteration ', itmin, '≤',totalIter, '≤', itmax
2019-04-15 19:23:46 +05:30
if (iand(debug_level(debug_spectral),debug_spectralRotation) /= 0) &
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
2020-02-25 23:23:12 +05:30
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
2019-04-15 19:23:46 +05:30
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
' deformation gradient aim =', transpose(F_aim)
flush(6)
endif newIteration
!--------------------------------------------------------------------------------------------------
2020-02-25 23:23:12 +05:30
!
2019-04-15 19:23:46 +05:30
tensorField_real = 0.0_pReal
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
tensorField_real(1:3,1:3,i,j,k) = &
polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
polarAlpha*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))
enddo; enddo; enddo
2020-02-25 23:23:12 +05:30
!--------------------------------------------------------------------------------------------------
2020-02-25 23:23:12 +05:30
! doing convolution in Fourier space
2019-04-15 19:23:46 +05:30
call utilities_FFTtensorForward
2020-02-25 23:23:12 +05:30
call utilities_fourierGammaConvolution(params%rotation_BC%rotate(polarBeta*F_aim,active=.true.))
2019-04-15 19:23:46 +05:30
call utilities_FFTtensorBackward
!--------------------------------------------------------------------------------------------------
2020-02-25 23:23:12 +05:30
! constructing residual
2019-04-15 19:23:46 +05:30
residual_F_tau = polarBeta*F - tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3)
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
2019-04-15 19:23:46 +05:30
call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, &
F - residual_F_tau/polarBeta,params%timeinc,params%rotation_BC)
2020-02-25 23:23:12 +05:30
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
!--------------------------------------------------------------------------------------------------
! stress BC handling
2019-04-15 19:23:46 +05:30
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 &
2020-02-25 23:23:12 +05:30
-params%rotation_BC%rotate(F_av)) + &
2019-04-15 19:23:46 +05:30
params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc
! calculate divergence
2019-04-15 19:23:46 +05:30
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
call utilities_FFTtensorForward
err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress
!--------------------------------------------------------------------------------------------------
! constructing residual
2019-04-15 19:23:46 +05:30
e = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
e = e + 1
residual_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - 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))) &
+ residual_F_tau(1:3,1:3,i,j,k)
enddo; enddo; enddo
2020-02-25 23:23:12 +05:30
!--------------------------------------------------------------------------------------------------
! calculating curl
2019-04-15 19:23:46 +05:30
tensorField_real = 0.0_pReal
tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F
call utilities_FFTtensorForward
err_curl = Utilities_curlRMS()
end subroutine formResidual
2019-03-12 23:26:28 +05:30
end module grid_mech_spectral_polarisation