Re: [firedrake] Solving this system of nonlinear equations
If the two roots are always distinct, there is also Weierstrass iteration. cheers --cjc On 3 September 2015 at 12:42, McRae, Andrew <a.mcrae12@imperial.ac.uk> wrote:
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).
-- http://www.imperial.ac.uk/people/colin.cotter www.cambridge.org/9781107663916
Hi everyone, Thanks for the input. So I was leaning more towards employing the newton method for solving this. And yes, all the vectors live in the same function space. If the system of equations were solved at each node i, the residual algebraic vector of size 3x1 would be: r = [c_A(i)+c_C(i)-psi_A(i); c_B(i)+c_C(i)-psi_B(i); k1*c_A(i)*c_B(i)-c_C(i)] and the Jacobian 3x3 matrix would be: J = [1, 0, 1; 0, 1, 1; k1*c_B(i), k1*c_A(i), -1] Ideally I would like to have these expressed as a PETSc matrix and vector because I would like to employ PETSc’s variational inequality solver to enforce lower/upper bounds for c_A/B/C. As for the pyop2 option:
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?
I am guessing void solve_system(…) is a separate C file called solve_system.c but how exactly do I get the components of each variable at the node? So here’s how I did this with scipy: ================================ # Advection-diffusion firedrake code # Compute psi_A, psi_B, and k1 ... … from scipy.optimize import solve import math # Solve system of nonlinear reactions def equations(p,psi_A,psi_B,k_1): c_A,c_B,c_C = p return (c_A+c_C-psi_A, c_B+c_C-psi_B, c_C-k_1*c_A*c_B) # define vectors for primary species c_A_vec = np.zeros(len(psi_A.vector())) c_B_vec = np.zeros(len(psi_A.vector())) # define vector for secondary species c_C_vec = np.zeros(len(psi_A.vector())) # Iterate through each node for i in range(len(psi_A.vector())): c_A_i,c_B_i,c_C_i = fsolve(equations,(0.0,0.0,0.0),args=(psi_A.vector()[i],psi_B.vector()[i],k1) c_A_vec[i] = c_A_i c_B_vec[i] = c_B_i c_C_vec[i] = c_C_i # Output c_A = Function(Q) c_B = Function(Q) c_C = Function(Q) c_A.vector()[:] = c_A_vec c_B.vector()[:] = c_B_vec c_C.vector()[:] = c_C_vec File(“Solution_c_A.pvd”) << c_A File(“Solution_c_B.pvd”) << c_B File(“Solution_c_C.pvd”) << c_C =================================== The PyOP2 equivalent of the above would be a good starting point I suppose. Thanks all, Justin
participants (2)
-
Colin Cotter
-
Justin Chang