Set Null space on Poisson matrices

Author: pbartholomew08

src/mesh.f90

@@ -13,7 +13,7 @@ module m_mesh
   type :: geo_t
     real(dp), dimension(3) :: d ! size of a cell in each direction (=edge length, distance between centers, distance between vertices)
     real(dp), dimension(3) :: L ! Global dimensions of the domain in each direction
-  end type
+  end type geo_t
 
   ! Stores parallel domain related information
   type :: parallel_t
@@ -26,7 +26,7 @@ module m_mesh
     integer, dimension(3) :: pprev ! rank ID of the next rank in each direction
   contains
     procedure :: is_root ! returns if the current rank is the root rank
-  end type
+  end type parallel_t
 
   ! The mesh class stores all the information about the global and local (due to domain decomposition) mesh
   ! It also includes getter functions to access some of its parameters

src/petsc/poisson_cg.f90

@@ -216,6 +216,7 @@ subroutine create_operator(self, n)
     class(petsc_poisson_cg_t) :: self
     integer, intent(in) :: n
 
+    type(tMatNullSpace) :: nsp
     integer :: ierr
 
     call MatCreateShell(PETSC_COMM_WORLD, n, n, PETSC_DETERMINE, &
@@ -227,13 +228,18 @@ subroutine create_operator(self, n)
     call MatAssemblyBegin(self%Amat, MAT_FINAL_ASSEMBLY, ierr)
     call MatAssemblyEnd(self%Amat, MAT_FINAL_ASSEMBLY, ierr)
 
+    call MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL_VEC, nsp, ierr)
+    call MatSetnullSpace(self%Amat, nsp, ierr)
+    call MatNullSpaceDestroy(nsp, ierr)
+
   end subroutine create_operator
 
   subroutine create_preconditioner(self, mesh, n)
     class(petsc_poisson_cg_t) :: self
     type(mesh_t), intent(in) :: mesh
     integer, intent(in) :: n
 
+    type(tMatNullSpace) :: nsp
     integer :: ierr
 
     integer, parameter :: nnb = 6 ! Number of neighbours (7-point star has 6 neighbours)
@@ -324,6 +330,10 @@ subroutine create_preconditioner(self, mesh, n)
     call MatAssemblyBegin(self%Pmat, MAT_FINAL_ASSEMBLY, ierr)
     call MatAssemblyEnd(self%Pmat, MAT_FINAL_ASSEMBLY, ierr)
 
+    call MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL_VEC, nsp, ierr)
+    call MatSetnullSpace(self%Pmat, nsp, ierr)
+    call MatNullSpaceDestroy(nsp, ierr)
+
   end subroutine create_preconditioner
 
   subroutine create_vectors(self, n)
@@ -356,7 +366,7 @@ subroutine create_solver(self)
     integer :: ierr
 
     call KSPCreate(PETSC_COMM_WORLD, self%ksp, ierr)
-    call KSPSetOperators(self%ksp, self%Amat, self%Pmat, ierr)
+    call KSPSetOperators(self%ksp, self%Pmat, self%Pmat, ierr)
     call KSPSetFromOptions(self%ksp, ierr)
 
   end subroutine create_solver