-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 On 17/09/15 10:17, Justin Chang wrote:
Hi again,
So I actually am not satisfied with solving those systems of equations as a weak form - I am getting negative concentrations which I suspect is caused by the Galerkin projection and the fact that I am not actually solving these equations. I had no such issues when I solved these problems with SciPy's fsolve() function.
I want to pursue the PyOP2 option, but would still appreciate a little more guidance.
The approach i want to take is exactly the way I did the variational form of this problem, but now my question is now do I ensure that the mapping locality (with respect to processes) between the solutions psi_A = Function(W.sub(0)) and psi_B = Function(W.sub(1)) will correspond to the pyop2 vectors/matrices?
OK, so there are two steps here. Let's recall the system you want to solve: psi_A = c_A + c_C psi_B = c_B + c_C c_C = k1*c_A*c_B with psi_{A,B} and k1 known and c_{A,B,C} unknown. k1 is scalar, everything else is a vector of dimension N (although in your previous code it looked like everything was scalar-valued?). So presumably you actually want to satisfy these systems "component-wise"? psi_A[i] = c_A[i] + c_C[i] psi_B[i] = c_B[i] + c_C[i] c_C[i] = k1*c_A[i]*c_B[i] for i running from 0, ..., N-1 the components of the vector. So the first thing we need to do is write some C code that solves this small system. As Andrew points out (via http://www.wolframalpha.com/input/?i=solve+x+%2B+z+%3D+a%2C+y+%2B+z+%3D+b%2C...) you can just write down the answer for this system I think. I presume that these concentrations should be positive, so I'm a little worried looking at those answers (because I can't convince myself that they always remain positive). void solve_system(const double psi_A, const double psi_B, const double k1, double *c_A, double *c_B, double *c_C) { // If k1 == 0, we're done: if (k1 == 0.0) { *c_C = 0.0; *c_A = psi_A; *c_B = psi_B; return; } // Otherwise, use formulas from above (or solve with Newton?) // You could, for example, just build a little SNES // on PETSC_COMM_SELF and provide the formfunction and solve // with finite differencing for the jacobian. ... } Now we need to write a PyOP2 kernel that can call this function: void solve_kernel(const double *psi_A, const double *psi_B, const double *k1, double *c_A, double *c_B, double *c_C) { // All the psis and cs have Nvals components (SET TO CORRECT VALUE!) 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])); } } Now in firedrake code do something like: species_kernel = op2.Kernel(""" PASTE solve_system PASTE solve_kernel """, "solve_kernel") # I assume now that everything lives on the same function space, Q. op2.par_loop(special_kernel, Q.dof_dset.set, 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)) So the idea here is that you write your solver kernel that operates on /local/ data, and then you have to tell pyop2 how to pack/unpack it (so it iterates over the global data). PyOP2's kernel API is described (in some detail, hopefully understandably) here: http://op2.github.io/PyOP2/kernels.html Cheers, Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iQEcBAEBAgAGBQJV+pBwAAoJECOc1kQ8PEYvIJAIAPDKFAPP+NCfp3Ru1EhqzeQ+ wmjVxlDTMatW5VlBOqoH3t6PDf6LabaMo4rJETZzLP0X4lbQ0wY2k8mGUc2Phv0Y Dmu8yOnM+uWQiZG+rp8RdpYUahO5pRry2XlKDQfCegmW4YYTMQugFnZfkxNaeiSA Q2PyT3wkc7a7LAOt7XDZwzvGGmUi7UQw+HdCcHkaTjPG8M9M4ODrq+0jXvhU9cIs WYvUHJeTroUD0baamx2mFD5hSFF+3MlK+1PVl+Jq6aCtr9glNmhG8np6pdj15mYm nz+oa7O+8jno12TgxgneaZqXmA0Til7IUeiizrQOMvmZtHzhF+TQV/+0uIwisXs= =sNXV -----END PGP SIGNATURE-----