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: [email protected]
tel: +64 (0)9 923 4611
