685 lines
32 KiB
Fortran
685 lines
32 KiB
Fortran
!--------------------------------------------------------------------------------------------------
|
|
!> @author Arko Jyoti Bhattacharjee, Max-Planck-Institut für Eisenforschung GmbH
|
|
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
|
|
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
|
|
!> @brief Grid solver for mechanics: FEM
|
|
!--------------------------------------------------------------------------------------------------
|
|
module grid_mech_FEM
|
|
#include <petsc/finclude/petscsnes.h>
|
|
#include <petsc/finclude/petscdmda.h>
|
|
use PETScdmda
|
|
use PETScsnes
|
|
|
|
use prec
|
|
use parallelization
|
|
use DAMASK_interface
|
|
use IO
|
|
use HDF5_utilities
|
|
use math
|
|
use spectral_utilities
|
|
use FEsolving
|
|
use config
|
|
use homogenization
|
|
use discretization
|
|
use discretization_grid
|
|
|
|
implicit none
|
|
private
|
|
|
|
type(tSolutionParams) :: params
|
|
|
|
type :: tNumerics
|
|
integer :: &
|
|
itmin, & !< minimum number of iterations
|
|
itmax !< maximum number of iterations
|
|
real(pReal) :: &
|
|
eps_div_atol, & !< absolute tolerance for equilibrium
|
|
eps_div_rtol, & !< relative tolerance for equilibrium
|
|
eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC
|
|
eps_stress_rtol !< relative tolerance for fullfillment of stress BC
|
|
end type tNumerics
|
|
|
|
type(tNumerics) :: num ! numerics parameters. Better name?
|
|
|
|
logical :: debugRotation
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! PETSc data
|
|
DM :: mech_grid
|
|
SNES :: mech_snes
|
|
Vec :: solution_current, solution_lastInc, solution_rate
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! common pointwise data
|
|
real(pReal), dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc
|
|
real(pReal) :: detJ
|
|
real(pReal), dimension(3) :: delta
|
|
real(pReal), dimension(3,8) :: BMat
|
|
real(pReal), dimension(8,8) :: HGMat
|
|
PetscInt :: xstart,ystart,zstart,xend,yend,zend
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! stress, stiffness and compliance average etc.
|
|
real(pReal), 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
|
|
P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress
|
|
P_aim = 0.0_pReal
|
|
character(len=:), allocatable :: incInfo !< time and increment information
|
|
real(pReal), dimension(3,3,3,3) :: &
|
|
C_volAvg = 0.0_pReal, & !< current volume average stiffness
|
|
C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness
|
|
S = 0.0_pReal !< current compliance (filled up with zeros)
|
|
|
|
real(pReal) :: &
|
|
err_BC !< deviation from stress BC
|
|
|
|
integer :: &
|
|
totalIter = 0 !< total iteration in current increment
|
|
|
|
public :: &
|
|
grid_mech_FEM_init, &
|
|
grid_mech_FEM_solution, &
|
|
grid_mech_FEM_forward, &
|
|
grid_mech_FEM_updateCoords, &
|
|
grid_mech_FEM_restartWrite
|
|
|
|
contains
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief allocates all necessary fields and fills them with data, potentially from restart info
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine grid_mech_FEM_init
|
|
|
|
real(pReal) :: HGCoeff = 0.0e-2_pReal
|
|
real(pReal), dimension(3,3) :: &
|
|
temp33_Real = 0.0_pReal
|
|
real(pReal), dimension(4,8) :: &
|
|
HGcomp = reshape([ 1.0_pReal, 1.0_pReal, 1.0_pReal,-1.0_pReal, &
|
|
1.0_pReal,-1.0_pReal,-1.0_pReal, 1.0_pReal, &
|
|
-1.0_pReal, 1.0_pReal,-1.0_pReal, 1.0_pReal, &
|
|
-1.0_pReal,-1.0_pReal, 1.0_pReal,-1.0_pReal, &
|
|
-1.0_pReal,-1.0_pReal, 1.0_pReal, 1.0_pReal, &
|
|
-1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, &
|
|
1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, &
|
|
1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8])
|
|
real(pReal), dimension(3,3,3,3) :: devNull
|
|
PetscErrorCode :: ierr
|
|
PetscScalar, pointer, dimension(:,:,:,:) :: &
|
|
u_current,u_lastInc
|
|
PetscInt, dimension(0:worldsize-1) :: localK
|
|
integer(HID_T) :: fileHandle, groupHandle
|
|
character(len=pStringLen) :: &
|
|
fileName
|
|
class(tNode), pointer :: &
|
|
num_grid, &
|
|
debug_grid
|
|
|
|
print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(IO_STDOUT)
|
|
|
|
!-------------------------------------------------------------------------------------------------
|
|
! debugging options
|
|
debug_grid => config_debug%get('grid', defaultVal=emptyList)
|
|
debugRotation = debug_grid%contains('rotation')
|
|
|
|
!-------------------------------------------------------------------------------------------------
|
|
! read numerical parameters and do sanity checks
|
|
num_grid => config_numerics%get('grid',defaultVal=emptyDict)
|
|
num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal)
|
|
num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal)
|
|
num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal)
|
|
num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal)
|
|
num%itmin = num_grid%get_asInt ('itmin', defaultVal=1)
|
|
num%itmax = num_grid%get_asInt ('itmax', defaultVal=250)
|
|
|
|
if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol')
|
|
if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol')
|
|
if (num%eps_stress_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_stress_atol')
|
|
if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol')
|
|
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
|
|
if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin')
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! set default and user defined options for PETSc
|
|
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type newtonls -mech_ksp_type fgmres &
|
|
&-mech_ksp_max_it 25 -mech_pc_type ml -mech_mg_levels_ksp_type chebyshev',ierr)
|
|
CHKERRQ(ierr)
|
|
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_grid%get_asString('petsc_options',defaultVal=''),ierr)
|
|
CHKERRQ(ierr)
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! allocate global fields
|
|
allocate(F (3,3,grid(1),grid(2),grid3),source = 0.0_pReal)
|
|
allocate(P_current (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)
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! initialize solver specific parts of PETSc
|
|
call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr); CHKERRQ(ierr)
|
|
call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr)
|
|
localK = 0
|
|
localK(worldrank) = grid3
|
|
call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr)
|
|
call DMDACreate3d(PETSC_COMM_WORLD, &
|
|
DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, &
|
|
DMDA_STENCIL_BOX, &
|
|
grid(1),grid(2),grid(3), &
|
|
1, 1, worldsize, &
|
|
3, 1, &
|
|
[grid(1)],[grid(2)],localK, &
|
|
mech_grid,ierr)
|
|
CHKERRQ(ierr)
|
|
call SNESSetDM(mech_snes,mech_grid,ierr); CHKERRQ(ierr)
|
|
call DMsetFromOptions(mech_grid,ierr); CHKERRQ(ierr)
|
|
call DMsetUp(mech_grid,ierr); CHKERRQ(ierr)
|
|
call DMDASetUniformCoordinates(mech_grid,0.0_pReal,geomSize(1),0.0_pReal,geomSize(2),0.0_pReal,geomSize(3),ierr)
|
|
CHKERRQ(ierr)
|
|
call DMCreateGlobalVector(mech_grid,solution_current,ierr); CHKERRQ(ierr)
|
|
call DMCreateGlobalVector(mech_grid,solution_lastInc,ierr); CHKERRQ(ierr)
|
|
call DMCreateGlobalVector(mech_grid,solution_rate ,ierr); CHKERRQ(ierr)
|
|
call DMSNESSetFunctionLocal(mech_grid,formResidual,PETSC_NULL_SNES,ierr)
|
|
CHKERRQ(ierr)
|
|
call DMSNESSetJacobianLocal(mech_grid,formJacobian,PETSC_NULL_SNES,ierr)
|
|
CHKERRQ(ierr)
|
|
call SNESSetConvergenceTest(mech_snes,converged,PETSC_NULL_SNES,PETSC_NULL_FUNCTION,ierr) ! specify custom convergence check function "_converged"
|
|
CHKERRQ(ierr)
|
|
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) ! ignore linear solve failures
|
|
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) ! pull it all together with additional cli arguments
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! init fields
|
|
call VecSet(solution_current,0.0_pReal,ierr);CHKERRQ(ierr)
|
|
call VecSet(solution_lastInc,0.0_pReal,ierr);CHKERRQ(ierr)
|
|
call VecSet(solution_rate ,0.0_pReal,ierr);CHKERRQ(ierr)
|
|
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
|
|
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
|
|
|
|
call DMDAGetCorners(mech_grid,xstart,ystart,zstart,xend,yend,zend,ierr) ! local grid extent
|
|
CHKERRQ(ierr)
|
|
xend = xstart+xend-1
|
|
yend = ystart+yend-1
|
|
zend = zstart+zend-1
|
|
delta = geomSize/real(grid,pReal) ! grid spacing
|
|
detJ = product(delta) ! cell volume
|
|
|
|
BMat = reshape(real([-1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), &
|
|
1.0_pReal/delta(1),-1.0_pReal/delta(2),-1.0_pReal/delta(3), &
|
|
-1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), &
|
|
1.0_pReal/delta(1), 1.0_pReal/delta(2),-1.0_pReal/delta(3), &
|
|
-1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), &
|
|
1.0_pReal/delta(1),-1.0_pReal/delta(2), 1.0_pReal/delta(3), &
|
|
-1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3), &
|
|
1.0_pReal/delta(1), 1.0_pReal/delta(2), 1.0_pReal/delta(3)],pReal), [3,8])/4.0_pReal ! shape function derivative matrix
|
|
|
|
HGMat = matmul(transpose(HGcomp),HGcomp) &
|
|
* HGCoeff*(delta(1)*delta(2) + delta(2)*delta(3) + delta(3)*delta(1))/16.0_pReal ! hourglass stabilization matrix
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! init fields
|
|
restartRead: if (interface_restartInc > 0) then
|
|
print'(/,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)
|
|
groupHandle = HDF5_openGroup(fileHandle,'solver')
|
|
|
|
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,u_current, 'u')
|
|
call HDF5_read(groupHandle,u_lastInc, 'u_lastInc')
|
|
|
|
elseif (interface_restartInc == 0) then restartRead
|
|
F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity
|
|
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
|
|
endif restartRead
|
|
|
|
homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
|
|
call utilities_updateCoords(F)
|
|
call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
|
|
F, & ! target F
|
|
0.0_pReal) ! time increment
|
|
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr)
|
|
CHKERRQ(ierr)
|
|
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr)
|
|
CHKERRQ(ierr)
|
|
|
|
restartRead2: if (interface_restartInc > 0) then
|
|
print'(a,i0,a)', ' reading more restart data of increment ', interface_restartInc, ' from file'
|
|
call HDF5_read(groupHandle,C_volAvg, 'C_volAvg')
|
|
call HDF5_read(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
|
|
|
|
call HDF5_closeGroup(groupHandle)
|
|
call HDF5_closeFile(fileHandle)
|
|
|
|
endif restartRead2
|
|
|
|
end subroutine grid_mech_FEM_init
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief solution for the FEM scheme with internal iterations
|
|
!--------------------------------------------------------------------------------------------------
|
|
function grid_mech_FEM_solution(incInfoIn) result(solution)
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! input data for solution
|
|
character(len=*), intent(in) :: &
|
|
incInfoIn
|
|
type(tSolutionState) :: &
|
|
solution
|
|
!--------------------------------------------------------------------------------------------------
|
|
! PETSc Data
|
|
PetscErrorCode :: ierr
|
|
SNESConvergedReason :: reason
|
|
|
|
incInfo = incInfoIn
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! update stiffness (and gamma operator)
|
|
S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg)
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! solve BVP
|
|
call SNESsolve(mech_snes,PETSC_NULL_VEC,solution_current,ierr); CHKERRQ(ierr)
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! check convergence
|
|
call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr)
|
|
|
|
solution%converged = reason > 0
|
|
solution%iterationsNeeded = totalIter
|
|
solution%termIll = terminallyIll
|
|
terminallyIll = .false.
|
|
|
|
end function grid_mech_FEM_solution
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @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_FEM_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, &
|
|
loadCaseTime !< remaining time of current load case
|
|
type(tBoundaryCondition), intent(in) :: &
|
|
stress_BC, &
|
|
deformation_BC
|
|
type(rotation), intent(in) :: &
|
|
rotation_BC
|
|
PetscErrorCode :: ierr
|
|
PetscScalar, pointer, dimension(:,:,:,:) :: &
|
|
u_current,u_lastInc
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! set module wide available data
|
|
params%stress_mask = stress_BC%mask
|
|
params%rotation_BC = rotation_BC
|
|
params%timeinc = timeinc
|
|
params%timeincOld = timeinc_old
|
|
|
|
|
|
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
|
|
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
|
|
|
|
if (cutBack) then
|
|
C_volAvg = C_volAvgLastInc
|
|
else
|
|
C_volAvgLastInc = C_volAvg
|
|
|
|
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess)
|
|
F_aim_lastInc = F_aim
|
|
|
|
!-----------------------------------------------------------------------------------------------
|
|
! calculate rate for aim
|
|
if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F
|
|
F_aimDot = F_aimDot &
|
|
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask)
|
|
elseif(deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
|
|
F_aimDot = F_aimDot &
|
|
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask)
|
|
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
|
|
F_aimDot = F_aimDot &
|
|
+ merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask)
|
|
endif
|
|
|
|
if (guess) then
|
|
call VecWAXPY(solution_rate,-1.0_pReal,solution_lastInc,solution_current,ierr)
|
|
CHKERRQ(ierr)
|
|
call VecScale(solution_rate,1.0_pReal/timeinc_old,ierr); CHKERRQ(ierr)
|
|
else
|
|
call VecSet(solution_rate,0.0_pReal,ierr); CHKERRQ(ierr)
|
|
endif
|
|
call VecCopy(solution_current,solution_lastInc,ierr); CHKERRQ(ierr)
|
|
|
|
F_lastInc = F
|
|
|
|
homogenization_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
|
|
if (stress_BC%myType=='P') then
|
|
P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask)
|
|
elseif (stress_BC%myType=='dot_P') then !UNTESTED
|
|
P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask)
|
|
endif
|
|
|
|
call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr)
|
|
|
|
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr)
|
|
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr)
|
|
|
|
end subroutine grid_mech_FEM_forward
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief Age
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine grid_mech_FEM_updateCoords
|
|
|
|
call utilities_updateCoords(F)
|
|
|
|
end subroutine grid_mech_FEM_updateCoords
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief Write current solver and constitutive data for restart to file
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine grid_mech_FEM_restartWrite
|
|
|
|
PetscErrorCode :: ierr
|
|
integer(HID_T) :: fileHandle, groupHandle
|
|
PetscScalar, dimension(:,:,:,:), pointer :: u_current,u_lastInc
|
|
character(len=pStringLen) :: fileName
|
|
|
|
call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr)
|
|
call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr)
|
|
|
|
print*, 'writing solver data required for restart to file'; flush(IO_STDOUT)
|
|
|
|
write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5'
|
|
fileHandle = HDF5_openFile(fileName,'w')
|
|
groupHandle = HDF5_addGroup(fileHandle,'solver')
|
|
|
|
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,u_current, 'u')
|
|
call HDF5_write(groupHandle,u_lastInc, 'u_lastInc')
|
|
|
|
call HDF5_write(groupHandle,C_volAvg, 'C_volAvg')
|
|
call HDF5_write(groupHandle,C_volAvgLastInc,'C_volAvgLastInc')
|
|
|
|
call HDF5_closeGroup(groupHandle)
|
|
call HDF5_closeFile(fileHandle)
|
|
|
|
call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr)
|
|
call DMDAVecRestoreArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr);CHKERRQ(ierr)
|
|
|
|
end subroutine grid_mech_FEM_restartWrite
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief convergence check
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,ierr)
|
|
|
|
SNES :: snes_local
|
|
PetscInt, intent(in) :: PETScIter
|
|
PetscReal, intent(in) :: &
|
|
devNull1, &
|
|
devNull2, &
|
|
fnorm
|
|
SNESConvergedReason :: reason
|
|
PetscObject :: dummy
|
|
PetscErrorCode :: ierr
|
|
real(pReal) :: &
|
|
err_div, &
|
|
divTol, &
|
|
BCTol
|
|
|
|
err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ
|
|
divTol = max(maxval(abs(P_av))*num%eps_div_rtol ,num%eps_div_atol)
|
|
BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
|
|
|
|
if ((totalIter >= num%itmin .and. &
|
|
all([ err_div/divTol, &
|
|
err_BC /BCTol ] < 1.0_pReal)) &
|
|
.or. terminallyIll) then
|
|
reason = 1
|
|
elseif (totalIter >= num%itmax) then
|
|
reason = -1
|
|
else
|
|
reason = 0
|
|
endif
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! report
|
|
print'(1/,a)', ' ... reporting .............................................................'
|
|
print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
|
|
err_div/divTol, ' (',err_div,' / m, tol = ',divTol,')'
|
|
print'(a,f12.2,a,es8.2,a,es9.2,a)', ' error stress BC = ', &
|
|
err_BC/BCTol, ' (',err_BC, ' Pa, tol = ',BCTol,')'
|
|
print'(/,a)', ' ==========================================================================='
|
|
flush(IO_STDOUT)
|
|
|
|
end subroutine converged
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief forms the residual vector
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine formResidual(da_local,x_local, &
|
|
f_local,dummy,ierr)
|
|
|
|
DM :: da_local
|
|
Vec :: x_local, f_local
|
|
PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal
|
|
PetscScalar, dimension(8,3) :: x_elem, f_elem
|
|
PetscInt :: i, ii, j, jj, k, kk, ctr, ele
|
|
PetscInt :: &
|
|
PETScIter, &
|
|
nfuncs
|
|
PetscObject :: dummy
|
|
PetscErrorCode :: ierr
|
|
real(pReal), dimension(3,3,3,3) :: devNull
|
|
|
|
call SNESGetNumberFunctionEvals(mech_snes,nfuncs,ierr); CHKERRQ(ierr)
|
|
call SNESGetIterationNumber(mech_snes,PETScIter,ierr); CHKERRQ(ierr)
|
|
|
|
if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! begin of new iteration
|
|
newIteration: if (totalIter <= PETScIter) then
|
|
totalIter = totalIter + 1
|
|
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax
|
|
if (debugRotation) &
|
|
write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
|
' deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
|
|
write(IO_STDOUT,'(/,a,/,3(3(f12.7,1x)/))',advance='no') &
|
|
' deformation gradient aim =', transpose(F_aim)
|
|
flush(IO_STDOUT)
|
|
endif newIteration
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! get deformation gradient
|
|
call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
|
|
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
|
|
ctr = 0
|
|
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
|
|
ctr = ctr + 1
|
|
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
|
|
enddo; enddo; enddo
|
|
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
|
|
F(1:3,1:3,ii,jj,kk) = params%rotation_BC%rotate(F_aim,active=.true.) + transpose(matmul(BMat,x_elem))
|
|
enddo; enddo; enddo
|
|
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! evaluate constitutive response
|
|
call utilities_constitutiveResponse(P_current,&
|
|
P_av,C_volAvg,devNull, &
|
|
F,params%timeinc,params%rotation_BC)
|
|
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! stress BC handling
|
|
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
|
|
err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask)))
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! constructing residual
|
|
call VecSet(f_local,0.0_pReal,ierr);CHKERRQ(ierr)
|
|
call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr)
|
|
call DMDAVecGetArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
|
|
ele = 0
|
|
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
|
|
ctr = 0
|
|
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
|
|
ctr = ctr + 1
|
|
x_elem(ctr,1:3) = x_scal(0:2,i+ii,j+jj,k+kk)
|
|
enddo; enddo; enddo
|
|
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
|
|
ele = ele + 1
|
|
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + &
|
|
matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,1,ele) + &
|
|
homogenization_dPdF(2,2,2,2,1,ele) + &
|
|
homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
|
|
ctr = 0
|
|
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
|
|
ctr = ctr + 1
|
|
f_scal(0:2,i+ii,j+jj,k+kk) = f_scal(0:2,i+ii,j+jj,k+kk) + f_elem(ctr,1:3)
|
|
enddo; enddo; enddo
|
|
enddo; enddo; enddo
|
|
call DMDAVecRestoreArrayF90(da_local,x_local,x_scal,ierr);CHKERRQ(ierr)
|
|
call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr)
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! applying boundary conditions
|
|
call DMDAVecGetArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr)
|
|
if (zstart == 0) then
|
|
f_scal(0:2,xstart,ystart,zstart) = 0.0
|
|
f_scal(0:2,xend+1,ystart,zstart) = 0.0
|
|
f_scal(0:2,xstart,yend+1,zstart) = 0.0
|
|
f_scal(0:2,xend+1,yend+1,zstart) = 0.0
|
|
endif
|
|
if (zend + 1 == grid(3)) then
|
|
f_scal(0:2,xstart,ystart,zend+1) = 0.0
|
|
f_scal(0:2,xend+1,ystart,zend+1) = 0.0
|
|
f_scal(0:2,xstart,yend+1,zend+1) = 0.0
|
|
f_scal(0:2,xend+1,yend+1,zend+1) = 0.0
|
|
endif
|
|
call DMDAVecRestoreArrayF90(da_local,f_local,f_scal,ierr);CHKERRQ(ierr)
|
|
|
|
end subroutine formResidual
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief forms the FEM stiffness matrix
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
|
|
|
|
DM :: da_local
|
|
Vec :: x_local, coordinates
|
|
Mat :: Jac_pre, Jac
|
|
MatStencil,dimension(4,24) :: row, col
|
|
PetscScalar,pointer,dimension(:,:,:,:) :: x_scal
|
|
PetscScalar,dimension(24,24) :: K_ele
|
|
PetscScalar,dimension(9,24) :: BMatFull
|
|
PetscInt :: i, ii, j, jj, k, kk, ctr, ele
|
|
PetscInt,dimension(3),parameter :: rows = [0, 1, 2]
|
|
PetscScalar :: diag
|
|
PetscObject :: dummy
|
|
MatNullSpace :: matnull
|
|
PetscErrorCode :: ierr
|
|
|
|
BMatFull = 0.0
|
|
BMatFull(1:3,1 :8 ) = BMat
|
|
BMatFull(4:6,9 :16) = BMat
|
|
BMatFull(7:9,17:24) = BMat
|
|
call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr)
|
|
call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr)
|
|
call MatZeroEntries(Jac,ierr); CHKERRQ(ierr)
|
|
ele = 0
|
|
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
|
|
ctr = 0
|
|
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
|
|
ctr = ctr + 1
|
|
col(MatStencil_i,ctr ) = i+ii
|
|
col(MatStencil_j,ctr ) = j+jj
|
|
col(MatStencil_k,ctr ) = k+kk
|
|
col(MatStencil_c,ctr ) = 0
|
|
col(MatStencil_i,ctr+8 ) = i+ii
|
|
col(MatStencil_j,ctr+8 ) = j+jj
|
|
col(MatStencil_k,ctr+8 ) = k+kk
|
|
col(MatStencil_c,ctr+8 ) = 1
|
|
col(MatStencil_i,ctr+16) = i+ii
|
|
col(MatStencil_j,ctr+16) = j+jj
|
|
col(MatStencil_k,ctr+16) = k+kk
|
|
col(MatStencil_c,ctr+16) = 2
|
|
enddo; enddo; enddo
|
|
row = col
|
|
ele = ele + 1
|
|
K_ele = 0.0
|
|
K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + &
|
|
homogenization_dPdF(2,2,2,2,1,ele) + &
|
|
homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
|
|
K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + &
|
|
homogenization_dPdF(2,2,2,2,1,ele) + &
|
|
homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
|
|
K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + &
|
|
homogenization_dPdF(2,2,2,2,1,ele) + &
|
|
homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
|
|
K_ele = K_ele + &
|
|
matmul(transpose(BMatFull), &
|
|
matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,1,ele), &
|
|
shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ
|
|
call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr)
|
|
CHKERRQ(ierr)
|
|
enddo; enddo; enddo
|
|
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
|
|
call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
|
|
call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
|
|
call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
! applying boundary conditions
|
|
diag = (C_volAvg(1,1,1,1)/delta(1)**2.0_pReal + &
|
|
C_volAvg(2,2,2,2)/delta(2)**2.0_pReal + &
|
|
C_volAvg(3,3,3,3)/delta(3)**2.0_pReal)*detJ
|
|
call MatZeroRowsColumns(Jac,size(rows),rows,diag,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr)
|
|
CHKERRQ(ierr)
|
|
call DMGetGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr)
|
|
call DMDAVecGetArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr)
|
|
ele = 0
|
|
do k = zstart, zend; do j = ystart, yend; do i = xstart, xend
|
|
ele = ele + 1
|
|
x_scal(0:2,i,j,k) = discretization_IPcoords(1:3,ele)
|
|
enddo; enddo; enddo
|
|
call DMDAVecRestoreArrayF90(da_local,coordinates,x_scal,ierr);CHKERRQ(ierr) ! initialize to undeformed coordinates (ToDo: use ip coordinates)
|
|
call MatNullSpaceCreateRigidBody(coordinates,matnull,ierr);CHKERRQ(ierr) ! get rigid body deformation modes
|
|
call DMRestoreGlobalVector(da_local,coordinates,ierr);CHKERRQ(ierr)
|
|
call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
|
|
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
|
|
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)
|
|
|
|
end subroutine formJacobian
|
|
|
|
end module grid_mech_FEM
|