[petsc-dev] PetscSF in Fortran

[Adrian Croucher]

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.

- 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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20171013/8017437e/attachment-0001.html>