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

adopting to PETSc >3.10 

thanks to Matthew Knepley from the PETSc team
This commit is contained in:
Martin Diehl 2020-02-06 18:44:36 +01:00
parent cab000f4b4
commit 1ceba73d31
1 changed files with 82 additions and 75 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

157
src/mesh/FEM_mech.f90
View File

@ -13,7 +13,7 @@ module FEM_mech
use PETScDM use PETScDM
use PETScDMplex use PETScDMplex
use PETScDT use PETScDT
use prec use prec
use FEM_utilities use FEM_utilities
use mesh use mesh
@ -24,18 +24,18 @@ module FEM_mech
use FEsolving use FEsolving
use homogenization use homogenization
use math use math
implicit none implicit none
private private
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type tSolutionParams type tSolutionParams
type(tFieldBC) :: fieldBC type(tFieldBC) :: fieldBC
real(pReal) :: timeinc real(pReal) :: timeinc
real(pReal) :: timeincOld real(pReal) :: timeincOld
end type tSolutionParams end type tSolutionParams
type(tSolutionParams) :: params type(tSolutionParams) :: params
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -53,7 +53,7 @@ module FEM_mech
P_av = 0.0_pReal P_av = 0.0_pReal
logical :: ForwardData logical :: ForwardData
real(pReal), parameter :: eps = 1.0e-18_pReal real(pReal), parameter :: eps = 1.0e-18_pReal
public :: & public :: &
FEM_mech_init, & FEM_mech_init, &
FEM_mech_solution, & FEM_mech_solution, &
@ -75,12 +75,12 @@ subroutine FEM_mech_init(fieldBC)
PetscDualSpace :: mechDualSpace PetscDualSpace :: mechDualSpace
DMLabel, dimension(:),pointer :: pBCLabel DMLabel, dimension(:),pointer :: pBCLabel
DMLabel :: BCLabel DMLabel :: BCLabel
PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint
PetscInt :: numBC, bcSize, nc, & PetscInt :: numBC, bcSize, nc, &
field, faceSet, topologDim, nNodalPoints, & field, faceSet, topologDim, nNodalPoints, &
cellStart, cellEnd, cell, basis cellStart, cellEnd, cell, basis
IS :: bcPoint IS :: bcPoint
IS, dimension(:), pointer :: pBcComps, pBcPoints IS, dimension(:), pointer :: pBcComps, pBcPoints
PetscSection :: section PetscSection :: section
@ -89,13 +89,13 @@ subroutine FEM_mech_init(fieldBC)
nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ
PetscReal :: detJ PetscReal :: detJ
PetscReal, allocatable, target :: cellJMat(:,:) PetscReal, allocatable, target :: cellJMat(:,:)
PetscScalar, pointer :: px_scal(:) PetscScalar, pointer :: px_scal(:)
PetscScalar, allocatable, target :: x_scal(:) PetscScalar, allocatable, target :: x_scal(:)
character(len=*), parameter :: prefix = 'mechFE_' character(len=*), parameter :: prefix = 'mechFE_'
PetscErrorCode :: ierr PetscErrorCode :: ierr
write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>' write(6,'(/,a)') ' <<<+- FEM_mech init -+>>>'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -120,8 +120,15 @@ subroutine FEM_mech_init(fieldBC)
call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr)
call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr)
nBasis = nBasis/nc nBasis = nBasis/nc
#if (PETSC_VERSION_MINOR > 10)
call DMAddField(mech_mesh,PETSC_NULL_DMLABEL,mechFE,ierr); CHKERRQ(ierr)
call DMCreateDS(mech_mesh,ierr); CHKERRQ(ierr)
#endif
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
#if (PETSC_VERSION_MINOR < 11)
call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr) call DMGetDS(mech_mesh,mechDS,ierr); CHKERRQ(ierr)
call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr) call PetscDSAddDiscretization(mechDS,mechFE,ierr); CHKERRQ(ierr)
#endif
call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr) call PetscDSGetTotalDimension(mechDS,cellDof,ierr); CHKERRQ(ierr)
call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr) call PetscFEDestroy(mechFE,ierr); CHKERRQ(ierr)
call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr) call PetscQuadratureDestroy(mechQuad,ierr); CHKERRQ(ierr)
@ -171,8 +178,8 @@ subroutine FEM_mech_init(fieldBC)
else else
call ISCreateGeneral(PETSC_COMM_WORLD,0,[0],PETSC_COPY_VALUES,pbcPoints(numBC),ierr) call ISCreateGeneral(PETSC_COMM_WORLD,0,[0],PETSC_COPY_VALUES,pbcPoints(numBC),ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
endif endif
endif endif
enddo; enddo enddo; enddo
#if (PETSC_VERSION_MINOR < 11) #if (PETSC_VERSION_MINOR < 11)
call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, & call DMPlexCreateSection(mech_mesh,dimPlex,1,pNumComp,pNumDof, &
@ -188,11 +195,11 @@ subroutine FEM_mech_init(fieldBC)
do faceSet = 1, numBC do faceSet = 1, numBC
call ISDestroy(pbcPoints(faceSet),ierr); CHKERRQ(ierr) call ISDestroy(pbcPoints(faceSet),ierr); CHKERRQ(ierr)
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize solver specific parts of PETSc ! initialize solver specific parts of PETSc
call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr);CHKERRQ(ierr) call SNESCreate(PETSC_COMM_WORLD,mech_snes,ierr);CHKERRQ(ierr)
call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr) call SNESSetOptionsPrefix(mech_snes,'mech_',ierr);CHKERRQ(ierr)
call SNESSetDM(mech_snes,mech_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver call SNESSetDM(mech_snes,mech_mesh,ierr); CHKERRQ(ierr) !< set the mesh for non-linear solver
call DMCreateGlobalVector(mech_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs call DMCreateGlobalVector(mech_mesh,solution ,ierr); CHKERRQ(ierr) !< locally owned displacement Dofs
call DMCreateGlobalVector(mech_mesh,solution_rate ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step call DMCreateGlobalVector(mech_mesh,solution_rate ,ierr); CHKERRQ(ierr) !< locally owned velocity Dofs to guess solution at next load step
@ -201,13 +208,13 @@ subroutine FEM_mech_init(fieldBC)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMSNESSetJacobianLocal(mech_mesh,FEM_mech_formJacobian,PETSC_NULL_VEC,ierr) !< function to evaluate stiffness matrix call DMSNESSetJacobianLocal(mech_mesh,FEM_mech_formJacobian,PETSC_NULL_VEC,ierr) !< function to evaluate stiffness matrix
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures call SNESSetMaxLinearSolveFailures(mech_snes, huge(1), ierr); CHKERRQ(ierr) !< ignore linear solve failures
call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr) call SNESSetConvergenceTest(mech_snes,FEM_mech_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetTolerances(mech_snes,1.0,0.0,0.0,itmax,itmax,ierr) call SNESSetTolerances(mech_snes,1.0,0.0,0.0,itmax,itmax,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr) call SNESSetFromOptions(mech_snes,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init fields ! init fields
call VecSet(solution ,0.0,ierr); CHKERRQ(ierr) call VecSet(solution ,0.0,ierr); CHKERRQ(ierr)
@ -226,9 +233,9 @@ subroutine FEM_mech_init(fieldBC)
call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr) call PetscFEGetDualSpace(mechFE,mechDualSpace,ierr); CHKERRQ(ierr)
call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr) call DMPlexGetHeightStratum(mech_mesh,0,cellStart,cellEnd,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
x_scal = 0.0_pReal x_scal = 0.0_pReal
call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(mech_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex])
do basis = 0, nBasis*dimPlex-1, dimPlex do basis = 0, nBasis*dimPlex-1, dimPlex
@ -240,11 +247,11 @@ subroutine FEM_mech_init(fieldBC)
enddo enddo
px_scal => x_scal px_scal => x_scal
call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr) call DMPlexVecSetClosure(mech_mesh,section,solution_local,cell,px_scal,INSERT_ALL_VALUES,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
end subroutine FEM_mech_init end subroutine FEM_mech_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the FEM load step !> @brief solution for the FEM load step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -260,23 +267,23 @@ type(tSolutionState) function FEM_mech_solution( &
fieldBC fieldBC
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
PetscErrorCode :: ierr PetscErrorCode :: ierr
SNESConvergedReason :: reason SNESConvergedReason :: reason
incInfo = incInfoIn incInfo = incInfoIn
FEM_mech_solution%converged =.false. FEM_mech_solution%converged =.false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%timeinc = timeinc
params%timeincOld = timeinc_old params%timeincOld = timeinc_old
params%fieldBC = fieldBC params%fieldBC = fieldBC
call SNESSolve(mech_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mech_snes based on solution guess (result in solution) call SNESSolve(mech_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mech_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? call SNESGetConvergedReason(mech_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
terminallyIll = .false. terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
FEM_mech_solution%converged = .false. FEM_mech_solution%converged = .false.
FEM_mech_solution%iterationsNeeded = itmax FEM_mech_solution%iterationsNeeded = itmax
@ -285,9 +292,9 @@ type(tSolutionState) function FEM_mech_solution( &
call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr) call SNESGetIterationNumber(mech_snes,FEM_mech_solution%iterationsNeeded,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
endif endif
write(6,'(/,a)') ' ===========================================================================' write(6,'(/,a)') ' ==========================================================================='
flush(6) flush(6)
end function FEM_mech_solution end function FEM_mech_solution
@ -300,12 +307,12 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
DM :: dm_local DM :: dm_local
PetscObject,intent(in) :: dummy PetscObject,intent(in) :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscDS :: prob PetscDS :: prob
Vec :: x_local, f_local, xx_local Vec :: x_local, f_local, xx_local
PetscSection :: section PetscSection :: section
PetscScalar, dimension(:), pointer :: x_scal, pf_scal PetscScalar, dimension(:), pointer :: x_scal, pf_scal
PetscScalar, target :: f_scal(cellDof) PetscScalar, target :: f_scal(cellDof)
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex) PetscReal :: detJ, IcellJMat(dimPlex,dimPlex)
PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
@ -316,7 +323,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
PetscInt :: bcSize PetscInt :: bcSize
IS :: bcPoints IS :: bcPoints
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2)) allocate(pcellJ(dimPlex**2))
allocate(pinvcellJ(dimPlex**2)) allocate(pinvcellJ(dimPlex**2))
@ -342,10 +349,10 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate field derivatives ! evaluate field derivatives
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr) CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
do qPt = 0, nQuadrature-1 do qPt = 0, nQuadrature-1
@ -360,19 +367,19 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
enddo enddo
materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo enddo
if (BBarStabilisation) then if (BBarStabilisation) then
detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature))
do qPt = 1, nQuadrature do qPt = 1, nQuadrature
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = & materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = &
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1)* & materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1)* &
(detFAvg/math_det33(materialpoint_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex)) (detFAvg/math_det33(materialpoint_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex))
enddo enddo
endif endif
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr) call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
@ -381,11 +388,11 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
ForwardData = .false. ForwardData = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! integrating residual ! integrating residual
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr) CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
f_scal = 0.0 f_scal = 0.0
@ -402,17 +409,17 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
f_scal = f_scal + & f_scal = f_scal + &
matmul(transpose(BMat), & matmul(transpose(BMat), &
reshape(transpose(materialpoint_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), & reshape(transpose(materialpoint_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), &
shape=[dimPlex*dimPlex]))*qWeights(qPt+1) shape=[dimPlex*dimPlex]))*qWeights(qPt+1)
enddo enddo
f_scal = f_scal*abs(detJ) f_scal = f_scal*abs(detJ)
pf_scal => f_scal pf_scal => f_scal
call DMPlexVecSetClosure(dm_local,section,f_local,cell,pf_scal,ADD_VALUES,ierr) call DMPlexVecSetClosure(dm_local,section,f_local,cell,pf_scal,ADD_VALUES,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr) call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_formResidual end subroutine FEM_mech_formResidual
@ -426,12 +433,12 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
Mat :: Jac_pre, Jac Mat :: Jac_pre, Jac
PetscObject, intent(in) :: dummy PetscObject, intent(in) :: dummy
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscDS :: prob PetscDS :: prob
Vec :: x_local, xx_local Vec :: x_local, xx_local
PetscSection :: section, gSection PetscSection :: section, gSection
PetscReal, dimension(1, cellDof) :: MatB PetscReal, dimension(1, cellDof) :: MatB
PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA
PetscReal, dimension(3,3) :: F, FAvg, FInv PetscReal, dimension(3,3) :: F, FAvg, FInv
@ -442,19 +449,19 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
PetscScalar, dimension(:), pointer :: pK_e, x_scal PetscScalar, dimension(:), pointer :: pK_e, x_scal
PetscScalar,dimension(cellDOF,cellDOF), target :: K_e PetscScalar,dimension(cellDOF,cellDOF), target :: K_e
PetscScalar,dimension(cellDOF,cellDOF) :: K_eA , & PetscScalar,dimension(cellDOF,cellDOF) :: K_eA , &
K_eB K_eB
PetscInt :: cellStart, cellEnd, cell, field, face, & PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx,bcSize qPt, basis, comp, cidx,bcSize
IS :: bcPoints IS :: bcPoints
allocate(pV0(dimPlex)) allocate(pV0(dimPlex))
allocate(pcellJ(dimPlex**2)) allocate(pcellJ(dimPlex**2))
allocate(pinvcellJ(dimPlex**2)) allocate(pinvcellJ(dimPlex**2))
call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,ierr); CHKERRQ(ierr) 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 MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,ierr); CHKERRQ(ierr)
call MatZeroEntries(Jac,ierr); CHKERRQ(ierr) call MatZeroEntries(Jac,ierr); CHKERRQ(ierr)
@ -462,7 +469,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr) call PetscDSGetTabulation(prob,0,basisField,basisFieldDer,ierr)
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr) call DMGetGlobalSection(dm_local,gSection,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,ierr); CHKERRQ(ierr) call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,ierr); CHKERRQ(ierr)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries do field = 1, dimPlex; do face = 1, mesh_Nboundaries
@ -478,10 +485,10 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
endif endif
enddo; enddo enddo; enddo
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr) call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
do cell = cellStart, cellEnd-1 !< loop over all elements do cell = cellStart, cellEnd-1 !< loop over all elements
call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,ierr) !< get Dofs belonging to element
CHKERRQ(ierr) CHKERRQ(ierr)
call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
K_eA = 0.0 K_eA = 0.0
K_eB = 0.0 K_eB = 0.0
@ -501,33 +508,33 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
enddo enddo
MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
if (BBarStabilisation) then if (BBarStabilisation) then
F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex]) F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex])
FInv = math_inv33(F) FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
K_eB = K_eB - & K_eB = K_eB - &
matmul(transpose(matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), & matmul(transpose(matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex*dimPlex,1]), & shape=[dimPlex*dimPlex,1]), &
matmul(reshape(FInv(1:dimPlex,1:dimPlex), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA) shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA)
MatB = MatB + & MatB = MatB + &
matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA) matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA)
FAvg = FAvg + F FAvg = FAvg + F
BMatAvg = BMatAvg + BMat BMatAvg = BMatAvg + BMat
else else
K_eA = K_eA + matmul(transpose(BMat),MatA) K_eA = K_eA + matmul(transpose(BMat),MatA)
endif endif
enddo enddo
if (BBarStabilisation) then if (BBarStabilisation) then
FInv = math_inv33(FAvg) FInv = math_inv33(FAvg)
K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + & K_e = K_eA*math_det33(FAvg/real(nQuadrature))**(1.0/real(dimPlex)) + &
(matmul(matmul(transpose(BMatAvg), & (matmul(matmul(transpose(BMatAvg), &
reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + & reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex*dimPlex,1],order=[2,1])),MatB) + &
K_eB)/real(dimPlex) K_eB)/real(dimPlex)
else else
K_e = K_eA K_e = K_eA
endif endif
K_e = (K_e + eps*math_identity2nd(cellDof)) * abs(detJ) K_e = (K_e + eps*math_identity2nd(cellDof)) * abs(detJ)
#ifndef __INTEL_COMPILER #ifndef __INTEL_COMPILER
pK_e(1:cellDOF**2) => K_e pK_e(1:cellDOF**2) => K_e
@ -538,16 +545,16 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr) call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr) call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr) call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyEnd(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 MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,ierr); CHKERRQ(ierr)
call MatAssemblyEnd(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) call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! apply boundary conditions ! apply boundary conditions
call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr) call DMPlexCreateRigidBody(dm_local,matnull,ierr); CHKERRQ(ierr)
call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) call MatSetNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr) call MatSetNearNullSpace(Jac,matnull,ierr); CHKERRQ(ierr)
@ -565,7 +572,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
fieldBC fieldBC
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc_old, & timeinc_old, &
timeinc timeinc
logical, intent(in) :: & logical, intent(in) :: &
guess guess
@ -574,11 +581,11 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
Vec :: x_local Vec :: x_local
PetscSection :: section PetscSection :: section
IS :: bcPoints IS :: bcPoints
PetscErrorCode :: ierr PetscErrorCode :: ierr
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forward last inc ! forward last inc
if (guess .and. .not. cutBack) then if (guess .and. .not. cutBack) then
ForwardData = .True. ForwardData = .True.
materialpoint_F0 = materialpoint_F materialpoint_F0 = materialpoint_F
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local
@ -611,7 +618,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
endif endif
call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr) call VecCopy(solution_rate,solution,ierr); CHKERRQ(ierr)
call VecScale(solution,timeinc,ierr); CHKERRQ(ierr) call VecScale(solution,timeinc,ierr); CHKERRQ(ierr)
end subroutine FEM_mech_forward end subroutine FEM_mech_forward
@ -639,7 +646,7 @@ subroutine FEM_mech_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dumm
write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',& write(6,'(/,a,/,3(3(2x,f12.4,1x)/))',advance='no') ' Piola--Kirchhoff stress / MPa =',&
transpose(P_av)*1.e-6_pReal transpose(P_av)*1.e-6_pReal
flush(6) flush(6)
end subroutine FEM_mech_converged end subroutine FEM_mech_converged
end module FEM_mech end module FEM_mech