Or, if you don't want to do Newton iterations <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>... (it may well be the case the Newton iterations are faster, I don't know. But you do have two solutions, so be careful). This still assumes that everything lives in the same functionspace, of course. On 3 September 2015 at 12:29, Mitchell, Lawrence < lawrence.mitchell@imperial.ac.uk> wrote:

...

On 3 Sep 2015, at 01:20, Justin Chang <jychang48@gmail.com> wrote:

Hi all,

I need to solve the following system of chemical reaction equations:

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

where psi_A, psi_B, and k1 are known. c_A, c_B, and c_C are the unknowns I need to solve.

k1 is a scalar and psi_A, psi_B, c_A, c_B, and c_C are vectors of size N.

psi_A and psi_B have been obtained through advection-diffusion equations. I need to "post-process" these solutions such that I get c_A, c_B, and c_C. The above system of equations are on a node-to-node basis (i.e., the reactions of node i are independent of node i+1).

I assume that all the variables live in the same functionspace. So effectively you need to do a little nonlinear iteration on a group of dofs to determine c_{A,B,C}. One approach you could use is to use the PyOP2-level kernel iteration to express this, and use a direct loop over the dofs.

This would look something like:

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

// In here, you can dereference these variables to get the components // of each variable at the node, e.g.

// Set first component of c_A to sum of first two components of psi_A. c_A[0] = psi_A[0] + psi_A[1]; }

I guess you can either write the newton (or fixed point) iteration out by hand in the kernel (it shouldn't be too difficult), or else make library calls (getting the correct link lines is possible I think).

Now in your code you do:

solve_kernel = op2.Kernel(kernel_code, "solve_system")

op2.par_loop(solve_kernel, c_A.dof_dset.node_set, c_A.dat(op2.WRITE), c_B.dat(op2.WRITE), c_C.dat(op2.WRITE), k1.dat(op2.READ), psi_A.dat(op2.READ), psi_B.dat(op2.READ))

Does this make sense?

Lawrence

...

I could use scipy's fsolve to solve the above equations but I was wondering if there's an easy way to express and solve the above using the features of firedrake and/or its dependencies (e.g., petsc, petsc4py, etc).