[petsc-users] Migrating to parallel PETSc

[Bojan Niceno]

Dear PETSc users,


I have an in-house computational fluid dynamics (CFD) solver, written in
Fortran 2008, parallelized with MPI with its own home-grown suite of linear
solvers.  The code is unstructured, performs domain decomposition with
METIS and all communication buffers, I mean connectivity between
processors, has been properly worked out.

A couple of weeks back, I decided to try out the PETSc suite of solvers.
After some initial setbacks, I managed to compile my code with PETSc and
have the sequential version working fine :-)

I have essentially using the following PETSc routines to get the code
solving linear systems with PETSc:

 I set up the working space as follows:

  call PetscInitialize(PETSC_NULL_CHARACTER, Pet % petsc_err)
  call MatCreateSeqAij(PETSC_COMM_SELF, ...
  call MatSeqAijSetColumnIndices(....
  call VecCreateSeq(PETSC_COMM_SELF,  ...  ! Create Vector x
  call VecCreateSeq(PETSC_COMM_SELF,  ...  ! Create Vector b
  call KSPCreate(PETSC_COMM_SELF,  ...

Then in order to solve a system, I do:

  call MatSetValue(Pet % petsc_A,   ! Inside a loop through matrix entries
                   i-PETSC_ONE,
                   k-PETSC_ONE, ...
  call MatAssemblyBegin(Pet % petsc_A, MAT_FINAL_ASSEMBLY, Pet % petsc_err)
  call MatAssemblyEnd  (Pet % petsc_A, MAT_FINAL_ASSEMBLY, Pet % petsc_err)

  call VecSetValue(Pet % petsc_x, ... ! Fill up x
  call VecSetValue(Pet % petsc_b, ... ! Fill up b

  call KSPSetType(Pet % petsc_ksp ... ! Set solver
  call KSPGetPC(Pet % petsc_ksp, ...  ! get preconditioner context
  call PCSetType(Pet % petsc_pc, ...  ! Set preconditioner

  call KSPSetFromOptions(Pet % petsc_ksp, Pet % petsc_err)

  call KSPSetUp         (Pet % petsc_ksp, Pet % petsc_err)


  ! Finally solve
  call KSPSolve(Pet % petsc_ksp, ...

Once this was up and running, I thought that in order to have the parallel
version I will merely have to replace the "Seq" versions of the above
functions, with their parallel counterparts.  I was expecting to find the
red function (MatSeqAijSetColumnIndices) for parallel runs, but it doesn't
seem to exist.  I have found non-seq versions of some other functions (
MatCreateAij, VecCreateSeq), but not something like MatAijSetColumnIndices,
which surprised me a bit, because I have this information in my code.

Is there a parallel counterpart of this function, and if there is none,
what should it be replaced with?  I understand that I will have to provide
non-zeros in buffers (o_nnz), which is not a big issue, but how to provide
information on columns for parallel version is not clear to me.  In a
nutshell, I would need a hint on which of the above functions could remain
the same in parallel, and which should be replaced and with what?

    Cheers,

    Bojan
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20220215/a594c196/attachment-0001.html>