-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 17/09/15 10:17, Justin Chang wrote:
> Hi again,
>
> So I actually am not satisfied with solving those systems of
> equations as a weak form - I am getting negative concentrations
> which I suspect is caused by the Galerkin projection and the fact
> that I am not actually solving these equations. I had no such
> issues when I solved these problems with SciPy's fsolve()
> function.
>
> I want to pursue the PyOP2 option, but would still appreciate a
> little more guidance.
>
> The approach i want to take is exactly the way I did the
> variational form of this problem, but now my question is now do I
> ensure that the mapping locality (with respect to processes)
> between the solutions psi_A = Function(W.sub(0)) and psi_B =
> Function(W.sub(1)) will correspond to the pyop2 vectors/matrices?

OK, so there are two steps here.  Let's recall the system you want to
solve:

psi_A = c_A + c_C
psi_B = c_B + c_C
c_C = k1*c_A*c_B

with psi_{A,B} and k1 known and c_{A,B,C} unknown.

k1 is scalar, everything else is a vector of dimension N (although in
your previous code it looked like everything was scalar-valued?).  So
presumably you actually want to satisfy these systems "component-wise"?

psi_A[i] = c_A[i] + c_C[i]
psi_B[i] = c_B[i] + c_C[i]
c_C[i] = k1*c_A[i]*c_B[i]

for i running from 0, ..., N-1 the components of the vector.

So the first thing we need to do is write some C code that solves this
small system.  As Andrew points out (via
http://www.wolframalpha.com/input/?i=solve+x+%2B+z+%3D+a%2C+y+%2B+z+%3D+b%2C+z+%3D+c*x*y+for+x%2C+y%2C+z)
you can just write down the answer for this system I think.

I presume that these concentrations should be positive, so I'm a
little worried looking at those answers (because I can't convince
myself that they always remain positive).

void solve_system(const double psi_A, const double psi_B,
                  const double k1,
                  double *c_A, double *c_B, double *c_C) {

    // If k1 == 0, we're done:
    if (k1 == 0.0) {
        *c_C = 0.0;
        *c_A = psi_A;
        *c_B = psi_B;
        return;
    }
    // Otherwise, use formulas from above (or solve with Newton?)
    // You could, for example, just build a little SNES
    // on PETSC_COMM_SELF and provide the formfunction and solve
    // with finite differencing for the jacobian.
    ...
}

Now we need to write a PyOP2 kernel that can call this function:

void solve_kernel(const double *psi_A,
                  const double *psi_B,
                  const double *k1,
                  double *c_A,
                  double *c_B,
                  double *c_C) {

    // All the psis and cs have Nvals components (SET TO CORRECT VALUE!)
    for ( int i = 0; i < Nvals; i++ ) {
        solve_system(psi_A[i], psi_B[i], k1[0],
                     &(c_A[i]), &(c_B[i]), &(c_C[i]));
    }
}

Now in firedrake code do something like:

species_kernel = op2.Kernel("""
PASTE solve_system
PASTE solve_kernel
""", "solve_kernel")

# I assume now that everything lives on the same function space, Q.

op2.par_loop(special_kernel, Q.dof_dset.set,
             psi_A.dat(op2.READ),
             psi_B.dat(op2.READ),
             k1.dat(op2.READ),
             c_A.dat(op2.WRITE),
             c_B.dat(op2.WRITE),
             c_C.dat(op2.WRITE))

So the idea here is that you write your solver kernel that operates on
/local/ data, and then you have to tell pyop2 how to pack/unpack it
(so it iterates over the global data).  PyOP2's kernel API is
described (in some detail, hopefully understandably) here:
http://op2.github.io/PyOP2/kernels.html

Cheers,

Lawrence

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1

iQEcBAEBAgAGBQJV+pBwAAoJECOc1kQ8PEYvIJAIAPDKFAPP+NCfp3Ru1EhqzeQ+
wmjVxlDTMatW5VlBOqoH3t6PDf6LabaMo4rJETZzLP0X4lbQ0wY2k8mGUc2Phv0Y
Dmu8yOnM+uWQiZG+rp8RdpYUahO5pRry2XlKDQfCegmW4YYTMQugFnZfkxNaeiSA
Q2PyT3wkc7a7LAOt7XDZwzvGGmUi7UQw+HdCcHkaTjPG8M9M4ODrq+0jXvhU9cIs
WYvUHJeTroUD0baamx2mFD5hSFF+3MlK+1PVl+Jq6aCtr9glNmhG8np6pdj15mYm
nz+oa7O+8jno12TgxgneaZqXmA0Til7IUeiizrQOMvmZtHzhF+TQV/+0uIwisXs=
=sNXV
-----END PGP SIGNATURE-----

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake