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}
…
===================

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?

Thanks,
Justin

On Sep 18, 2015, at 12:06 PM, Justin Chang <jychang48@gmail.com> wrote:

Ah I understand now. Got it to work now, thank you very much!

Justin

On Fri, Sep 18, 2015 at 3:41 AM, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:

> On 18 Sep 2015, at 10:29, Justin Chang <jychang48@gmail.com> wrote:
>
> Hi Lawrence,
>
> So I attempted what you had suggested but I am getting these nasty errors,
> which I do not understand at all:

...

> ===========
>
> When I look at the error log, it says "error: use of undeclared identifier
> 'Nvals'". I attached the working code. How do I declare that value?

Ah, I see.  I intended you to replace that with the number of components your reaction system has!  Since these appear to be scalars, you could replace Nvals with 1.  Or else, if you want to be a little more generic do something like this:

species_kernel = op2.Kernel("""

  // Solve reactions at each nodee
  void solve_system(const double psi_A, const double psi_B, const double k1,
    double *c_A, double *c_B, double *c_C) {
    *c_A = (sqrt(pow(psi_A*k1 - psi_B*k1 - 1,2) + 4*psi_A*k1) + psi_A*k1 - psi_B*k1 - 1)/(2*k1);
    *c_B = (sqrt(pow(psi_A*k1 - psi_B*k1 - 1,2) + 4*psi_A*k1) - psi_A*k1 + psi_B*k1     - 1)/(2*k1);
    *c_C = (sqrt(pow(psi_A*k1 - psi_B*k1 - 1,2) + 4*psi_A*k1) + psi_A*k1 + psi_B*k1     + 1)/(2*k1);
  }

  // Iterate through the node
  void solve_kernel(const double *psi_A, const double *psi_B, const double *k1,
    double *c_A, double *c_B, double *c_C) {
    const int Nvals = %(nvals)d;
    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]));
    }
  }

""" % {'nvals': Q.dof_dset.cdim}, "solve_kernel")


So all I've done here is declare:

const int Nvals = %(nvals)d;

in the kernel, and use the format spec to replace the assignment by whatever the correct value is (Q.dof_dset.cdim).  So the code the compiler sees will now be:

const int Nvals = 1;
for (int i = 0; i < Nvals; i++) {
   ...;
}

Rather than previously where Nvals was never declared.

Hope this works!

Lawrence

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