[petsc-users] Having trouble with basic Fortran-PETSc interoperability

[Michael Whitten via petsc-users]

Hi PETSc mailing list users,

I have managed to install PETSc and run some PETSc examples and little test 
codes of my own in Fortran. I am now trying to make PETSc work with my existing 
Fortran code. I have tried to build little test examples of the functionality 
that I can then incorporate into my larger code base. However, I am having 
trouble just passing vectors back and forth between PETSc and Fortran.

I have attached a minimum semi-working example that can be compiled with the 
standard 'Makefile.user'. It throws an error when I try to copy the PETSc 
vector back to a Fortran vector using VecGetValues(). I get that it can only 
access values of the array on the local process but how do I fix this? Is this 
even the right approach?

In the final implementation I want to be able to assemble my matrix and vector, 
convert them to PETSc data structures, use PETSc to solve, and then convert the 
solution vector back to Fortran and return. I want to be able to do this with 
both the linear and nonlinear solvers. It seems like this is what PETSc is, in 
part, built to do. Is this a reasonable expectation to achieve? Is this a 
reasonable use case for PETSc?

Thanks in advance for any help you can offer.

best,
Michael

program main

#include <petsc/finclude/petscksp.h>
    use petscvec
    implicit none

    PetscErrorCode ierr
    PetscInt i
    Vec v1p

    integer :: n
    real(8), allocatable :: v1(:), v2(:)

    n = 10

    allocate(v1(n),v2(n),source=0.0d0)

    v1 = 1.0d0

    PetscCallA(PetscInitialize(ierr))

    !Create PETSc vector
    PetscCallA(VecCreate(PETSC_COMM_WORLD,v1p,ierr))
    PetscCallA(VecSetSizes(v1p,PETSC_DECIDE,n,ierr))
    PetscCallA(VecSetFromOptions(v1p,ierr))

    !Pass Fortran Vector to PETSc
    PetscCallA(VecSetValues(v1p,n,[(i,i=0,n-1)],v1,INSERT_VALUES,ierr))

    !Send PETSc vector to each processor
    PetscCallA(VecAssemblyBegin(v1p,ierr))
    PetscCallA(VecAssemblyEnd(v1p,ierr))

    !Print PETSc vector
    PetscCallA(VecView(v1p,PETSC_VIEWER_STDOUT_WORLD,ierr))

    !***ERROR HERE***
    !Pass PETSc vector to Fortran
    PetscCallA(VecGetValues(v1p,n,[(i,i=1,n)],v2,ierr))

    !Print Fortran vector
    print *, "v2: ", v2

    !PETSc clean up
    PetscCallA(VecDestroy(v1p,ierr))
    PetscCallA(PetscFinalize(ierr))


end program main