PETSc 3.18.0 compatibility

CHKERR macros are now longer, should be replaced with PetscCall in the
future

CMakeLists.txt

@@ -11,7 +11,7 @@ endif()
 project(Prerequisites LANGUAGES)
 set(ENV{PKG_CONFIG_PATH} "$ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig:$ENV{PKG_CONFIG_PATH}")
 pkg_check_modules(PETSC_MIN REQUIRED PETSc>=3.12.0 QUIET) #CMake does not support version range
-pkg_check_modules(PETSC REQUIRED PETSc<3.18.0)
+pkg_check_modules(PETSC REQUIRED PETSc<3.19.0)
 
 pkg_get_variable(CMAKE_Fortran_COMPILER PETSc fcompiler)
 pkg_get_variable(CMAKE_C_COMPILER PETSc ccompiler)

src/CLI.f90

@@ -7,7 +7,7 @@
 !--------------------------------------------------------------------------------------------------
 #define PETSC_MAJOR 3
 #define PETSC_MINOR_MIN 12
-#define PETSC_MINOR_MAX 17
+#define PETSC_MINOR_MAX 18
 
 module CLI
   use, intrinsic :: ISO_fortran_env

src/CMakeLists.txt

@@ -1,7 +1,7 @@
 # special flags for some files
 if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
     # long lines for interaction matrix
-    SET_SOURCE_FILES_PROPERTIES("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240")
+    set_source_files_properties("lattice.f90" PROPERTIES COMPILE_FLAGS "-ffree-line-length-240")
 endif()
 
 file(GLOB damask-sources CONFIGURE_DEPENDS *.f90 *.c)
@@ -14,6 +14,15 @@ elseif(PROJECT_NAME STREQUAL "damask-mesh")
     file(GLOB solver-sources CONFIGURE_DEPENDS mesh/*.f90)
 endif()
 
+foreach(solver-source ${solver-sources})
+    file(READ ${solver-source} content)
+    string(FIND "${content}" "CHKERR" found)
+    if(NOT ${found} EQUAL -1)
+      set_source_files_properties(${solver-source} PROPERTIES COMPILE_FLAGS "-ffree-line-length-160")
+    endif()
+endforeach()
+
+
 if(NOT CMAKE_BUILD_TYPE STREQUAL "SYNTAXONLY")
     add_executable(${executable-name} ${damask-sources} ${solver-sources})
     install(TARGETS ${executable-name} RUNTIME DESTINATION bin)

src/mesh/FEM_utilities.f90

@@ -192,7 +192,8 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
   call ISGetSize(bcPointsIS,nBcPoints,err_PETSc)
   CHKERRQ(err_PETSc)
   if (nBcPoints > 0) call ISGetIndicesF90(bcPointsIS,bcPoints,err_PETSc)
-  call VecGetArrayF90(localVec,localArray,err_PETSc); CHKERRQ(err_PETSc)
+  call VecGetArrayF90(localVec,localArray,err_PETSc)
+  CHKERRQ(err_PETSc)
   do point = 1, nBcPoints
     call PetscSectionGetFieldDof(section,bcPoints(point),field,numDof,err_PETSc)
     CHKERRQ(err_PETSc)
@@ -202,9 +203,12 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
       localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
     end do
   end do
-  call VecRestoreArrayF90(localVec,localArray,err_PETSc); CHKERRQ(err_PETSc)
-  call VecAssemblyBegin(localVec, err_PETSc); CHKERRQ(err_PETSc)
-  call VecAssemblyEnd  (localVec, err_PETSc); CHKERRQ(err_PETSc)
+  call VecRestoreArrayF90(localVec,localArray,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call VecAssemblyBegin(localVec, err_PETSc)
+  CHKERRQ(err_PETSc)
+  call VecAssemblyEnd  (localVec, err_PETSc)
+  CHKERRQ(err_PETSc)
   if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,err_PETSc)
   CHKERRQ(err_PETSc)
 

src/mesh/discretization_mesh.f90

@@ -158,7 +158,8 @@ subroutine discretization_mesh_init(restart)
   call MPI_Bcast(mesh_boundaries,mesh_Nboundaries,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
   if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
 
-  call DMDestroy(globalMesh,err_PETSc); CHKERRQ(err_PETSc)
+  call DMDestroy(globalMesh,err_PETSc)
+  CHKERRQ(err_PETSc)
 
   call DMGetStratumSize(geomMesh,'depth',dimPlex,Nelems,err_PETSc)
   CHKERRQ(err_PETSc)

src/mesh/mesh_mech_FEM.f90

@@ -146,8 +146,10 @@ subroutine FEM_mechanical_init(fieldBC)
 
 !--------------------------------------------------------------------------------------------------
 ! Setup FEM mech mesh
-  call DMClone(geomMesh,mechanical_mesh,err_PETSc); CHKERRQ(err_PETSc)
-  call DMGetDimension(mechanical_mesh,dimPlex,err_PETSc); CHKERRQ(err_PETSc)
+  call DMClone(geomMesh,mechanical_mesh,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call DMGetDimension(mechanical_mesh,dimPlex,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 !--------------------------------------------------------------------------------------------------
 ! Setup FEM mech discretization
@@ -162,24 +164,34 @@ subroutine FEM_mechanical_init(fieldBC)
   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)
+                            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)
+  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)
+  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
@@ -206,11 +218,14 @@ subroutine FEM_mechanical_init(fieldBC)
       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 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)
+        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)
@@ -220,14 +235,17 @@ subroutine FEM_mechanical_init(fieldBC)
   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)
+  call DMSetSection(mechanical_mesh,section,err_PETSc)
+  CHKERRQ(err_PETSc)
   do faceSet = 1, numBC
-    call ISDestroy(pbcPoints(faceSet),err_PETSc); CHKERRQ(err_PETSc)
+    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 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
@@ -248,12 +266,15 @@ subroutine FEM_mechanical_init(fieldBC)
   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)
+  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)
+  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))
@@ -263,7 +284,8 @@ subroutine FEM_mechanical_init(fieldBC)
   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 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
@@ -365,8 +387,10 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
   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 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)
@@ -383,7 +407,8 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,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)
+        call ISDestroy(bcPoints,err_PETSc)
+        CHKERRQ(err_PETSc)
       endif
     endif
   enddo; enddo
@@ -464,7 +489,8 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,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)
+  call DMRestoreLocalVector(dm_local,x_local,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 end subroutine FEM_mechanical_formResidual
 
@@ -531,7 +557,8 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
         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)
+        call ISDestroy(bcPoints,err_PETSc)
+        CHKERRQ(err_PETSc)
       endif
     endif
   enddo; enddo
@@ -598,23 +625,32 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
     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)
+  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)
+  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)
+  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
 
@@ -644,15 +680,20 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
   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 SNESGetDM(mechanical_snes,dm_local,err_PETSc)                                              !< retrieve mesh info from mechanical_snes into dm_local
+    CHKERRQ(err_PETSc)
+    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)
+    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)
@@ -661,19 +702,25 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
           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)
+          call ISDestroy(bcPoints,err_PETSc)
+          CHKERRQ(err_PETSc)
         endif
       endif
     enddo; enddo
-    call DMRestoreLocalVector(dm_local,x_local,err_PETSc); CHKERRQ(err_PETSc)
+    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)
+    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)
+  call VecCopy(solution_rate,solution,err_PETSc)
+  CHKERRQ(err_PETSc)
+  call VecScale(solution,timeinc,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 end subroutine FEM_mechanical_forward
 
@@ -732,24 +779,32 @@ subroutine FEM_mechanical_updateCoords()
     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)
+  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)
+  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)
+    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)
+  call VecRestoreArrayF90(x_local,nodeCoords_linear,err_PETSc)
+  CHKERRQ(err_PETSc)
 
   ! write ip displacements
   call DMPlexGetHeightStratum(dm_local,0_pPETSCINT,cellStart,cellEnd,err_PETSc)
@@ -779,7 +834,8 @@ subroutine FEM_mechanical_updateCoords()
     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)
+  call DMRestoreLocalVector(dm_local,x_local,err_PETSc)
+  CHKERRQ(err_PETSc)
 
 end subroutine FEM_mechanical_updateCoords