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