DAMASK_EICMD/src/mesh/mesh_mech_FEM.f90

706 lines
32 KiB
Fortran
Raw Normal View History

2018-08-17 03:44:25 +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 FEM PETSc solver
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module mesh_mechanical_FEM
#include <petsc/finclude/petscdmplex.h>
#include <petsc/finclude/petscdm.h>
2018-08-18 17:35:57 +05:30
#include <petsc/finclude/petsc.h>
2019-06-11 13:18:07 +05:30
use PETScsnes
use PETScDM
use PETScDMplex
use PETScDT
2019-06-11 13:18:07 +05:30
use prec
use FEM_utilities
use discretization
2020-03-20 19:38:07 +05:30
use discretization_mesh
2019-06-11 13:18:07 +05:30
use DAMASK_interface
2020-08-15 19:32:10 +05:30
use config
use IO
2020-03-20 19:25:10 +05:30
use FEM_quadrature
2019-06-11 13:18:07 +05:30
use homogenization
use math
2019-06-11 13:18:07 +05:30
implicit none
private
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! derived types
type tSolutionParams
2019-06-11 13:18:07 +05:30
type(tFieldBC) :: fieldBC
real(pReal) :: timeinc
end type tSolutionParams
2019-06-11 13:18:07 +05:30
type(tSolutionParams) :: params
2018-08-17 03:44:25 +05:30
2020-06-24 20:15:13 +05:30
type, private :: tNumerics
integer :: &
integrationOrder, & !< order of quadrature rule required
itmax
logical :: &
BBarStabilisation
real(pReal) :: &
eps_struct_atol, & !< absolute tolerance for mechanical equilibrium
eps_struct_rtol !< relative tolerance for mechanical equilibrium
end type tNumerics
type(tNumerics), private :: num
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! PETSc data
2021-02-09 03:51:53 +05:30
SNES :: mechanical_snes
2019-06-11 13:18:07 +05:30
Vec :: solution, solution_rate, solution_local
PetscInt :: dimPlex, cellDof, nQuadrature, nBasis
PetscReal, allocatable, target :: qPoints(:), qWeights(:)
MatNullSpace :: matnull
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! stress, stiffness and compliance average etc.
2020-01-26 16:28:13 +05:30
character(len=pStringLen) :: incInfo
2019-06-11 13:18:07 +05:30
real(pReal), dimension(3,3) :: &
P_av = 0.0_pReal
logical :: ForwardData
real(pReal), parameter :: eps = 1.0e-18_pReal
2019-06-11 13:18:07 +05:30
public :: &
2021-02-09 03:51:53 +05:30
FEM_mechanical_init, &
FEM_mechanical_solution, &
FEM_mechanical_forward, &
FEM_mechanical_updateCoords
2019-06-11 13:18:07 +05:30
2018-08-17 03:44:25 +05:30
contains
!--------------------------------------------------------------------------------------------------
2019-10-24 02:08:17 +05:30
!> @brief allocates all neccessary fields and fills them with data
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
subroutine FEM_mechanical_init(fieldBC)
2019-06-11 13:18:07 +05:30
type(tFieldBC), intent(in) :: fieldBC
2021-02-09 03:51:53 +05:30
DM :: mechanical_mesh
2019-06-11 13:18:07 +05:30
PetscFE :: mechFE
PetscQuadrature :: mechQuad, functional
PetscDS :: mechDS
PetscDualSpace :: mechDualSpace
2020-05-13 16:28:39 +05:30
DMLabel, dimension(:),pointer :: nolabel=> NULL()
2019-06-11 13:18:07 +05:30
DMLabel :: BCLabel
2019-06-11 13:18:07 +05:30
PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint
PetscInt :: numBC, bcSize, nc, &
field, faceSet, topologDim, nNodalPoints, &
cellStart, cellEnd, cell, basis
2019-06-11 13:18:07 +05:30
IS :: bcPoint
IS, dimension(:), pointer :: pBcComps, pBcPoints
2019-06-11 13:18:07 +05:30
PetscSection :: section
2019-06-11 13:18:07 +05:30
PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, &
nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ
2019-06-11 13:18:07 +05:30
PetscReal :: detJ
PetscReal, allocatable, target :: cellJMat(:,:)
PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:)
character(len=*), parameter :: prefix = 'mechFE_'
2019-06-11 13:18:07 +05:30
PetscErrorCode :: ierr
2020-06-16 22:45:01 +05:30
class(tNode), pointer :: &
2020-06-19 06:02:33 +05:30
num_mesh
2020-06-18 00:16:03 +05:30
2020-09-22 16:39:12 +05:30
print'(/,a)', ' <<<+- FEM_mech init -+>>>'; flush(IO_STDOUT)
2018-08-17 03:44:25 +05:30
2020-06-18 00:16:03 +05:30
!-----------------------------------------------------------------------------
! read numerical parametes and do sanity checks
2020-09-13 14:47:49 +05:30
num_mesh => config_numerics%get('mesh',defaultVal=emptyDict)
2020-06-24 20:15:13 +05:30
num%integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2)
num%itmax = num_mesh%get_asInt('itmax',defaultVal=250)
num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
num%eps_struct_atol = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal)
num%eps_struct_rtol = num_mesh%get_asFloat('eps_struct_rtol', defaultVal = 1.0e-4_pReal)
2020-06-18 00:16:03 +05:30
2020-06-24 20:15:13 +05:30
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
if (num%eps_struct_rtol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_rtol')
if (num%eps_struct_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_struct_atol')
2020-06-16 22:45:01 +05:30
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! Setup FEM mech mesh
2021-02-09 03:51:53 +05:30
call DMClone(geomMesh,mechanical_mesh,ierr); CHKERRQ(ierr)
call DMGetDimension(mechanical_mesh,dimPlex,ierr); CHKERRQ(ierr)
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! Setup FEM mech discretization
2020-06-24 20:15:13 +05:30
qPoints = FEM_quadrature_points( dimPlex,num%integrationOrder)%p
qWeights = FEM_quadrature_weights(dimPlex,num%integrationOrder)%p
nQuadrature = FEM_nQuadrature( dimPlex,num%integrationOrder)
2020-03-20 19:25:10 +05:30
qPointsP => qPoints
2019-06-11 13:18:07 +05:30
qWeightsP => qWeights
call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr)
CHKERRQ(ierr)
nc = dimPlex
call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr)
CHKERRQ(ierr)
call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, &
2020-06-24 20:15:13 +05:30
num%integrationOrder,mechFE,ierr); CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
nBasis = nBasis/nc
2021-02-09 03:51:53 +05:30
call DMAddField(mechanical_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr)
call DMCreateDS(mechanical_mesh,ierr); CHKERRQ(ierr)
call DMGetDS(mechanical_mesh,mechDS,ierr); CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr)
call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr)
call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr)
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! Setup FEM mech boundary conditions
2021-02-09 03:51:53 +05:30
call DMGetLabel(mechanical_mesh,'Face Sets',BCLabel,ierr); CHKERRQ(ierr)
call DMPlexLabelComplete(mechanical_mesh,BCLabel,ierr); CHKERRQ(ierr)
call DMGetLocalSection(mechanical_mesh,section,ierr); CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
allocate(pnumComp(1), source=dimPlex)
2020-01-21 19:28:35 +05:30
allocate(pnumDof(0:dimPlex), source = 0)
2019-06-11 13:18:07 +05:30
do topologDim = 0, dimPlex
2021-02-09 03:51:53 +05:30
call DMPlexGetDepthStratum(mechanical_mesh,topologDim,cellStart,cellEnd,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
2020-01-21 19:28:35 +05:30
call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
enddo
numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1
enddo; enddo
allocate(pbcField(numBC), source=0)
allocate(pbcComps(numBC))
allocate(pbcPoints(numBC))
numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) then
numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),ierr)
CHKERRQ(ierr)
2021-02-09 03:51:53 +05:30
call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
if (bcSize > 0) then
2021-02-09 03:51:53 +05:30
call DMGetStratumIS(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
call ISGetIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,pbcPoints(numBC),ierr)
CHKERRQ(ierr)
call ISRestoreIndicesF90(bcPoint,pBcPoint,ierr); CHKERRQ(ierr)
call ISDestroy(bcPoint,ierr); CHKERRQ(ierr)
else
call ISCreateGeneral(PETSC_COMM_WORLD,0,[0],PETSC_COPY_VALUES,pbcPoints(numBC),ierr)
CHKERRQ(ierr)
endif
endif
2019-06-11 13:18:07 +05:30
enddo; enddo
2021-02-09 03:51:53 +05:30
call DMPlexCreateSection(mechanical_mesh,nolabel,pNumComp,pNumDof, &
numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
2021-02-09 03:51:53 +05:30
call DMSetSection(mechanical_mesh,section,ierr); CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
do faceSet = 1, numBC
call ISDestroy(pbcPoints(faceSet),ierr); CHKERRQ(ierr)
enddo
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc
2021-02-09 03:51:53 +05:30
call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,ierr);CHKERRQ(ierr)
call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',ierr);CHKERRQ(ierr)
call SNESSetDM(mechanical_snes,mechanical_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver
call DMCreateGlobalVector(mechanical_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs
call DMCreateGlobalVector(mechanical_mesh,solution_rate ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMCreateLocalVector (mechanical_mesh,solution_local ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
call DMSNESSetFunctionLocal(mechanical_mesh,FEM_mechanical_formResidual,PETSC_NULL_VEC,ierr) !< function to evaluate residual forces
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
2021-02-09 03:51:53 +05:30
call DMSNESSetJacobianLocal(mechanical_mesh,FEM_mechanical_formJacobian,PETSC_NULL_VEC,ierr) !< function to evaluate stiffness matrix
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
2021-02-09 03:51:53 +05:30
call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures
call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
2021-02-09 03:51:53 +05:30
call SNESSetTolerances(mechanical_snes,1.0,0.0,0.0,num%itmax,num%itmax,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
2021-02-09 03:51:53 +05:30
call SNESSetFromOptions(mechanical_snes,ierr); CHKERRQ(ierr)
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! init fields
2019-06-11 13:18:07 +05:30
call VecSet(solution ,0.0,ierr); CHKERRQ(ierr)
call VecSet(solution_rate ,0.0,ierr); CHKERRQ(ierr)
allocate(x_scal(cellDof))
allocate(nodalWeightsP(1))
allocate(nodalPointsP(dimPlex))
2019-06-11 13:18:07 +05:30
allocate(pv0(dimPlex))
allocate(pcellJ(dimPlex**2))
allocate(pinvcellJ(dimPlex**2))
2019-06-11 13:18:07 +05:30
allocate(cellJMat(dimPlex,dimPlex))
call PetscDSGetDiscretization(mechDS,0,mechFE,ierr)
CHKERRQ(ierr)
call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr)
2021-02-09 03:51:53 +05:30
call DMPlexGetHeightStratum(mechanical_mesh,0,cellStart,cellEnd,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0_pReal
2021-02-09 03:51:53 +05:30
call DMPlexComputeCellGeometryAffineFEM(mechanical_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
do basis = 0, nBasis*dimPlex-1, dimPlex
call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,ierr)
CHKERRQ(ierr)
call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,ierr)
CHKERRQ(ierr)
x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal)
2019-06-11 13:18:07 +05:30
enddo
px_scal => x_scal
2021-02-09 03:51:53 +05:30
call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr)
CHKERRQ(ierr)
enddo
2018-08-17 03:44:25 +05:30
2021-02-09 03:51:53 +05:30
end subroutine FEM_mechanical_init
2020-02-07 23:04:00 +05:30
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM load step
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
type(tSolutionState) function FEM_mechanical_solution( &
2018-08-17 03:44:25 +05:30
incInfoIn,timeinc,timeinc_old,fieldBC)
!--------------------------------------------------------------------------------------------------
! input data for solution
2019-06-11 13:18:07 +05:30
real(pReal), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment
type(tFieldBC), intent(in) :: &
fieldBC
character(len=*), intent(in) :: &
incInfoIn
PetscErrorCode :: ierr
2019-06-11 13:18:07 +05:30
SNESConvergedReason :: reason
2018-08-17 03:44:25 +05:30
2019-06-11 13:18:07 +05:30
incInfo = incInfoIn
2021-02-09 03:51:53 +05:30
FEM_mechanical_solution%converged =.false.
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
2019-06-11 13:18:07 +05:30
params%timeinc = timeinc
params%fieldBC = fieldBC
2021-02-09 03:51:53 +05:30
call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
2019-06-11 13:18:07 +05:30
terminallyIll = .false.
2019-06-11 13:18:07 +05:30
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
2021-02-09 03:51:53 +05:30
FEM_mechanical_solution%converged = .false.
FEM_mechanical_solution%iterationsNeeded = num%itmax
2019-06-11 13:18:07 +05:30
else ! >= 1 proper convergence (or terminally ill)
2021-02-09 03:51:53 +05:30
FEM_mechanical_solution%converged = .true.
call SNESGetIterationNumber(mechanical_snes,FEM_mechanical_solution%iterationsNeeded,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
endif
2020-09-19 11:50:29 +05:30
print'(/,a)', ' ==========================================================================='
2020-09-22 16:39:12 +05:30
flush(IO_STDOUT)
2018-08-17 03:44:25 +05:30
2021-02-09 03:51:53 +05:30
end function FEM_mechanical_solution
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief forms the FEM residual vector
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr)
2019-06-11 13:18:07 +05:30
DM :: dm_local
PetscObject,intent(in) :: dummy
PetscErrorCode :: ierr
2019-06-11 13:18:07 +05:30
PetscDS :: prob
Vec :: x_local, f_local, xx_local
PetscSection :: section
PetscScalar, dimension(:), pointer :: x_scal, pf_scal
PetscScalar, dimension(cellDof), target :: f_scal
PetscReal :: IcellJMat(dimPlex,dimPlex)
PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer
2019-06-11 13:18:07 +05:30
PetscInt :: cellStart, cellEnd, cell, field, face, &
2020-02-08 13:47:44 +05:30
qPt, basis, comp, cidx, &
numFields, &
bcSize,m
PetscReal :: detFAvg, detJ
PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat
2019-06-11 13:18:07 +05:30
IS :: bcPoints
2019-06-11 13:18:07 +05:30
allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2))
allocate(pinvcellJ(dimPlex**2))
2020-02-08 13:47:44 +05:30
allocate(x_scal(cellDof))
2020-02-08 12:49:06 +05:30
call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
CHKERRQ(ierr)
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0,xx_local,solution_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
2019-06-11 13:18:07 +05:30
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
enddo; enddo
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! evaluate field derivatives
do cell = cellStart, cellEnd-1 !< loop over all elements
2020-02-08 13:47:44 +05:30
call PetscSectionGetNumFields(section,numFields,ierr)
CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1
2019-06-11 13:18:07 +05:30
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: &
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo
enddo
homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo
2020-06-24 20:15:13 +05:30
if (num%BBarStabilisation) then
detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature))
do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1
homogenization_F(1:dimPlex,1:dimPlex,m) = homogenization_F(1:dimPlex,1:dimPlex,m) &
* (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0/real(dimPlex))
enddo
2019-06-11 13:18:07 +05:30
endif
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr)
enddo
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! evaluate constitutive response
2019-06-11 13:18:07 +05:30
call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,PETSC_COMM_WORLD,ierr)
ForwardData = .false.
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! integrating residual
do cell = cellStart, cellEnd-1 !< loop over all elements
2019-06-11 13:18:07 +05:30
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
f_scal = 0.0
do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1
2019-06-11 13:18:07 +05:30
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul(IcellJMat,basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: &
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo
enddo
f_scal = f_scal + &
matmul(transpose(BMat), &
reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), &
shape=[dimPlex*dimPlex]))*qWeights(qPt+1)
enddo
f_scal = f_scal*abs(detJ)
2019-06-11 13:18:07 +05:30
pf_scal => f_scal
call DMPlexVecSetClosure(dm_local,section,f_local,cell,pf_scal,ADD_VALUES,ierr)
CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
enddo
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
2021-02-09 03:51:53 +05:30
end subroutine FEM_mechanical_formResidual
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief forms the FEM stiffness matrix
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
2019-06-11 13:18:07 +05:30
DM :: dm_local
Mat :: Jac_pre, Jac
PetscObject, intent(in) :: dummy
PetscErrorCode :: ierr
2019-06-11 13:18:07 +05:30
PetscDS :: prob
Vec :: x_local, xx_local
2019-06-11 13:18:07 +05:30
PetscSection :: section, gSection
PetscReal, dimension(1, cellDof) :: MatB
PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA
PetscReal, dimension(3,3) :: F, FAvg, FInv
2019-06-11 13:18:07 +05:30
PetscReal :: detJ
PetscReal, dimension(:), pointer :: basisField, basisFieldDer, &
pV0, pCellJ, pInvcellJ
PetscScalar, dimension(:), pointer :: pK_e, x_scal
2020-01-21 11:38:02 +05:30
PetscScalar,dimension(cellDOF,cellDOF), target :: K_e
PetscScalar,dimension(cellDOF,cellDOF) :: K_eA , &
K_eB
PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx,bcSize, m
2019-06-11 13:18:07 +05:30
IS :: bcPoints
2019-06-11 13:18:07 +05:30
allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2))
allocate(pinvcellJ(dimPlex**2))
2019-06-11 13:18:07 +05:30
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)
call DMGetDS(dm_local,prob,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
2020-02-08 13:38:52 +05:30
call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,ierr); CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(x_local,section,0,field-1,bcPoints, &
0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc)
2019-06-11 13:18:07 +05:30
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
enddo; enddo
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements
2019-06-11 13:18:07 +05:30
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
2019-06-11 13:18:07 +05:30
CHKERRQ(ierr)
K_eA = 0.0
K_eB = 0.0
MatB = 0.0
FAvg = 0.0
BMatAvg = 0.0
do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt + 1
2019-06-11 13:18:07 +05:30
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
cidx = basis*dimPlex+comp
BMat(comp*dimPlex+1:(comp+1)*dimPlex,basis*dimPlex+comp+1) = &
matmul( reshape(pInvcellJ, shape = [dimPlex,dimPlex]),&
basisFieldDer((((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp )*dimPlex+1: &
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo
enddo
MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), &
2019-06-11 13:18:07 +05:30
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
2020-06-24 20:15:13 +05:30
if (num%BBarStabilisation) then
2019-06-11 13:18:07 +05:30
F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex])
FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
2019-06-11 13:18:07 +05:30
K_eB = K_eB - &
matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex*dimPlex,1]), &
2019-06-11 13:18:07 +05:30
matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA)
MatB = MatB &
+ matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1,dimPlex*dimPlex]),MatA)
FAvg = FAvg + F
2019-06-11 13:18:07 +05:30
BMatAvg = BMatAvg + BMat
else
K_eA = K_eA + matmul(transpose(BMat),MatA)
endif
enddo
2020-06-24 20:15:13 +05:30
if (num%BBarStabilisation) then
2019-06-11 13:18:07 +05:30
FInv = math_inv33(FAvg)
K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + &
(matmul(matmul(transpose(BMatAvg), &
reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + &
K_eB)/real(dimPlex)
2019-06-11 13:18:07 +05:30
else
K_e = K_eA
endif
2020-09-13 14:47:49 +05:30
K_e = (K_e + eps*math_eye(cellDof)) * abs(detJ)
2020-01-21 13:13:16 +05:30
#ifndef __INTEL_COMPILER
2020-01-21 11:38:02 +05:30
pK_e(1:cellDOF**2) => K_e
2020-01-21 13:13:16 +05:30
#else
! https://software.intel.com/en-us/forums/intel-fortran-compiler/topic/782230 (bug)
allocate(pK_e(cellDOF**2),source = reshape(K_e,[cellDOF**2]))
#endif
2019-06-11 13:18:07 +05:30
call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr)
CHKERRQ(ierr)
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr)
enddo
2019-06-11 13:18:07 +05:30
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)
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! apply boundary conditions
2020-10-06 10:42:04 +05:30
#if (PETSC_VERSION_MINOR < 14)
2019-06-11 13:18:07 +05:30
call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr)
2020-10-06 10:42:04 +05:30
#else
call DMPlexCreateRigidBody(dm_local,0,matnull,ierr); CHKERRQ(ierr)
#endif
2019-06-11 13:18:07 +05:30
call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatNullSpaceDestroy(matnull,ierr); CHKERRQ(ierr)
2018-08-17 03:44:25 +05:30
2021-02-09 03:51:53 +05:30
end subroutine FEM_mechanical_formJacobian
2018-08-17 03:44:25 +05:30
2019-05-12 12:59:51 +05:30
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
2019-06-11 13:18:07 +05:30
type(tFieldBC), intent(in) :: &
fieldBC
real(pReal), intent(in) :: &
2019-06-11 13:18:07 +05:30
timeinc_old, &
timeinc
logical, intent(in) :: &
2019-06-11 13:18:07 +05:30
guess
PetscInt :: field, face, bcSize
DM :: dm_local
Vec :: x_local
PetscSection :: section
IS :: bcPoints
PetscErrorCode :: ierr
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! forward last inc
if (guess .and. .not. cutBack) then
2019-06-11 13:18:07 +05:30
ForwardData = .True.
homogenization_F0 = homogenization_F
2021-02-09 03:51:53 +05:30
call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mechanical_snes into dm_local
2019-06-11 13:18:07 +05:30
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
2020-01-21 11:38:02 +05:30
call VecSet(x_local,0.0_pReal,ierr); CHKERRQ(ierr)
2019-06-11 13:18:07 +05:30
call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,ierr) !< retrieve my partition of global solution vector
CHKERRQ(ierr)
call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,ierr)
CHKERRQ(ierr)
call VecAXPY(solution_local,1.0,x_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,ierr)
if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,ierr)
CHKERRQ(ierr)
call utilities_projectBCValues(solution_local,section,0,field-1,bcPoints, &
0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old)
2019-06-11 13:18:07 +05:30
call ISDestroy(bcPoints,ierr); CHKERRQ(ierr)
endif
endif
enddo; enddo
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! update rate and forward last inc
2019-06-11 13:18:07 +05:30
call VecCopy(solution,solution_rate,ierr); CHKERRQ(ierr)
call VecScale(solution_rate,1.0/timeinc_old,ierr); CHKERRQ(ierr)
endif
call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr)
call VecScale(solution,timeinc,ierr); CHKERRQ(ierr)
2021-02-09 03:51:53 +05:30
end subroutine FEM_mechanical_forward
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief reporting
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,ierr)
2019-06-11 13:18:07 +05:30
SNES :: snes_local
PetscInt :: PETScIter
PetscReal :: xnorm,snorm,fnorm,divTol
SNESConvergedReason :: reason
PetscObject :: dummy
PetscErrorCode :: ierr
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
! report
2020-06-24 20:15:13 +05:30
divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*num%eps_struct_rtol,num%eps_struct_atol)
2019-06-11 13:18:07 +05:30
call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,ierr)
CHKERRQ(ierr)
if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN
2020-09-19 11:50:29 +05:30
print'(/,1x,a,a,i0,a,i0,f0.3)', trim(incInfo), &
2019-06-11 13:18:07 +05:30
' @ Iteration ',PETScIter,' mechanical residual norm = ', &
int(fnorm/divTol),fnorm/divTol-int(fnorm/divTol)
2020-12-15 03:07:14 +05:30
print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
' Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal
2020-09-22 16:39:12 +05:30
flush(IO_STDOUT)
2021-02-09 03:51:53 +05:30
end subroutine FEM_mechanical_converged
2018-08-17 03:44:25 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Calculate current coordinates (FEM nodal coordinates only at the moment)
!--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_updateCoords()
real(pReal), pointer, dimension(:) :: &
nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes)
real(pReal), pointer, dimension(:,:) :: &
nodeCoords !< nodal coordinates (3,Nnodes)
DM :: dm_local
Vec :: x_local
PetscErrorCode :: ierr
PetscInt :: dimPlex, pStart, pEnd, p, s, e
PetscSection :: section
call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr)
call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call DMGetDimension(dm_local,dimPlex,ierr); CHKERRQ(ierr)
call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,ierr); CHKERRQ(ierr)
allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal)
call VecGetArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
do p=pStart, pEnd-1
call DMPlexGetPointLocal(dm_local, p, s, e, ierr); CHKERRQ(ierr)
nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e)
end do
call discretization_setNodeCoords(nodeCoords)
call VecRestoreArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
end subroutine FEM_mechanical_updateCoords
2021-02-09 03:51:53 +05:30
end module mesh_mechanical_FEM