On 2 Oct 2015, at 09:43, Justin Chang <jychang48@gmail.com> wrote:
Hi Lawrence,
Yes I would like to use Eigen for this. How would I go about this then? Because I installed Eigen onto my computer but I am guessing I would need to put the #include <Eigen> headers into the kernel code?
Yes, you can just write: kernel_code = """ #include <Eigen/Dense> using namespace Eigen; void kernel(...) { ...; } """ At the moment, we don't support C++ kernels, but I have a PyOP2 branch (cpp-kernels) which does this. You write: kernel = op2.Kernel(kernel_code, kernel_name, cpp=True) and then hopefully things works. Lawrence
Thanks, Justin
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 <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake<https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake