<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
</head>
<body text="#000000" bgcolor="#FFFFFF">
<p>hi Jed<br>
</p>
<p>I have had a go at writing a Fortran binding for the
PetscSFGetGraph() function, along the lines of what you suggested:</p>
<p>
<blockquote type="cite">On 21/09/17 15:15, Jed Brown wrote:
<br>
Unfortunately PetscSFSetGraph/PetscSFGetGraph need custom
bindings for
<br>
Fortran and they have not been written. Shouldn't be difficult,
but
<br>
will take a little work/testing. I would probably make the
array of
<br>
PetscSFNode be a PetscInt array of length 2n (or of shape (2,n))
--
<br>
that's an equivalent memory representation to what we're using
now.
</blockquote>
</p>
<p>From looking at various other existing custom Fortran bindings
that return arrays, I tried the following:<br>
</p>
<br>
<tt>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))</tt><tt><br>
</tt><tt>{</tt><tt><br>
</tt><tt> const PetscInt *ilocal;</tt><tt><br>
</tt><tt> const PetscSFNode *iremote;</tt><tt><br>
</tt><tt><br>
</tt><tt> *ierr = PetscSFGetGraph(*sf, nroots, nleaves,
&ilocal, &iremote); if (*ierr) return;</tt><tt><br>
</tt><tt><br>
</tt><tt> *ierr = F90Array1dCreate((void*) ilocal, PETSC_INT, 1,
*nleaves, lptr PETSC_F90_2PTR_PARAM(lptrd));</tt><tt><br>
</tt><tt> *ierr = F90Array2dCreate((void*) iremote, PETSC_INT, 1,
2, 1, *nleaves, rptr PETSC_F90_2PTR_PARAM(rptrd));</tt><tt><br>
</tt><tt>}</tt><tt><br>
</tt><tt><br>
</tt>I put this in a new file zsff90.c, in a new directory
/vec/is/sf/interface/f90-custom/, together with an appropriate
makefile.<br>
<br>
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.<br>
<br>
In my test calling program I declare the arrays for ilocal and
iremote as PetscInt, pointer :: ilocal(:), iremote(:,:) and pass
them in.<br>
<br>
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.<br>
<br>
- Adrian<br>
<tt><br>
</tt>
<div class="moz-cite-prefix">On 28/09/17 15:34, Adrian Croucher
wrote:<br>
</div>
<blockquote type="cite"
cite="mid:38a9273c-b45d-74c5-bb44-ad325042c41c@auckland.ac.nz">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<p>hi<br>
</p>
<div class="moz-cite-prefix">On 28/09/17 04:18, Matthew Knepley
wrote:<br>
</div>
<blockquote type="cite"
cite="mid:CAMYG4Gn_aDBc-Mwr=tynPE_WEdQ3b5ai9ign0xGKHMbN0TGC6Q@mail.gmail.com">
<meta http-equiv="Content-Type" content="text/html;
charset=utf-8">
<div dir="ltr">
<div class="gmail_extra">
<div class="gmail_quote"><br>
<div>Okay, I think this should be easy to solve.</div>
<div><br>
</div>
<div>First a little bit about SF. There are two parts to
the specification. You have the communication part,
which maps</div>
<div>a certain location p on this process to another
location q on another process. This might not change for
you. The</div>
<div>second part just tells it how big the data array is
(numRoots), which is the thing that the locations in the
communication</div>
<div>part index into.</div>
<div><br>
</div>
<div>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?</div>
<div>I think you can easily create the new SF using the
info from SFGetGraph().</div>
</div>
</div>
</div>
</blockquote>
<br>
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:<br>
<br>
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).<br>
<br>
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).<br>
<br>
<br>
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).<br>
<br>
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?<br>
<br>
Thanks,<br>
Adrian </blockquote>
<br>
<pre class="moz-signature" cols="72">--
Dr Adrian Croucher
Senior Research Fellow
Department of Engineering Science
University of Auckland, New Zealand
email: <a class="moz-txt-link-abbreviated" href="mailto:a.croucher@auckland.ac.nz">a.croucher@auckland.ac.nz</a>
tel: +64 (0)9 923 4611
</pre>
</body>
</html>