[petsc-dev] PetscSF in Fortran

[Matthew Knepley]

On Fri, Oct 13, 2017 at 12:25 AM, Adrian Croucher <a.croucher at auckland.ac.nz
> wrote:

> hi Jed
>
> I have had a go at writing a Fortran binding for the PetscSFGetGraph()
> function, along the lines of what you suggested:
>
> On 21/09/17 15:15, Jed Brown wrote:
> Unfortunately PetscSFSetGraph/PetscSFGetGraph need custom bindings for
> Fortran and they have not been written.  Shouldn't be difficult, but
> will take a little work/testing.  I would probably make the array of
> PetscSFNode be a PetscInt array of length 2n (or of shape (2,n)) --
> that's an equivalent memory representation to what we're using now.
>
> From looking at various other existing custom Fortran bindings that return
> arrays, I tried the following:
>
> PETSC_EXTERN void PETSC_STDCALL petscsfgetgraph_(PetscSF *sf, PetscInt
> *nroots, PetscInt *nleaves, F90Array1d *lptr, F90Array2d *rptr, int *ierr
> PETSC_F90_2PTR_PROTO(lptrd) PETSC_F90_2PTR_PROTO(rptrd))
> {
>   const PetscInt *ilocal;
>   const PetscSFNode *iremote;
>
>   *ierr = PetscSFGetGraph(*sf, nroots, nleaves, &ilocal, &iremote); if
> (*ierr) return;
>
>   *ierr = F90Array1dCreate((void*) ilocal, PETSC_INT, 1, *nleaves, lptr
> PETSC_F90_2PTR_PARAM(lptrd));
>   *ierr = F90Array2dCreate((void*) iremote, PETSC_INT, 1, 2, 1, *nleaves,
> rptr PETSC_F90_2PTR_PARAM(rptrd));
> }
>
> I put this in a new file zsff90.c, in a new directory
> /vec/is/sf/interface/f90-custom/, together with an appropriate makefile.
>
> However it crashes with a segmentation violation in the call to
> F90Array1dCreate(), when that goes into the Fortran routine
> F90Array1dCreateInt() and tries to assign the pointer.
>
> In my test calling program I declare the arrays for ilocal and iremote as
> PetscInt, pointer :: ilocal(:), iremote(:,:) and pass them in.
>
> Any clues? I have tried various variations on the above but I suspect I've
> either made some basic mistake or I just don't understand how this should
> be done.
>

If you want the wrapper to take in F90 pointer arguments, then you have to
declare it or you get an SEGV. These
get autogenerated in include/petsc/finclude/ftn-auto/*.h90 when you run
'make allfortranstubs'. Did this happen
for your function?

  Thanks,

     Matt


>
> - Adrian
>
> On 28/09/17 15:34, Adrian Croucher wrote:
>
> hi
> On 28/09/17 04:18, Matthew Knepley wrote:
>
>
> Okay, I think this should be easy to solve.
>
> First a little bit about SF. There are two parts to the specification. You
> have the communication part, which maps
> a certain location p on this process to another location q on another
> process. This might not change for you. The
> second part just tells it how big the data array is (numRoots), which is
> the thing that the locations in the communication
> part index into.
>
> So my question is, where did you put the numbers for the points you added?
> Are they after all points, or only after the old cells?
> I think you can easily create the new SF using the info from SFGetGraph().
>
>
> I've had a closer look at how the SF stuff works (partly by viewing the
> point SF from the original DM) and I think I understand it now. There were
> two things I hadn't realised:
>
> 1) It looks like all DM points are always considered potential roots for
> leaves on another process, which is why it complained with an error message
> when the number of roots was less than pEnd . I don't think this really
> makes so much sense in the dual-porosity mesh case I'm working on, because
> the new points I'm adding are all 'inside' the original cells and have no
> possible connection with any other points. So they can't be roots for any
> off-process leaf points. But I guess it won't hurt to tell it that they can
> (by increasing the number of roots passed to PetscSFSetGraph so it's equal
> to the new pEnd).
>
> 2) Because I'm doing finite volume I had only thought about ghost cells,
> but the SF needs to include leaf points in other height strata as well
> (vertices, edges and faces). My new points in each stratum have been added
> after the partition ghost points, so the leaf cells won't have changed
> their point indices. However the leaf points in other strata will have been
> shifted (because of points added into preceding strata).
>
>
> So I think I will need to use PetscSFSetGraph() after all, so I can
> increase the number of roots, update the leaf point indices, and also
> update the remote root index for each leaf (presumably I can use the
> original SF and PetscSFBcastBegin/ PetscSFBcastEnd to do that).
>
> If you agree, then I will need working Fortran interfaces to the
> PetscSFGetGraph/ PetscSFSetGraph functions, which are missing at present.
> Are you able to add those easily?
>
> Thanks,
> Adrian
>
>
> --
> Dr Adrian Croucher
> Senior Research Fellow
> Department of Engineering Science
> University of Auckland, New Zealand
> email: a.croucher at auckland.ac.nz
> tel: +64 (0)9 923 4611 <+64%209-923%204611>
>
>


-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.caam.rice.edu/~mk51/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20171013/8dd6b24a/attachment.html>