Wrapping PetscDS into petsc4py

[rbeucher]

I am trying to wrap the PetscDS functionalities into Petsc4py by completing the
cython class. The aim is to be be able to reproduce finite element pb such as
the Poisson example (petsc/src/snes/examples/tutorials/ex12.c) and ultimately
a Stokes problem such as in (petsc/src/snes/examples/tutorials/ex12.c).

The main hurdle I am facing is the wrapping of the PetscDSSetResidual (and
PetscDSSetJacobian) function which set the Residual integrand f0 and f1 as
callback functions.

see PETSc doc here

The signature of those callback functions is as follow:

void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScala
      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScala
      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])

Ideally I would like those functions to be written in python in a way similar to
what has been done for the SNESSetFunction, SNESSetJacobian etc.

It would be great if a pointer to the DS could be added to the list of arguments
passed to the above function.

Another possibility would be to add a void pointer to the argument passed to the
f* function that could be used to provide some context information...

This requires some changes in PETSc itself and I am not sure how we should
proceed. Maybe contact Matt Knepley directly?

I can't see any other way around it and I think this is also going to be a
problem for passing underworld function (whether they are C or python function).
For now julian uses the auxiliar fields on which he casts the required
information... but he still has to write the C function...

[rbeucher]

In cython I could do something like that:
Note there are still some array size issues I have not figured out...
Maybe we could get the dimension from the DS itself...

cdef void F0_Function(
    PetscInt dim,
    PetscInt Nf,
    PetscInt NfAux,
    const PetscInt uOff[],
    const PetscInt uOff_x[],
    const PetscScalar u[],
    const PetscScalar u_t[],
    const PetscScalar u_x[],
    const PetscInt aOff[],
    const PetscInt aOff_x[],
    const PetscScalar a[],
    const PetscScalar a_t[],
    const PetscScalar a_x[],
    PetscReal t,
    const PetscReal x[],
    PetscInt numConstants,
    const PetscScalar constants[],
    PetscScalar f[]):

    cdef DS Ds = ref_DS(ds) # retrieve the DS object from the C pointer
    cdef object f0 = getattr(Ds, '__f0__') # retrieve the python funtion from the DS
    cdef double[:] retv

    py_uOff = array_i(Nf, uOff)
    py_uOff_x = array_i(Nf, uOff_x)

    # Lets just make it big enough for now
    py_u = array_r(10, u)
    py_u_t = array_r(10, u_t)
    py_u_x = array_r(10, u_x)

    py_aOff = array_i(Nf, uOff)
    py_aOff_x = array_i(Nf, uOff_x)

    # Lets just make it big enough for now
    py_a = array_r(10, u)
    py_a_t = array_r(10, u_t)
    py_a_x = array_r(10, u_x)

    py_x = array_r(dim, x)
    py_constants = array_r(numConstants, constants)

    retv = f0(dim, Nf, NfAux, py_uOff, py_uOff_x, py_u, py_u_t, py_u_x, py_aOff, py_aOff_x,
       py_a, py_a_t, py_a_x, t, py_x, numConstants, py_constants)

    f = &retv[0]
    return

[julesghub]

Without the DS being passed into the F0_Function wrapper we're incapable of building the memory blocks required to construct the wrapper.
2 possible solution come to mind. We modify petsc to pass either 1) the ds or 2) a void* that could be the ds.

[lmoresi]

That sounds suspiciously like a void * context solution. Prof Louis Moresi louis.moresi@unimelb.edu.au<mailto:louis.moresi@unimelb.edu.au> (w) +61 3 8344 1217 (m) +61 4 0333 1413 (us) +1 505 349 4425 www.moresi.info<http://www.moresi.info/> www.facebook.com/underworldcode<http://www.facebook.com/underworldcode> @LouisMoresi<https://twitter.com/LouisMoresi> On 14 Jun 2019, 4:46 AM +0100, Julian Giordani <notifications@github.com>, wrote: Without the DS being passed into the F0_Function wrapper we're incapable of building the memory blocks required to construct the wrapper. 2 possible solution come to mind. We modify petsc to pass either 1) the ds or 2) a void* that could be the ds. — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#9?email_source=notifications&email_token=ADABPIYB5GS5TPMZZ7R247LP2MIAFA5CNFSM4HRQNNZKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGODXVUEWA#issuecomment-501957208>, or mute the thread<https://github.com/notifications/unsubscribe-auth/ADABPI4R6DBARVI5RDXUT6DP2MIAFANCNFSM4HRQNNZA>.

[rbeucher]

Might be good to have both 1 and 2...the DS and the void * which we could use to pass UW functions?

[rbeucher]

OK I am copying the discussion we had with Matt Knepley a while ago as it is not available anymore on bitbucket

[rbeucher]

If you look at a system like FIredrake, they compile down a
specification of the pointfunction to C and compile it, so that they
Python interacts with this at the level of full assembly. At least this
is my interpretation. Thus, what I would like to expose to Python are
functions like DMPlexComputeResidual_Internal() or
DMPlexComputeResidual_Patch_Internal(). Does that seem right to you?

[rbeucher]

This is still ongoing... I have not found a solution to this pb.