Re: [firedrake] Solving this system of nonlinear equations
I'm not totally sure that's obvious. I agree that it's not the theoretically optimal solution, however it is well engineered which may make it a convenient solution, and if it's fast enough, then that's just fine. On Fri, 4 Sep 2015 at 11:22 Cotter, Colin J <colin.cotter@imperial.ac.uk> wrote:
PETSc shouldn't be used for this, as the problem is solved independently at each node.
--cjc
On 4 September 2015 at 10:40, Ham, David A <david.ham@imperial.ac.uk> wrote:
OK. I think I understand. In essence, assembling PETSc matrices is a core PyOP2 competence, so it *can* do this, although it'll be a little clunky because we don't have the full glory of Firedrake on top.
I think it would help if you try to read some of the PyOP2 documentation for a bit of background:
http://op2.github.io/PyOP2/concepts.html http://op2.github.io/PyOP2/kernels.html http://op2.github.io/PyOP2/linear_algebra.html http://op2.github.io/PyOP2/mixed.html
The last one is important in this case because you're actually solving a mixed system over the mixed function space (A, B, C).
Inconveniently, that means we need 3 kernels to assemble the residual and 7 kernels (one for each nonzero entry) to assemble the Jacobian.
Once you have the code for assembling, I think it should be straightforward to provide PETSc4py with the right callbacks so that its variational inequality solver works.
(Another possibility would be to bypass PyOP2 and write your own code accessing PETSc4py directly. That might actually be easier but assembly will be slower - which you may or may not care about).
I don't have time right this instant to start writing out the Kernels (I have to teach on Monday and I'm not ready yet.) I'll try to get back to this later, or maybe someone else on the list will dive in and save me ;)
Cheers,
David
On Fri, 4 Sep 2015 at 09:37 Justin Chang <jychang48@gmail.com> wrote:
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
www.cambridge.org/9781107663916
Hi everyone, Attached is a tentative working code I have for an advective-diffusive-reactive system using operator splitting. I tried to express the system of equations as a mixed finite element form (not sure if this is the right approach). When I run the attached code, i get an error saying "Nonlinear solver failed to converge after 0 nonlinear iterations. Reason: DIVERGED_LINE_SEARCH". Can you guys have a look and see if this is even a logical approach? The parts in question begin at line 133, and the solver is at line 217 (if you comment out line 217, the advection-diffusion part runs well). If this doesn't work, I will explore the pyop2 option. Thanks, Justin On Fri, Sep 4, 2015 at 4:35 AM, David Ham <David.Ham@imperial.ac.uk> wrote:
I'm not totally sure that's obvious. I agree that it's not the theoretically optimal solution, however it is well engineered which may make it a convenient solution, and if it's fast enough, then that's just fine.
On Fri, 4 Sep 2015 at 11:22 Cotter, Colin J <colin.cotter@imperial.ac.uk> wrote:
PETSc shouldn't be used for this, as the problem is solved independently at each node.
--cjc
On 4 September 2015 at 10:40, Ham, David A <david.ham@imperial.ac.uk> wrote:
OK. I think I understand. In essence, assembling PETSc matrices is a core PyOP2 competence, so it *can* do this, although it'll be a little clunky because we don't have the full glory of Firedrake on top.
I think it would help if you try to read some of the PyOP2 documentation for a bit of background:
http://op2.github.io/PyOP2/concepts.html http://op2.github.io/PyOP2/kernels.html http://op2.github.io/PyOP2/linear_algebra.html http://op2.github.io/PyOP2/mixed.html
The last one is important in this case because you're actually solving a mixed system over the mixed function space (A, B, C).
Inconveniently, that means we need 3 kernels to assemble the residual and 7 kernels (one for each nonzero entry) to assemble the Jacobian.
Once you have the code for assembling, I think it should be straightforward to provide PETSc4py with the right callbacks so that its variational inequality solver works.
(Another possibility would be to bypass PyOP2 and write your own code accessing PETSc4py directly. That might actually be easier but assembly will be slower - which you may or may not care about).
I don't have time right this instant to start writing out the Kernels (I have to teach on Monday and I'm not ready yet.) I'll try to get back to this later, or maybe someone else on the list will dive in and save me ;)
Cheers,
David
On Fri, 4 Sep 2015 at 09:37 Justin Chang <jychang48@gmail.com> wrote:
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
www.cambridge.org/9781107663916
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hi Justin,
On 11 Sep 2015, at 00:56, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
Attached is a tentative working code I have for an advective-diffusive-reactive system using operator splitting. I tried to express the system of equations as a mixed finite element form (not sure if this is the right approach).
When I run the attached code, i get an error saying "Nonlinear solver failed to converge after 0 nonlinear iterations. Reason: DIVERGED_LINE_SEARCH". Can you guys have a look and see if this is even a logical approach? The parts in question begin at line 133, and the solver is at line 217 (if you comment out line 217, the advection-diffusion part runs well).
If this doesn't work, I will explore the pyop2 option.
You've been bitten by the most subtle part of the fenics interface! You're not the first. https://github.com/firedrakeproject/firedrake/issues/407 and https://bitbucket.org/fenics-project/dolfin/issues/194/split-u-versus-usplit Specifically, there's an unfortunate difference in the behaviour of the results of: foo = Function(MixedFunctionSpace) a, b = foo.split() as opposed to: a, b = split(foo) The former is useable in outputs (and for pointwise assignments), however, they don't function as fully-fledged "coefficients" in forms. The latter is /not/ useable in outputs (and pointless assignments). But does function correctly in forms. I was initially confused as to why you'd hand-coded the Jacobian for the species reaction at all, so I tried taking it out (in case it was wrong), and got this error: /Users/lmitche1/Documents/work/fenics/ufl/ufl/log.pyc in error(self, *message) 149 "Write error message and raise an exception." 150 self._log.error(*message) --> 151 raise self._exception_type(self._format_raw(*message)) 152 153 def begin(self, *message): UFLException: Expecting all integrals in a form to share geometric dimension, got (). Which you probably saw as well. This occurs because you've used the first split form above, and so UFL doesn't know about the split coefficients in the form and things it can't differentiate the residual with respect to the full state c. If I go ahead and change the version of split I use, everything works nicely (last least I don't get convergence errors). W_R = Q*Q*Q dc_A,dc_B,dc_C = TrialFunctions(W_R) q_A,q_B,q_C = TestFunctions(W_R) c = Function(W_R) # Use UFL's "split" here so that the residual knows about where # the components come from c_A,c_B,c_C = split(c) # Residual and Jacobian form F = (q_A*(c_A + c_C - u0_A) + q_B*(c_B + c_C - u0_B) + q_C*(c_A*c_B - c_C))*dx problem_R = NonlinearVariationalProblem(F,c) solver_R = NonlinearVariationalSolver(problem_R,solver_parameters=R_parameters,options_prefix="R_") # Now go back to the firedrake "split" so that we can do pointless # assignments later on. c_A,c_B,c_C = c.split() Cheers, Lawrence
Hi all, Hmm strangely that did the trick. Thanks! As for the "accuracy" of these concentrations going into the reactions, that's something I will need to test thoroughly since I am doing operator splitting of A and D and not solving an actual AD equation. But it is known that anisotropy and advection violates the discrete maximum principles (i.e., negative concentrations) so one of our ways to address this issue is through optimization. Justin On Fri, Sep 11, 2015 at 3:55 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,
On 11 Sep 2015, at 00:56, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
Attached is a tentative working code I have for an advective-diffusive-reactive system using operator splitting. I tried to express the system of equations as a mixed finite element form (not sure if this is the right approach).
When I run the attached code, i get an error saying "Nonlinear solver failed to converge after 0 nonlinear iterations. Reason: DIVERGED_LINE_SEARCH". Can you guys have a look and see if this is even a logical approach? The parts in question begin at line 133, and the solver is at line 217 (if you comment out line 217, the advection-diffusion part runs well).
If this doesn't work, I will explore the pyop2 option.
You've been bitten by the most subtle part of the fenics interface! You're not the first. https://github.com/firedrakeproject/firedrake/issues/407 and https://bitbucket.org/fenics-project/dolfin/issues/194/split-u-versus-usplit
Specifically, there's an unfortunate difference in the behaviour of the results of:
foo = Function(MixedFunctionSpace)
a, b = foo.split()
as opposed to:
a, b = split(foo)
The former is useable in outputs (and for pointwise assignments), however, they don't function as fully-fledged "coefficients" in forms.
The latter is /not/ useable in outputs (and pointless assignments). But does function correctly in forms.
I was initially confused as to why you'd hand-coded the Jacobian for the species reaction at all, so I tried taking it out (in case it was wrong), and got this error:
/Users/lmitche1/Documents/work/fenics/ufl/ufl/log.pyc in error(self, *message) 149 "Write error message and raise an exception." 150 self._log.error(*message) --> 151 raise self._exception_type(self._format_raw(*message)) 152 153 def begin(self, *message):
UFLException: Expecting all integrals in a form to share geometric dimension, got ().
Which you probably saw as well.
This occurs because you've used the first split form above, and so UFL doesn't know about the split coefficients in the form and things it can't differentiate the residual with respect to the full state c.
If I go ahead and change the version of split I use, everything works nicely (last least I don't get convergence errors).
W_R = Q*Q*Q dc_A,dc_B,dc_C = TrialFunctions(W_R) q_A,q_B,q_C = TestFunctions(W_R) c = Function(W_R) # Use UFL's "split" here so that the residual knows about where # the components come from c_A,c_B,c_C = split(c)
# Residual and Jacobian form F = (q_A*(c_A + c_C - u0_A) + q_B*(c_B + c_C - u0_B) + q_C*(c_A*c_B - c_C))*dx problem_R = NonlinearVariationalProblem(F,c) solver_R = NonlinearVariationalSolver(problem_R,solver_parameters=R_parameters,options_prefix="R_") # Now go back to the firedrake "split" so that we can do pointless # assignments later on. c_A,c_B,c_C = c.split()
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
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? Thanks, Justin On Fri, Sep 11, 2015 at 5:52 AM, Justin Chang <jychang48@gmail.com> wrote:
Hi all,
Hmm strangely that did the trick. Thanks!
As for the "accuracy" of these concentrations going into the reactions, that's something I will need to test thoroughly since I am doing operator splitting of A and D and not solving an actual AD equation. But it is known that anisotropy and advection violates the discrete maximum principles (i.e., negative concentrations) so one of our ways to address this issue is through optimization.
Justin
On Fri, Sep 11, 2015 at 3:55 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,
On 11 Sep 2015, at 00:56, Justin Chang <jychang48@gmail.com> wrote:
Hi everyone,
Attached is a tentative working code I have for an advective-diffusive-reactive system using operator splitting. I tried to express the system of equations as a mixed finite element form (not sure if this is the right approach).
When I run the attached code, i get an error saying "Nonlinear solver failed to converge after 0 nonlinear iterations. Reason: DIVERGED_LINE_SEARCH". Can you guys have a look and see if this is even a logical approach? The parts in question begin at line 133, and the solver is at line 217 (if you comment out line 217, the advection-diffusion part runs well).
If this doesn't work, I will explore the pyop2 option.
You've been bitten by the most subtle part of the fenics interface! You're not the first. https://github.com/firedrakeproject/firedrake/issues/407 and https://bitbucket.org/fenics-project/dolfin/issues/194/split-u-versus-usplit
Specifically, there's an unfortunate difference in the behaviour of the results of:
foo = Function(MixedFunctionSpace)
a, b = foo.split()
as opposed to:
a, b = split(foo)
The former is useable in outputs (and for pointwise assignments), however, they don't function as fully-fledged "coefficients" in forms.
The latter is /not/ useable in outputs (and pointless assignments). But does function correctly in forms.
I was initially confused as to why you'd hand-coded the Jacobian for the species reaction at all, so I tried taking it out (in case it was wrong), and got this error:
/Users/lmitche1/Documents/work/fenics/ufl/ufl/log.pyc in error(self, *message) 149 "Write error message and raise an exception." 150 self._log.error(*message) --> 151 raise self._exception_type(self._format_raw(*message)) 152 153 def begin(self, *message):
UFLException: Expecting all integrals in a form to share geometric dimension, got ().
Which you probably saw as well.
This occurs because you've used the first split form above, and so UFL doesn't know about the split coefficients in the form and things it can't differentiate the residual with respect to the full state c.
If I go ahead and change the version of split I use, everything works nicely (last least I don't get convergence errors).
W_R = Q*Q*Q dc_A,dc_B,dc_C = TrialFunctions(W_R) q_A,q_B,q_C = TestFunctions(W_R) c = Function(W_R) # Use UFL's "split" here so that the residual knows about where # the components come from c_A,c_B,c_C = split(c)
# Residual and Jacobian form F = (q_A*(c_A + c_C - u0_A) + q_B*(c_B + c_C - u0_B) + q_C*(c_A*c_B - c_C))*dx problem_R = NonlinearVariationalProblem(F,c) solver_R = NonlinearVariationalSolver(problem_R,solver_parameters=R_parameters,options_prefix="R_") # Now go back to the firedrake "split" so that we can do pointless # assignments later on. c_A,c_B,c_C = c.split()
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-----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-----
Hi Lawrence, So I attempted what you had suggested but I am getting these nasty errors, which I do not understand at all: ==================== Traceback (most recent call last): File "2D_plume_ADR_ex1.py", line 262, in <module> outFile_A << c_A File "/Users/justin/Software/firedrake/src/firedrake/firedrake/io.py", line 120, in __lshift__ self._file << data File "/Users/justin/Software/firedrake/src/firedrake/firedrake/io.py", line 534, in __lshift__ self._update_PVD(data) File "/Users/justin/Software/firedrake/src/firedrake/firedrake/io.py", line 552, in _update_PVD new_vtk << function File "/Users/justin/Software/firedrake/src/firedrake/firedrake/io.py", line 342, in __lshift__ data = output.dat.data_ro_with_halos.flatten() File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/base.py", line 1818, in data_ro_with_halos self.data_ro # force evaluation File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/base.py", line 1797, in data_ro _trace.evaluate(set([self]), set()) File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/base.py", line 169, in evaluate comp._run() File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/base.py", line 4031, in _run return self.compute() File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/base.py", line 4072, in compute fun = self._jitmodule File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/utils.py", line 64, in __get__ obj.__dict__[self.__name__] = result = self.fget(obj) File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/sequential.py", line 158, in _jitmodule direct=self.is_direct, iterate=self.iteration_region) File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line 203, in __new__ obj = make_obj() File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/caching.py", line 193, in make_obj obj.__init__(*args, **kwargs) File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/host.py", line 706, in __init__ self.compile() File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/host.py", line 804, in compile compiler=compiler.get('name')) File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/compilation.py", line 269, in load dll = compiler.get_so(src, extension) File "/Users/justin/Software/firedrake/src/PyOP2/pyop2/compilation.py", line 138, in get_so Original error: %s""" % (cc, logfile, errfile, e)) pyop2.exceptions.CompilationError: Command "['mpicc', '-std=c99', '-fPIC', '-Wall', '-framework', 'Accelerate', '-O3', '-I/Users/justin/Software/firedrake/lib/python2.7/site-packages/petsc/include', '-I/Users/justin/Software/firedrake/src/PyOP2/pyop2', '-msse', '-o', '/var/folders/jd/qh5zn6jn5kz_byz9gxz5kl2m0000gn/T/pyop2-cache-uid501/4e604d30ba50b457b02d97bff28da26e.so.tmp', '/var/folders/jd/qh5zn6jn5kz_byz9gxz5kl2m0000gn/T/pyop2-cache-uid501/4e604d30ba50b457b02d97bff28da26e.c', '-dynamiclib', '-L/Users/justin/Software/firedrake/lib/python2.7/site-packages/petsc/lib', '-Wl,-rpath,/Users/justin/Software/firedrake/lib/python2.7/site-packages/petsc/lib', '-lpetsc', '-lm']" returned with error. Unable to compile code Compile log in /var/folders/jd/qh5zn6jn5kz_byz9gxz5kl2m0000gn/T/pyop2-cache-uid501/4e604d30ba50b457b02d97bff28da26e.log Compile errors in /var/folders/jd/qh5zn6jn5kz_byz9gxz5kl2m0000gn/T/pyop2-cache-uid501/4e604d30ba50b457b02d97bff28da26e.err Original error: status 1 invoking 'mpicc -std=c99 -fPIC -Wall -framework Accelerate -O3 -I/Users/justin/Software/firedrake/lib/python2.7/site-packages/petsc/include -I/Users/justin/Software/firedrake/src/PyOP2/pyop2 -msse -o /var/folders/jd/qh5zn6jn5kz_byz9gxz5kl2m0000gn/T/pyop2-cache-uid501/4e604d30ba50b457b02d97bff28da26e.so.tmp /var/folders/jd/qh5zn6jn5kz_byz9gxz5kl2m0000gn/T/pyop2-cache-uid501/4e604d30ba50b457b02d97bff28da26e.c -dynamiclib -L/Users/justin/Software/firedrake/lib/python2.7/site-packages/petsc/lib -Wl,-rpath,/Users/justin/Software/firedrake/lib/python2.7/site-packages/petsc/lib -lpetsc -lm' =========== 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? Thanks, Justin On Thu, Sep 17, 2015 at 4:05 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
-----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-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
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
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
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 <mailto:lawrence.mitchell@imperial.ac.uk>> wrote:
On 18 Sep 2015, at 10:29, Justin Chang <jychang48@gmail.com <mailto: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 <mailto:firedrake@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/firedrake <https://mailman.ic.ac.uk/mailman/listinfo/firedrake>
-----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-----
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? 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>
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
participants (3)
-
David Ham
-
Justin Chang
-
Lawrence Mitchell