On Oct 1, 2015, at 3:39 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:

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

On 01/10/15 09:54, Justin Chang wrote:

Hello again :)

Now I need to do an even more complicated problem:

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

Unfortunately the online calculator gives some stupidly long and
convoluted answer, and my problems are only going to get more
complicated. Seems to me I need to implement the Newton’s method.

All five vectors live in Q = FunctionSpace(mesh,”CG”,1)
individually, whereas the aggregate psi_{A/B} lives in W = Q*Q and
c_{A/B/C} lives in W_R = Q*Q*Q. Basically I need to create a matrix
(sparsity?) that corresponds to the ordering/mapping of W_R.

Would I do something like this?

================= mds = W_R.dof_dset mmap = W_R.cell_node_map J =
op2.Mat(op2.Sparsity(mds,mmap)) F = Function(W_R)

op2.par_loop(kernel,mds.set,psi_A.dat(op2.READ),psi_B.dat(op2.READ),k1.dat(op2.READ),

c_A.dat(op2.READ), c_B.dat(op2.READ), c_C.dat(op2.READ),

F.dat(op2.WRITE,mmap), J(op2.INC,
(mmap[op2.i[0]],mmap[op2.i[0]])))

# Perform newton-rhason and update c_{A/B/C} … ===================

This is the approach to use if you want to form the global system.
However, you probably don't need to do that, since all your nodes are
decoupled (you can instead solve small local systems at each node).

I would instead do something like:

op2.par_loop(kernel, Q.dof_dset, 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))

The kernel is then.  Note how for direct loops all arguments are
passed as pointers.  You need to dereference them to read and write.

void kernel(const double *psi_A, const double *psi_B, const double *k1,
           double *c_A, double *c_B, double *c_C)
{
   // Assemble 3x3 Jacobian and 3x1 residual.
   // Do Newton-Raphson on resulting system and write results
   // to c_{A,B,C}
}

For such a small system, you're probably best inverting the matrix
explicitly.  With a very small amount of effort, we can support
kernels that use C++, rather than just C, so if you like you could use
Eigen to do some of this.

And kernel would be something like:

void kernel(const double psi_A, const double psi_B, const double
k1, const double c_A, const double c_B, const double c_C, double
**f, double j[3][3]) {

// Assemble 3x3 Jacobian j and 3x1 residual f

}

Is that the general idea here? Also I am looking through the pyop2
documentation and it seems it only does linear algebra, meaning I
would have to do my own Nonlinear/Newton iteration?

Yes, PyOP2 isn't a general purpose solver library, it's a mesh
iteration library.  There are some generic routines for pointwise
linear algebra operations, but nothing complicated.

Cheers,

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

iQEcBAEBAgAGBQJWDP9cAAoJECOc1kQ8PEYvxUAH/1Bp4iLeR13nlu4k55ciPvh4
SgJhiXg04C00q7lYqLSl5IZ6jdJ5etZ39ETDi6o8tT0UUQXd4I/wz5ovuaRmopVA
AV27/4lfSnNKocwmKeTltQIC60UwafMDGDmaur8/5dYrdwzTl11Rj2ygXkp2lCR/
z2jHmT0dSbOjhbGulObe/yFvU2upHspj4Ycc0zJX20/MpcSOfM5mFntubmnoab3B
LunILgifefFh2+OhLhy+aytIkcP7haQrlVAHv+61IwAN6InRUc/SjAC5qZAjNXfV
EtuYIDtT1EnUN0/haWpFxj8bf3/1AJd8QU0dQAzQxra9AQcx6/dcZQzJ4mUERjo=
=PrS+
-----END PGP SIGNATURE-----

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