!-------------------------------------------------------------------------------------------------- !> @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 !-------------------------------------------------------------------------------------------------- module mesh_mechanical_FEM #include #include #include use PETScSNES use PETScDM use PETScDMplex use PETScDT #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) use MPI_f08 #endif use prec use FEM_utilities use discretization use discretization_mesh use config use IO use FEM_quadrature use homogenization use math #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY) implicit none(type,external) #else implicit none #endif private !-------------------------------------------------------------------------------------------------- ! derived types type tSolutionParams type(tFieldBC) :: fieldBC real(pReal) :: timeinc end type tSolutionParams type(tSolutionParams) :: params type, private :: tNumerics PetscInt :: & p_i, & !< integration order (quadrature rule) 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 !-------------------------------------------------------------------------------------------------- ! PETSc data SNES :: mechanical_snes Vec :: solution, solution_rate, solution_local PetscInt :: dimPlex, cellDof, nBasis integer :: nQuadrature PetscReal, allocatable, target :: qPoints(:), qWeights(:) MatNullSpace :: matnull !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. character(len=pStringLen) :: incInfo real(pReal), dimension(3,3) :: & P_av = 0.0_pReal logical :: ForwardData real(pReal), parameter :: eps = 1.0e-18_pReal external :: & ! ToDo: write interfaces PetscSectionGetNumFields, & PetscFESetQuadrature, & PetscFEGetDimension, & PetscFEDestroy, & PetscSectionGetDof, & PetscFEGetDualSpace, & PetscDualSpaceGetFunctional public :: & FEM_mechanical_init, & FEM_mechanical_solution, & FEM_mechanical_forward, & FEM_mechanical_updateCoords contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data !-------------------------------------------------------------------------------------------------- subroutine FEM_mechanical_init(fieldBC) type(tFieldBC), intent(in) :: fieldBC DM :: mechanical_mesh PetscFE :: mechFE PetscQuadrature :: mechQuad, functional PetscDS :: mechDS PetscDualSpace :: mechDualSpace DMLabel, dimension(:),pointer :: nolabel=> NULL() DMLabel :: BCLabel PetscInt, dimension(:), pointer :: pNumComp, pNumDof, pBcField, pBcPoint PetscInt :: numBC, bcSize, nc, & field, faceSet, topologDim, nNodalPoints, & cellStart, cellEnd, cell, basis IS :: bcPoint IS, dimension(:), pointer :: pBcComps, pBcPoints PetscSection :: section PetscReal, dimension(:), pointer :: qPointsP, qWeightsP, & nodalPointsP, nodalWeightsP,pV0, pCellJ, pInvcellJ PetscReal :: detJ PetscReal, allocatable, target :: cellJMat(:,:) PetscScalar, pointer, dimension(:) :: px_scal PetscScalar, allocatable, target, dimension(:) :: x_scal character(len=*), parameter :: prefix = 'mechFE_' PetscErrorCode :: err_PETSc real(pReal), dimension(3,3) :: devNull class(tNode), pointer :: & num_mesh print'(/,1x,a)', '<<<+- FEM_mech init -+>>>'; flush(IO_STDOUT) !----------------------------------------------------------------------------- ! read numerical parametes and do sanity checks num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) num%p_i = int(num_mesh%get_asInt('p_i',defaultVal = 2),pPETSCINT) num%itmax = int(num_mesh%get_asInt('itmax',defaultVal=250),pPETSCINT) 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) 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') !-------------------------------------------------------------------------------------------------- ! Setup FEM mech mesh call DMClone(geomMesh,mechanical_mesh,err_PETSc); CHKERRQ(err_PETSc) call DMGetDimension(mechanical_mesh,dimPlex,err_PETSc); CHKERRQ(err_PETSc) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech discretization qPoints = FEM_quadrature_points( dimPlex,num%p_i)%p qWeights = FEM_quadrature_weights(dimPlex,num%p_i)%p nQuadrature = FEM_nQuadrature( dimPlex,num%p_i) qPointsP => qPoints qWeightsP => qWeights call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,err_PETSc) CHKERRQ(err_PETSc) nc = dimPlex call PetscQuadratureSetData(mechQuad,dimPlex,nc,int(nQuadrature,pPETSCINT),qPointsP,qWeightsP,err_PETSc) CHKERRQ(err_PETSc) call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, & num%p_i,mechFE,err_PETSc); CHKERRQ(err_PETSc) call PetscFESetQuadrature(mechFE,mechQuad,err_PETSc); CHKERRQ(err_PETSc) call PetscFEGetDimension(mechFE,nBasis,err_PETSc); CHKERRQ(err_PETSc) nBasis = nBasis/nc call DMAddField(mechanical_mesh,PETSC_NULL_DMLABEL,mechFE,err_PETSc) CHKERRQ(err_PETSc) call DMCreateDS(mechanical_mesh,err_PETSc); CHKERRQ(err_PETSc) call DMGetDS(mechanical_mesh,mechDS,err_PETSc); CHKERRQ(err_PETSc) call PetscDSGetTotalDimension(mechDS,cellDof,err_PETSc); CHKERRQ(err_PETSc) call PetscFEDestroy(mechFE,err_PETSc); CHKERRQ(err_PETSc) call PetscQuadratureDestroy(mechQuad,err_PETSc); CHKERRQ(err_PETSc) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech boundary conditions call DMGetLabel(mechanical_mesh,'Face Sets',BCLabel,err_PETSc) CHKERRQ(err_PETSc) call DMPlexLabelComplete(mechanical_mesh,BCLabel,err_PETSc); CHKERRQ(err_PETSc) call DMGetLocalSection(mechanical_mesh,section,err_PETSc); CHKERRQ(err_PETSc) allocate(pnumComp(1), source=dimPlex) allocate(pnumDof(0:dimPlex), source = 0_pPETSCINT) do topologDim = 0, dimPlex call DMPlexGetDepthStratum(mechanical_mesh,topologDim,cellStart,cellEnd,err_PETSc) CHKERRQ(err_PETSc) call PetscSectionGetDof(section,cellStart,pnumDof(topologDim),err_PETSc) CHKERRQ(err_PETSc) 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_pPETSCINT) 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_pPETSCINT,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc) CHKERRQ(err_PETSc) call DMGetStratumSize(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcSize,err_PETSc) CHKERRQ(err_PETSc) if (bcSize > 0) then call DMGetStratumIS(mechanical_mesh,'Face Sets',mesh_boundaries(faceSet),bcPoint,err_PETSc) CHKERRQ(err_PETSc) call ISGetIndicesF90(bcPoint,pBcPoint,err_PETSc); CHKERRQ(err_PETSc) call ISCreateGeneral(PETSC_COMM_WORLD,bcSize,pBcPoint,PETSC_COPY_VALUES,pbcPoints(numBC),err_PETSc) CHKERRQ(err_PETSc) call ISRestoreIndicesF90(bcPoint,pBcPoint,err_PETSc); CHKERRQ(err_PETSc) call ISDestroy(bcPoint,err_PETSc); CHKERRQ(err_PETSc) else call ISCreateGeneral(PETSC_COMM_WORLD,0_pPETSCINT,[0_pPETSCINT],PETSC_COPY_VALUES,pbcPoints(numBC),err_PETSc) CHKERRQ(err_PETSc) endif endif enddo; enddo call DMPlexCreateSection(mechanical_mesh,nolabel,pNumComp,pNumDof, & numBC,pBcField,pBcComps,pBcPoints,PETSC_NULL_IS,section,err_PETSc) CHKERRQ(err_PETSc) call DMSetSection(mechanical_mesh,section,err_PETSc); CHKERRQ(err_PETSc) do faceSet = 1, numBC call ISDestroy(pbcPoints(faceSet),err_PETSc); CHKERRQ(err_PETSc) enddo !-------------------------------------------------------------------------------------------------- ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,mechanical_snes,err_PETSc);CHKERRQ(err_PETSc) call SNESSetOptionsPrefix(mechanical_snes,'mechanical_',err_PETSc) CHKERRQ(err_PETSc) call SNESSetDM(mechanical_snes,mechanical_mesh,err_PETSc) ! set the mesh for non-linear solver CHKERRQ(err_PETSc) call DMCreateGlobalVector(mechanical_mesh,solution, err_PETSc) ! locally owned displacement Dofs CHKERRQ(err_PETSc) call DMCreateGlobalVector(mechanical_mesh,solution_rate, err_PETSc) ! locally owned velocity Dofs to guess solution at next load step CHKERRQ(err_PETSc) call DMCreateLocalVector (mechanical_mesh,solution_local,err_PETSc) ! locally owned velocity Dofs to guess solution at next load step CHKERRQ(err_PETSc) call DMSNESSetFunctionLocal(mechanical_mesh,FEM_mechanical_formResidual,PETSC_NULL_VEC,err_PETSc) ! function to evaluate residual forces CHKERRQ(err_PETSc) call DMSNESSetJacobianLocal(mechanical_mesh,FEM_mechanical_formJacobian,PETSC_NULL_VEC,err_PETSc) ! function to evaluate stiffness matrix CHKERRQ(err_PETSc) call SNESSetMaxLinearSolveFailures(mechanical_snes, huge(1_pPETSCINT), err_PETSc) ! ignore linear solve failures CHKERRQ(err_PETSc) call SNESSetConvergenceTest(mechanical_snes,FEM_mechanical_converged,PETSC_NULL_VEC,PETSC_NULL_FUNCTION,err_PETSc) CHKERRQ(err_PETSc) call SNESSetTolerances(mechanical_snes,1.0_pReal,0.0_pReal,0.0_pReal,num%itmax,num%itmax,err_PETSc) CHKERRQ(err_PETSc) call SNESSetFromOptions(mechanical_snes,err_PETSc); CHKERRQ(err_PETSc) !-------------------------------------------------------------------------------------------------- ! init fields call VecSet(solution ,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc) call VecSet(solution_rate,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc) allocate(x_scal(cellDof)) allocate(nodalWeightsP(1)) allocate(nodalPointsP(dimPlex)) allocate(pv0(dimPlex)) allocate(pcellJ(dimPlex**2)) allocate(pinvcellJ(dimPlex**2)) allocate(cellJMat(dimPlex,dimPlex)) call PetscDSGetDiscretization(mechDS,0_pPETSCINT,mechFE,err_PETSc) CHKERRQ(err_PETSc) call PetscFEGetDualSpace(mechFE,mechDualSpace,err_PETSc); CHKERRQ(err_PETSc) call DMPlexGetHeightStratum(mechanical_mesh,0_pPETSCINT,cellStart,cellEnd,err_PETSc) CHKERRQ(err_PETSc) do cell = cellStart, cellEnd-1 !< loop over all elements x_scal = 0.0_pReal call DMPlexComputeCellGeometryAffineFEM(mechanical_mesh,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc) CHKERRQ(err_PETSc) cellJMat = reshape(pCellJ,shape=[dimPlex,dimPlex]) do basis = 0, nBasis*dimPlex-1, dimPlex call PetscDualSpaceGetFunctional(mechDualSpace,basis,functional,err_PETSc) CHKERRQ(err_PETSc) call PetscQuadratureGetData(functional,dimPlex,nc,nNodalPoints,nodalPointsP,nodalWeightsP,err_PETSc) CHKERRQ(err_PETSc) x_scal(basis+1:basis+dimPlex) = pV0 + matmul(transpose(cellJMat),nodalPointsP + 1.0_pReal) enddo px_scal => x_scal call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,err_PETSc) CHKERRQ(err_PETSc) enddo call utilities_constitutiveResponse(0.0_pReal,devNull,.true.) end subroutine FEM_mechanical_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the FEM load step !-------------------------------------------------------------------------------------------------- type(tSolutionState) function FEM_mechanical_solution( & incInfoIn,timeinc,timeinc_old,fieldBC) !-------------------------------------------------------------------------------------------------- ! input data for solution 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 :: err_PETSc SNESConvergedReason :: reason incInfo = incInfoIn FEM_mechanical_solution%converged =.false. !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%timeinc = timeinc params%fieldBC = fieldBC call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,err_PETSc) ! solve mechanical_snes based on solution guess (result in solution) CHKERRQ(err_PETSc) call SNESGetConvergedReason(mechanical_snes,reason,err_PETSc) ! solution converged? CHKERRQ(err_PETSc) terminallyIll = .false. if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error FEM_mechanical_solution%converged = .false. FEM_mechanical_solution%iterationsNeeded = num%itmax else ! >= 1 proper convergence (or terminally ill) FEM_mechanical_solution%converged = .true. call SNESGetIterationNumber(mechanical_snes,FEM_mechanical_solution%iterationsNeeded,err_PETSc) CHKERRQ(err_PETSc) endif print'(/,1x,a)', '===========================================================================' flush(IO_STDOUT) end function FEM_mechanical_solution !-------------------------------------------------------------------------------------------------- !> @brief forms the FEM residual vector !-------------------------------------------------------------------------------------------------- subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc) DM :: dm_local PetscObject,intent(in) :: dummy PetscErrorCode :: err_PETSc integer(MPI_INTEGER_KIND) :: err_MPI 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 PetscInt :: cellStart, cellEnd, cell, field, face, & qPt, basis, comp, cidx, & numFields, & bcSize,m,i PetscReal :: detFAvg, detJ PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat IS :: bcPoints allocate(pV0(dimPlex)) allocate(pcellJ(dimPlex**2)) allocate(pinvcellJ(dimPlex**2)) allocate(x_scal(cellDof)) call DMGetLocalSection(dm_local,section,err_PETSc); CHKERRQ(err_PETSc) call DMGetDS(dm_local,prob,err_PETSc); CHKERRQ(err_PETSc) call PetscDSGetTabulation(prob,0_pPETSCINT,basisField,basisFieldDer,err_PETSc) CHKERRQ(err_PETSc) call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc) CHKERRQ(err_PETSc) call DMGetLocalVector(dm_local,x_local,err_PETSc) CHKERRQ(err_PETSc) call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,err_PETSc) CHKERRQ(err_PETSc) 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,err_PETSc) if (bcSize > 0) then call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) CHKERRQ(err_PETSc) call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc) call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc) endif endif enddo; enddo !-------------------------------------------------------------------------------------------------- ! evaluate field derivatives do cell = cellStart, cellEnd-1 !< loop over all elements call PetscSectionGetNumFields(section,numFields,err_PETSc) CHKERRQ(err_PETSc) call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,err_PETSc) !< get Dofs belonging to element CHKERRQ(err_PETSc) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc) CHKERRQ(err_PETSc) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) do qPt = 0, nQuadrature-1 m = cell*nQuadrature + qPt+1 BMat = 0.0_pReal do basis = 0, nBasis-1 do comp = 0, dimPlex-1 cidx = basis*dimPlex+comp i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = & matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex)) enddo enddo homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) enddo if (num%BBarStabilisation) then detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature,pReal)) 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_pReal/real(dimPlex,pReal)) enddo endif call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,err_PETSc) CHKERRQ(err_PETSc) enddo !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' ForwardData = .false. !-------------------------------------------------------------------------------------------------- ! integrating residual do cell = cellStart, cellEnd-1 !< loop over all elements call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,err_PETSc) !< get Dofs belonging to element CHKERRQ(err_PETSc) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc) CHKERRQ(err_PETSc) IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex]) f_scal = 0.0_pReal do qPt = 0, nQuadrature-1 m = cell*nQuadrature + qPt+1 BMat = 0.0_pReal do basis = 0, nBasis-1 do comp = 0, dimPlex-1 cidx = basis*dimPlex+comp i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = & matmul(IcellJMat,basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*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_pPETSCINT) enddo f_scal = f_scal*abs(detJ) pf_scal => f_scal call DMPlexVecSetClosure(dm_local,section,f_local,cell,pf_scal,ADD_VALUES,err_PETSc) CHKERRQ(err_PETSc) call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,err_PETSc) CHKERRQ(err_PETSc) enddo call DMRestoreLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc) end subroutine FEM_mechanical_formResidual !-------------------------------------------------------------------------------------------------- !> @brief forms the FEM stiffness matrix !-------------------------------------------------------------------------------------------------- subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_PETSc) DM :: dm_local Mat :: Jac_pre, Jac PetscObject, intent(in) :: dummy PetscErrorCode :: err_PETSc PetscDS :: prob Vec :: x_local, xx_local PetscSection :: section, gSection PetscReal, dimension(1, cellDof) :: MatB PetscReal, dimension(dimPlex**2,cellDof) :: BMat, BMatAvg, MatA PetscReal, dimension(3,3) :: F, FAvg, FInv PetscReal :: detJ PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & pV0, pCellJ, pInvcellJ PetscScalar, dimension(:), pointer :: pK_e, x_scal 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, i IS :: bcPoints allocate(pV0(dimPlex)) allocate(pcellJ(dimPlex**2)) allocate(pinvcellJ(dimPlex**2)) call MatSetOption(Jac,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE,err_PETSc) CHKERRQ(err_PETSc) call MatSetOption(Jac,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE,err_PETSc) CHKERRQ(err_PETSc) call MatZeroEntries(Jac,err_PETSc) CHKERRQ(err_PETSc) call DMGetDS(dm_local,prob,err_PETSc) CHKERRQ(err_PETSc) call PetscDSGetTabulation(prob,0_pPETSCINT,basisField,basisFieldDer,err_PETSc) call DMGetLocalSection(dm_local,section,err_PETSc) CHKERRQ(err_PETSc) call DMGetGlobalSection(dm_local,gSection,err_PETSc) CHKERRQ(err_PETSc) call DMGetLocalVector(dm_local,x_local,err_PETSc) CHKERRQ(err_PETSc) call VecWAXPY(x_local,1.0_pReal,xx_local,solution_local,err_PETSc) CHKERRQ(err_PETSc) 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,err_PETSc) if (bcSize > 0) then call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) CHKERRQ(err_PETSc) call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & 0.0_pReal,params%fieldBC%componentBC(field)%Value(face),params%timeinc) call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc) endif endif enddo; enddo call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc) CHKERRQ(err_PETSc) do cell = cellStart, cellEnd-1 !< loop over all elements call DMPlexVecGetClosure(dm_local,section,x_local,cell,x_scal,err_PETSc) !< get Dofs belonging to element CHKERRQ(err_PETSc) call DMPlexComputeCellGeometryAffineFEM(dm_local,cell,pV0,pCellJ,pInvcellJ,detJ,err_PETSc) CHKERRQ(err_PETSc) K_eA = 0.0_pReal K_eB = 0.0_pReal MatB = 0.0_pReal FAvg = 0.0_pReal BMatAvg = 0.0_pReal do qPt = 0, nQuadrature-1 m = cell*nQuadrature + qPt + 1 BMat = 0.0_pReal do basis = 0, nBasis-1 do comp = 0, dimPlex-1 cidx = basis*dimPlex+comp i = ((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp BMat(comp*dimPlex+1_pPETSCINT:(comp+1_pPETSCINT)*dimPlex,basis*dimPlex+comp+1_pPETSCINT) = & matmul(reshape(pInvcellJ,[dimPlex,dimPlex]),basisFieldDer(i*dimPlex+1_pPETSCINT:(i+1_pPETSCINT)*dimPlex)) enddo enddo MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT) if (num%BBarStabilisation) then 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_pReal/real(dimPlex,pReal)) K_eB = K_eB - & matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex**2,1_pPETSCINT]), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), & shape=[1_pPETSCINT,dimPlex**2],order=[2,1]),BMat))),MatA) MatB = MatB & + matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1_pPETSCINT,dimPlex**2]),MatA) FAvg = FAvg + F BMatAvg = BMatAvg + BMat else K_eA = K_eA + matmul(transpose(BMat),MatA) endif enddo if (num%BBarStabilisation) then FInv = math_inv33(FAvg) K_e = K_eA*math_det33(FAvg/real(nQuadrature,pReal))**(1.0_pReal/real(dimPlex,pReal)) + & (matmul(matmul(transpose(BMatAvg), & reshape(FInv(1:dimPlex,1:dimPlex),shape=[dimPlex**2,1_pPETSCINT],order=[2,1])),MatB) + & K_eB)/real(dimPlex,pReal) else K_e = K_eA endif K_e = (K_e + eps*math_eye(int(cellDof))) * abs(detJ) #ifndef __INTEL_COMPILER pK_e(1:cellDOF**2) => K_e #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 call DMPlexMatSetClosure(dm_local,section,gSection,Jac,cell,pK_e,ADD_VALUES,err_PETSc) CHKERRQ(err_PETSc) call DMPlexVecRestoreClosure(dm_local,section,x_local,cell,x_scal,err_PETSc) CHKERRQ(err_PETSc) enddo call MatAssemblyBegin(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc) call MatAssemblyEnd(Jac,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc) call MatAssemblyBegin(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc) call MatAssemblyEnd(Jac_pre,MAT_FINAL_ASSEMBLY,err_PETSc); CHKERRQ(err_PETSc) call DMRestoreLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc) !-------------------------------------------------------------------------------------------------- ! apply boundary conditions #if (PETSC_VERSION_MINOR < 14) call DMPlexCreateRigidBody(dm_local,matnull,err_PETSc); CHKERRQ(err_PETSc) #else call DMPlexCreateRigidBody(dm_local,0_pPETSCINT,matnull,err_PETSc) CHKERRQ(err_PETSc) #endif call MatSetNullSpace(Jac,matnull,err_PETSc); CHKERRQ(err_PETSc) call MatSetNearNullSpace(Jac,matnull,err_PETSc); CHKERRQ(err_PETSc) call MatNullSpaceDestroy(matnull,err_PETSc); CHKERRQ(err_PETSc) end subroutine FEM_mechanical_formJacobian !-------------------------------------------------------------------------------------------------- !> @brief forwarding routine !-------------------------------------------------------------------------------------------------- subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) type(tFieldBC), intent(in) :: & fieldBC real(pReal), intent(in) :: & timeinc_old, & timeinc logical, intent(in) :: & guess PetscInt :: field, face, bcSize DM :: dm_local Vec :: x_local PetscSection :: section IS :: bcPoints PetscErrorCode :: err_PETSc !-------------------------------------------------------------------------------------------------- ! forward last inc if (guess .and. .not. cutBack) then ForwardData = .True. homogenization_F0 = homogenization_F call SNESGetDM(mechanical_snes,dm_local,err_PETSc); CHKERRQ(err_PETSc) !< retrieve mesh info from mechanical_snes into dm_local call DMGetSection(dm_local,section,err_PETSc); CHKERRQ(err_PETSc) call DMGetLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc) call VecSet(x_local,0.0_pReal,err_PETSc); CHKERRQ(err_PETSc) call DMGlobalToLocalBegin(dm_local,solution,INSERT_VALUES,x_local,err_PETSc) !< retrieve my partition of global solution vector CHKERRQ(err_PETSc) call DMGlobalToLocalEnd(dm_local,solution,INSERT_VALUES,x_local,err_PETSc) CHKERRQ(err_PETSc) call VecAXPY(solution_local,1.0_pReal,x_local,err_PETSc); CHKERRQ(err_PETSc) 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,err_PETSc) if (bcSize > 0) then call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) CHKERRQ(err_PETSc) call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, & 0.0_pReal,fieldBC%componentBC(field)%Value(face),timeinc_old) call ISDestroy(bcPoints,err_PETSc); CHKERRQ(err_PETSc) endif endif enddo; enddo call DMRestoreLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc) !-------------------------------------------------------------------------------------------------- ! update rate and forward last inc call VecCopy(solution,solution_rate,err_PETSc); CHKERRQ(err_PETSc) call VecScale(solution_rate,timeinc_old**(-1),err_PETSc); CHKERRQ(err_PETSc) endif call VecCopy(solution_rate,solution,err_PETSc); CHKERRQ(err_PETSc) call VecScale(solution,timeinc,err_PETSc); CHKERRQ(err_PETSc) end subroutine FEM_mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief reporting !-------------------------------------------------------------------------------------------------- subroutine FEM_mechanical_converged(snes_local,PETScIter,xnorm,snorm,fnorm,reason,dummy,err_PETSc) SNES :: snes_local PetscInt :: PETScIter PetscReal :: xnorm,snorm,fnorm,divTol SNESConvergedReason :: reason PetscObject :: dummy PetscErrorCode :: err_PETSc !-------------------------------------------------------------------------------------------------- ! report divTol = max(maxval(abs(P_av(1:dimPlex,1:dimPlex)))*num%eps_struct_rtol,num%eps_struct_atol) call SNESConvergedDefault(snes_local,PETScIter,xnorm,snorm,fnorm/divTol,reason,dummy,err_PETSc) CHKERRQ(err_PETSc) if (terminallyIll) reason = SNES_DIVERGED_FUNCTION_DOMAIN print'(/,1x,a,a,i0,a,f0.3)', trim(incInfo), & ' @ Iteration ',PETScIter,' mechanical residual norm = ',fnorm/divTol print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & 'Piola--Kirchhoff stress / MPa =',transpose(P_av)*1.e-6_pReal flush(IO_STDOUT) end subroutine FEM_mechanical_converged !-------------------------------------------------------------------------------------------------- !> @brief Calculate current coordinates (both nodal and ip coordinates) !-------------------------------------------------------------------------------------------------- subroutine FEM_mechanical_updateCoords() PetscReal, pointer, dimension(:,:) :: & nodeCoords !< nodal coordinates (3,Nnodes) real(pReal), pointer, dimension(:,:,:) :: & ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems) integer :: & qPt, & comp, & qOffset, & nOffset DM :: dm_local Vec :: x_local PetscErrorCode :: err_PETSc PetscInt :: pStart, pEnd, p, s, e, q, & cellStart, cellEnd, c, n PetscSection :: section PetscQuadrature :: mechQuad PetscReal, dimension(:), pointer :: basisField, basisFieldDer, & nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes) PetscScalar, dimension(:), pointer :: x_scal call SNESGetDM(mechanical_snes,dm_local,err_PETSc); CHKERRQ(err_PETSc) call DMGetDS(dm_local,mechQuad,err_PETSc); CHKERRQ(err_PETSc) call DMGetLocalSection(dm_local,section,err_PETSc); CHKERRQ(err_PETSc) call DMGetLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc) call DMGetDimension(dm_local,dimPlex,err_PETSc); CHKERRQ(err_PETSc) ! write cell vertex displacements call DMPlexGetDepthStratum(dm_local,0_pPETSCINT,pStart,pEnd,err_PETSc) CHKERRQ(err_PETSc) allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal) call VecGetArrayF90(x_local,nodeCoords_linear,err_PETSc); CHKERRQ(err_PETSc) do p=pStart, pEnd-1 call DMPlexGetPointLocal(dm_local, p, s, e, err_PETSc); CHKERRQ(err_PETSc) nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e) enddo call discretization_setNodeCoords(nodeCoords) call VecRestoreArrayF90(x_local,nodeCoords_linear,err_PETSc); CHKERRQ(err_PETSc) ! write ip displacements call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc) CHKERRQ(err_PETSc) call PetscDSGetTabulation(mechQuad,0_pPETSCINT,basisField,basisFieldDer,err_PETSc) CHKERRQ(err_PETSc) allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal) do c=cellStart,cellEnd-1 qOffset=0 call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,err_PETSc) !< get nodal coordinates of each element CHKERRQ(err_PETSc) do qPt=0,nQuadrature-1 qOffset= qPt * (size(basisField)/nQuadrature) do comp=0,dimPlex-1 !< loop over components nOffset=0 q = comp do n=0,nBasis-1 ipCoords(comp+1,qPt+1,c+1)=ipCoords(comp+1,qPt+1,c+1)+& sum(basisField(qOffset+(q*dimPlex)+1:qOffset+(q*dimPlex)+dimPlex)*& x_scal(nOffset+1:nOffset+dimPlex)) q = q+dimPlex nOffset = nOffset+dimPlex enddo enddo enddo call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,err_PETSc) CHKERRQ(err_PETSc) end do call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature])) call DMRestoreLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc) end subroutine FEM_mechanical_updateCoords end module mesh_mechanical_FEM