Calling PETSc within PyOP2 kernels
Hi all, I haven't experimented with this yet, but is it possible (or rather, efficient) to call PETSc routines (e.g., SNES) within a PyOP2 kernel? Because I really would like to (or need to) use PETSc's Variational Inequality solver point-wise. I have done this via petsc4py, but it seems the PyOP2 kernel is much faster in terms of accessing point wise values (which iirc is what one of you guys had told me in another thread). My only concern is if additional headers/include/library files would be needed if I wanted the C/C++ version of PETSc inside my PyOP2 kernels. Thanks, Justin
Hi Justin,
On 30 Oct 2015, at 05:41, Justin Chang <jychang48@gmail.com> wrote:
Hi all,
I haven't experimented with this yet, but is it possible (or rather, efficient) to call PETSc routines (e.g., SNES) within a PyOP2 kernel? Because I really would like to (or need to) use PETSc's Variational Inequality solver point-wise. I have done this via petsc4py, but it seems the PyOP2 kernel is much faster in terms of accessing point wise values (which iirc is what one of you guys had told me in another thread).
My only concern is if additional headers/include/library files would be needed if I wanted the C/C++ version of PETSc inside my PyOP2 kernels.
Yes, you can definitely do this. In fact, for global matrix assembly we effectively do the same thing. So petsc.h is already included and the generated library is linked against libpetsc. So for a first attempt, see if you can do the naive thing, creating the SNES in the kernel. A sketch below: static void pointFormFunction(...) { ...; } static void pointFormJacobian(...) { ...; } static void solvePointwiseSystem(double *a, double *b, ...) { SNES snes; Mat jac; Vec x, u; SNESCreate(PETSC_COMM_SELF, &snes); MatCreate(PETSC_COMM_SELF, &jac); MatSetType(...); SNESSetComputeJacobian(snes, jac, ..., pointFormJacobian); etc... } You should then be able to wrap all that up in a PyOP2 Kernel object: kernel = op2.Kernel(..., name="solvePointwiseSystem") And execute it over all the data. This may be slow because you build and destroy a SNES, Mat and some Vecs for every unknown you want to solve for. I assume that the system is the same everywhere, and so it's only the residual (and hence the /values/ in the jacobian) that changes at each point. If this is the case, you can make the SNES, Mat and Vec objects once (in petsc4py) and then pass handles into the kernel (a bit like PETSc's AppCtx arguments). This latter approach is a bit hacky, so let's get the simpler one working first. Lawrence
Lawrence, So I have implemented PETSc into the PYOP2 kernel, but I am sometimes getting strange errors. Attached is the working code plus the reactions kernel file. Run the code by typing "python 1D_analytical_ADR_ex1.py <seed> <0/1> <0/1>". The last option is for using either Newtons Method (0) or Variational Inequality (1). The reason my errors are strange is because it only happens sometimes. Sometimes I get errors like this: python(11976,0x7fff7497a000) malloc: *** error for object 0x7fbb4b5e8ef0: incorrect checksum for freed object - object was probably modified after being freed. *** set a breakpoint in malloc_error_break to debug [Justins-MBP:11976] *** Process received signal *** [Justins-MBP:11976] Signal: Abort trap: 6 (6) [Justins-MBP:11976] Signal code: (0) [Justins-MBP:11976] [ 0] 0 libsystem_platform.dylib 0x00007fff94c5452a _sigtramp + 26 [Justins-MBP:11976] [ 1] 0 ??? 0x0000000000000000 0x0 + 0 [Justins-MBP:11976] [ 2] 0 libsystem_c.dylib 0x00007fff946f937b abort + 129 [Justins-MBP:11976] [ 3] 0 libsystem_malloc.dylib 0x00007fff9502b3a6 szone_error + 626 [Justins-MBP:11976] [ 4] 0 libsystem_malloc.dylib 0x00007fff95021604 tiny_free_list_remove_ptr + 289 [Justins-MBP:11976] [ 5] 0 libsystem_malloc.dylib 0x00007fff9501f956 szone_free_definite_size + 1480 [Justins-MBP:11976] [ 6] 0 Python 0x000000010228f236 dict_dealloc + 154 [Justins-MBP:11976] [ 7] 0 Python 0x00000001022af3f2 subtype_dealloc + 481 [Justins-MBP:11976] [ 8] 0 Python 0x00000001022859c5 list_dealloc + 105 [Justins-MBP:11976] [ 9] 0 Python 0x0000000102280216 frame_dealloc + 110 [Justins-MBP:11976] [10] 0 Python 0x00000001022dd1c4 PyEval_EvalCodeEx + 1915 [Justins-MBP:11976] [11] 0 Python 0x00000001022e3640 fast_function + 117 [Justins-MBP:11976] [12] 0 Python 0x00000001022e0553 PyEval_EvalFrameEx + 13165 [Justins-MBP:11976] [13] 0 Python 0x00000001022dcfb4 PyEval_EvalCodeEx + 1387 [Justins-MBP:11976] [14] 0 Python 0x0000000102281bf5 function_call + 352 [Justins-MBP:11976] [15] 0 Python 0x0000000102263ad7 PyObject_Call + 99 [Justins-MBP:11976] [16] 0 Python 0x0000000102263c65 call_function_tail + 72 [Justins-MBP:11976] [17] 0 Python 0x0000000102263c06 PyObject_CallFunction + 196 [Justins-MBP:11976] [18] 0 Python 0x0000000102294755 _PyObject_GenericGetAttrWithDict + 219 [Justins-MBP:11976] [19] 0 Python 0x00000001022df15d PyEval_EvalFrameEx + 8055 [Justins-MBP:11976] [20] 0 Python 0x00000001022dcfb4 PyEval_EvalCodeEx + 1387 [Justins-MBP:11976] [21] 0 Python 0x0000000102281bf5 function_call + 352 [Justins-MBP:11976] [22] 0 Python 0x0000000102263ad7 PyObject_Call + 99 [Justins-MBP:11976] [23] 0 Python 0x0000000102263c65 call_function_tail + 72 [Justins-MBP:11976] [24] 0 Python 0x0000000102263c06 PyObject_CallFunction + 196 [Justins-MBP:11976] [25] 0 Python 0x0000000102294755 _PyObject_GenericGetAttrWithDict + 219 [Justins-MBP:11976] [26] 0 Python 0x00000001022df15d PyEval_EvalFrameEx + 8055 [Justins-MBP:11976] [27] 0 Python 0x00000001022dcfb4 PyEval_EvalCodeEx + 1387 [Justins-MBP:11976] [28] 0 Python 0x0000000102281bf5 function_call + 352 [Justins-MBP:11976] [29] 0 Python 0x0000000102263ad7 PyObject_Call + 99 [Justins-MBP:11976] *** End of error message *** [1] 11976 abort python 1D_analytical_ADR_ex1.py 100 0 1 Other times I get errors like this: [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger [0]PETSC ERROR: or see http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind [0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors [0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run [0]PETSC ERROR: to get more information on the crash. -------------------------------------------------------------------------- MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD with errorcode 59. NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes. You may or may not see output from other processes, depending on exactly when Open MPI kills them. -------------------------------------------------------------------------- I am confused as to why these errors occur randomly, even if I use the same input arguments. I suspect it has sometime to do with how I have implemented PETSc within PYOP2. Can you tell what is going on? Thanks, Justin
Hi Justin,
On 3 Nov 2015, at 06:59, Justin Chang <jychang48@gmail.com> wrote:
Lawrence,
So I have implemented PETSc into the PYOP2 kernel, but I am sometimes getting strange errors.
Attached is the working code plus the reactions kernel file.
Run the code by typing "python 1D_analytical_ADR_ex1.py <seed> <0/1> <0/1>". The last option is for using either Newtons Method (0) or Variational Inequality (1).
I think the following is the problem. Your kernel does: void solve_kernel(const double *psiVar, double *cVar, const double *kVar) { ... for (i = 0; i < numPsi; i++) user.psi[i] = psiVar[i]; for (i = 0; i < numKeq; i++) user.keq[i] = kVar[i]; ... } Where numPsi is 2 (and later on copies out 3 values into cVar. However, the way data is passed to the kernel is like so: void wrap_solve_kernel(int start, int end, double *arg0_0, double *arg0_1, double *arg1_0, double *arg1_1, double *arg1_2, double *arg2_0) { for ( int n = start; n < end; n++ ) { int i = n; solve_kernel(arg0_0 + (i * 1), arg1_0 + (i * 1), arg2_0); } } Which comes from the par_loop invocation: op2.par_loop(reactions_kernel,Q.dof_dset, u0.dat(op2.READ), c.dat(op2.WRITE), k1.dat(op2.READ)) So as we can see, the mixed Function u0 (with two components, each scalar-valued) is passed to the calling wrapper in two parts (arg0_0 and arg0_1). Similarly c (with three parts) is passed in three bits). However, only the first components are passed on to the kernel itself (which then reads past the end of arrays). Unfortunately, it looks like the correct code generation doesn't work for this use case. Instead, you can make things work by explicitly passing the components to the par_loop and accepting them in the kernel. Something like: void solve_kernel(const double *psi0, const double *psi1, double *c0, double *c1, double *c2, ...) { ... } And then write the par_loop as: u0_A, u0_B = u0.split() c0, c1, c2 = c.split() op2.par_loop(reactions_kernel, Q.dof_dset, u0_A(op2.READ), u0_B(op2.READ), c0(op2.WRITE), c1(op2.WRITE), c2(op2.WRITE), ...) However, this is clearly not very natural, so we should probably fix the code generation for your use case (which I think is the natural thing). In fact, if you run in a more debug-friendly mode in firedrake, where PyOP2 does additional run-time type checking, your code throws a more friendly error telling you this is unsupported behaviour. This should probably be switched back to being the default in firedrake (see https://github.com/firedrakeproject/firedrake/issues/626) Cheers, Lawrence
On 3 Nov 2015, at 11:56, Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,
On 3 Nov 2015, at 06:59, Justin Chang <jychang48@gmail.com> wrote:
Lawrence,
So I have implemented PETSc into the PYOP2 kernel, but I am sometimes getting strange errors.
Attached is the working code plus the reactions kernel file.
Run the code by typing "python 1D_analytical_ADR_ex1.py <seed> <0/1> <0/1>". The last option is for using either Newtons Method (0) or Variational Inequality (1).
I note as well that if I remove the bad par_loop code and run with PETSc in debug mode, I get the following error (running with options 5 1 1): Time = 0.010000 Traceback (most recent call last): File "1D_analytical_ADR_ex1.py", line 187, in <module> u0 = ADsolve(solver_D,u0) File "1D_analytical_ADR_ex1.py", line 171, in ADsolve solver_O.solve(u0_O,b) File "/Users/lmitche1/Documents/work/mapdes/firedrake/firedrake/linear_solver.py", line 138, in solve self.ksp.solve(rhs, solution) File "PETSc/KSP.pyx", line 363, in petsc4py.PETSc.KSP.solve (src/petsc4py.PETSc.c:152597) petsc4py.PETSc.Error: error code 77 [0] KSPSolve() line 527 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/ksp/interface/itfunc.c [0] KSPSetUp() line 332 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 984 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/pc/interface/precon.c [0] PCSetUp_GAMG() line 566 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/pc/impls/gamg/gamg.c [0] PCGAMGCoarsen_AGG() line 997 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/pc/impls/gamg/agg.c [0] smoothAggs() line 353 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/pc/impls/gamg/agg.c [0] Petsc has generated inconsistent data [0] !(matA_1 && !matA_1->compressedrow.use) ***[0]PCReset_GAMG this should not happen, cleaned up in SetUp Lawrence
Hi Lawrence, So I had a nice working code where it runs point-wise SNES in pyop2, and it ran decently fast. However, with the latest firedrake update, it forced me to make a few additional changes to the pyop2 kernel, namely adding PetcInitialize(....) and PetscFinalize(), and my code now runs really slow after I run it after several time steps on a Unit Interval of 100 elements. Attached is the code as well as the supplemental 1D analytical expressions and pyop2 reactions kernel. Run the code as "python 1D_OS_analytical_ex1.py 100 0 0". If I comment out line 250 in the 1D_OS_analytical.ex1.py file, the program performs well, w.r.t. wall-clock time that I normally get, but now it has become a serious bottleneck. Is there a way to improve what I have written in the reactions kernel? Thanks, Justin On Tue, Nov 3, 2015 at 5:01 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
On 3 Nov 2015, at 11:56, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
Hi Justin,
On 3 Nov 2015, at 06:59, Justin Chang <jychang48@gmail.com> wrote:
Lawrence,
So I have implemented PETSc into the PYOP2 kernel, but I am sometimes getting strange errors.
Attached is the working code plus the reactions kernel file.
Run the code by typing "python 1D_analytical_ADR_ex1.py <seed> <0/1> <0/1>". The last option is for using either Newtons Method (0) or Variational Inequality (1).
I note as well that if I remove the bad par_loop code and run with PETSc in debug mode, I get the following error (running with options 5 1 1):
Time = 0.010000 Traceback (most recent call last): File "1D_analytical_ADR_ex1.py", line 187, in <module> u0 = ADsolve(solver_D,u0) File "1D_analytical_ADR_ex1.py", line 171, in ADsolve solver_O.solve(u0_O,b) File "/Users/lmitche1/Documents/work/mapdes/firedrake/firedrake/linear_solver.py", line 138, in solve self.ksp.solve(rhs, solution) File "PETSc/KSP.pyx", line 363, in petsc4py.PETSc.KSP.solve (src/petsc4py.PETSc.c:152597) petsc4py.PETSc.Error: error code 77 [0] KSPSolve() line 527 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/ksp/interface/itfunc.c [0] KSPSetUp() line 332 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/ksp/interface/itfunc.c [0] PCSetUp() line 984 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/pc/interface/precon.c [0] PCSetUp_GAMG() line 566 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/pc/impls/gamg/gamg.c [0] PCGAMGCoarsen_AGG() line 997 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/pc/impls/gamg/agg.c [0] smoothAggs() line 353 in /Users/lmitche1/Documents/work/mapdes/petsc/src/ksp/pc/impls/gamg/agg.c [0] Petsc has generated inconsistent data [0] !(matA_1 && !matA_1->compressedrow.use) ***[0]PCReset_GAMG this should not happen, cleaned up in SetUp
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Hi Justin, On 12/11/15 09:22, Justin Chang wrote:
Hi Lawrence,
So I had a nice working code where it runs point-wise SNES in pyop2, and it ran decently fast. However, with the latest firedrake update, it forced me to make a few additional changes to the pyop2 kernel, namely adding PetcInitialize(....) and PetscFinalize(), and my code now runs really slow after I run it after several time steps on a Unit Interval of 100 elements.
So you shouldn't have to do this, something is definitely not right here .
Attached is the code as well as the supplemental 1D analytical expressions and pyop2 reactions kernel. Run the code as "python 1D_OS_analytical_ex1.py 100 0 0".
I tried this with firedrake commit 6134b86, PyOP2 commit 64db689 (i.e. current master for both). And didn't run into this problem. Could it be that the compiled code from before the update was explicitly linked (with a baked-in rpath) against a different petsc to the one you have now? To assuage my fears, can you run firedrake/scripts/firedrake-clean to clear out the compiled kernel cache? And then run again without calling PetscInitialize/Finalize inside your kernel. Cheers, Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux) iQEcBAEBAgAGBQJWRF3NAAoJECOc1kQ8PEYvOuoIAMSymNuMusrzJOmYAFTnK+gR EPpuBzj3tCt5QIBSwGFuKgW8+IuWTmT4kY8t0gBfAx6Diw8RrQ6WSot9nKB/xdI2 EUTz9RuoJQBWbzgmvP1mxvnDfT+oulzj5arNEeFdHfxDrtaoKgwWU8PrLznegEGY e7Z9MJNH3S+DaamYeOF/gXVbNilrCHaKiKaKVsIMc2HYitBdOIwaJliOBTkFGwKC dnqLn2YCsX8CpPu9939UwxcLnG9BPfyFuEguSodMXVjZKAxLQh6DgIZ25U9F4q80 etmLHaubbgVTZfahMHH/5UxdLYr+ILzXmzikzQzkP20rrvZHPdauYhz5F4pGp4M= =0WNr -----END PGP SIGNATURE-----
Lawrence, So I did what you said, but still got the same error: PETSC ERROR: Logging has not been enabled. You might have forgotten to call PetscInitialize(). -------------------------------------------------------------------------- MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD with errorcode 56. NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes. You may or may not see output from other processes, depending on exactly when Open MPI kills them. -------------------------------------------------------------------------- Any ideas why this is happening? Justin On Thu, Nov 12, 2015 at 2:37 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
Hi Justin,
On 12/11/15 09:22, Justin Chang wrote:
Hi Lawrence,
So I had a nice working code where it runs point-wise SNES in pyop2, and it ran decently fast. However, with the latest firedrake update, it forced me to make a few additional changes to the pyop2 kernel, namely adding PetcInitialize(....) and PetscFinalize(), and my code now runs really slow after I run it after several time steps on a Unit Interval of 100 elements.
So you shouldn't have to do this, something is definitely not right here .
Attached is the code as well as the supplemental 1D analytical expressions and pyop2 reactions kernel. Run the code as "python 1D_OS_analytical_ex1.py 100 0 0".
I tried this with firedrake commit 6134b86, PyOP2 commit 64db689 (i.e. current master for both). And didn't run into this problem. Could it be that the compiled code from before the update was explicitly linked (with a baked-in rpath) against a different petsc to the one you have now? To assuage my fears, can you run firedrake/scripts/firedrake-clean to clear out the compiled kernel cache? And then run again without calling PetscInitialize/Finalize inside your kernel.
Cheers,
Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)
iQEcBAEBAgAGBQJWRF3NAAoJECOc1kQ8PEYvOuoIAMSymNuMusrzJOmYAFTnK+gR EPpuBzj3tCt5QIBSwGFuKgW8+IuWTmT4kY8t0gBfAx6Diw8RrQ6WSot9nKB/xdI2 EUTz9RuoJQBWbzgmvP1mxvnDfT+oulzj5arNEeFdHfxDrtaoKgwWU8PrLznegEGY e7Z9MJNH3S+DaamYeOF/gXVbNilrCHaKiKaKVsIMc2HYitBdOIwaJliOBTkFGwKC dnqLn2YCsX8CpPu9939UwxcLnG9BPfyFuEguSodMXVjZKAxLQh6DgIZ25U9F4q80 etmLHaubbgVTZfahMHH/5UxdLYr+ILzXmzikzQzkP20rrvZHPdauYhz5F4pGp4M= =0WNr -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Strange, I had no issue running this on my Ubuntu system, it only happened with my MacOSX (El Capitan). Perhaps I need to reinstall firedrake on my MacOSX and will holler if I get the same error again On Thu, Nov 12, 2015 at 2:46 AM, Justin Chang <jychang48@gmail.com> wrote:
Lawrence,
So I did what you said, but still got the same error:
PETSC ERROR: Logging has not been enabled.
You might have forgotten to call PetscInitialize().
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 56.
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.
--------------------------------------------------------------------------
Any ideas why this is happening?
Justin
On Thu, Nov 12, 2015 at 2:37 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1
Hi Justin,
On 12/11/15 09:22, Justin Chang wrote:
Hi Lawrence,
So I had a nice working code where it runs point-wise SNES in pyop2, and it ran decently fast. However, with the latest firedrake update, it forced me to make a few additional changes to the pyop2 kernel, namely adding PetcInitialize(....) and PetscFinalize(), and my code now runs really slow after I run it after several time steps on a Unit Interval of 100 elements.
So you shouldn't have to do this, something is definitely not right here .
Attached is the code as well as the supplemental 1D analytical expressions and pyop2 reactions kernel. Run the code as "python 1D_OS_analytical_ex1.py 100 0 0".
I tried this with firedrake commit 6134b86, PyOP2 commit 64db689 (i.e. current master for both). And didn't run into this problem. Could it be that the compiled code from before the update was explicitly linked (with a baked-in rpath) against a different petsc to the one you have now? To assuage my fears, can you run firedrake/scripts/firedrake-clean to clear out the compiled kernel cache? And then run again without calling PetscInitialize/Finalize inside your kernel.
Cheers,
Lawrence -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.22 (GNU/Linux)
iQEcBAEBAgAGBQJWRF3NAAoJECOc1kQ8PEYvOuoIAMSymNuMusrzJOmYAFTnK+gR EPpuBzj3tCt5QIBSwGFuKgW8+IuWTmT4kY8t0gBfAx6Diw8RrQ6WSot9nKB/xdI2 EUTz9RuoJQBWbzgmvP1mxvnDfT+oulzj5arNEeFdHfxDrtaoKgwWU8PrLznegEGY e7Z9MJNH3S+DaamYeOF/gXVbNilrCHaKiKaKVsIMc2HYitBdOIwaJliOBTkFGwKC dnqLn2YCsX8CpPu9939UwxcLnG9BPfyFuEguSodMXVjZKAxLQh6DgIZ25U9F4q80 etmLHaubbgVTZfahMHH/5UxdLYr+ILzXmzikzQzkP20rrvZHPdauYhz5F4pGp4M= =0WNr -----END PGP SIGNATURE-----
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 12/11/15 10:01, Justin Chang wrote:
Strange, I had no issue running this on my Ubuntu system, it only happened with my MacOSX (El Capitan).
Perhaps I need to reinstall firedrake on my MacOSX and will holler if I get the same error again
Hmm, I initially tried on an ubuntu system, but just also tried on an El Cap mac build and could not reproduce there either. Cheers, Lawrence
Okay, I have no problems after reinstalling. I must have unintentionally done something taboo to my previous MacOSX firedrake installation. Strange On Thu, Nov 12, 2015 at 3:06 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
On 12/11/15 10:01, Justin Chang wrote:
Strange, I had no issue running this on my Ubuntu system, it only happened with my MacOSX (El Capitan).
Perhaps I need to reinstall firedrake on my MacOSX and will holler if I get the same error again
Hmm, I initially tried on an ubuntu system, but just also tried on an El Cap mac build and could not reproduce there either.
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 12/11/15 10:38, Justin Chang wrote:
Okay, I have no problems after reinstalling. I must have unintentionally done something taboo to my previous MacOSX firedrake installation.
OK, good. I note two things about your code. 1. If VI is on, you never call VecDestroy on ub. 2. If you don't want to create and destroy a SNES and Mat for every dof, you can do the following, and noting that the FormJacobian and FormFunction calls are the same for every kernel (it's just the data you feed them that changes): import ctypes class CtxPointer(ctypes.Structure): _fields_ = [("snes", ctypes.c_voidp), ("J", ctypes.c_voidp), ("setup", ctypes.c_int)] reaction_snes = PETSc.SNES().create(comm=PETSc.COMM_SELF) reaction_snes.setOptionsPrefix("reactions_") reaction_snes.setFromOptions() reaction_mat = PETSc.Mat().create(comm=PETSc.COMM_SELF) reaction_mat.setOptionsPrefix("reactions_") reaction_mat.setType(PETSc.Mat.Type.SEQDENSE) reaction_mat.setSizes((3, 3)) reaction_mat.setUp() reaction_mat.setFromOptions() AppCtx = CtxPointer() AppCtx.snes = reaction_snes.handle AppCtx.J = reaction_mat.handle AppCtx.setup = 0 ctx = op2.Global(1, data=ctypes.addressof(AppCtx), dtype=np.uintp) Add an extra argument to the reactions par_loop: op2.par_loop(reactions3_kernel,Q.dof_dset,u0_A.dat(op2.READ),u0_ modify your solve_kernel:B.dat(op2.READ), c_A.dat(op2.WRITE),c_B.dat(op2.WRITE), c_C.dat(op2.WRITE),k1.dat(op2.READ), ctx(op2.READ)) Abusing op2.Globals for this purpose is a bit icky. We should probably add the ability to pass in an AppCtx-like object, but I don't know exactly how this should interface with the python-level stuff. Now you need to modify your solve kernel: typedef struct { SNES snes; Mat J; int setup; } AppCtx; void solve_kernel(const double *psiA, const double *psiB, double *cA, double *cB, double *cC, const double *k1, const uintptr_t *ctx_) { AppCtx *ctx = (AppCtx *)(ctx_[0]); SNES snes = ctx->snes; Mat J = ctx->J; Vec x,r,lb,ub; const int numC = 3; const int VI = %(VI)d; const PetscScalar *xx; PetscInt i; PetscErrorCode ierr; /* Initialize */ problem_type = %(problem)d; VecCreateSeq(PETSC_COMM_SELF, numC, &x); VecDuplicate(x,&r); /* Initialize values */ VecSet(x,1); psiValues[0] = psiA[0]; psiValues[1] = psiB[0]; keq[0] = k1[0]; /* Solver */ if (!ctx->setup) { SNESSetFunction(snes,r,FormResidual,NULL); SNESSetJacobian(snes,J,J,FormJacobian,NULL); ctx->setup = 1; } if (VI) { VecDuplicate(x,&lb); VecDuplicate(x,&ub); VecSet(lb,0); VecSet(ub,2); SNESSetType(snes,SNESVINEWTONRSLS); SNESVISetVariableBounds(snes,lb,ub); } SNESSolve(snes,NULL,x); /* Store solution */ VecGetArrayRead(x,&xx); cA[0] = xx[0]; cB[0] = xx[1]; cC[0] = xx[2]; VecRestoreArrayRead(x,&xx); /* Destroy */ VecDestroy(&x); VecDestroy(&r); if (VI) { VecDestroy(&lb); VecDestroy(&ub); } } """ Now you reuse the same SNES and Mat every call, this way you don't build and destroy thousands of SNESes. And you can set up the solver in the usual way using the -reactions_ options prefix. Cheers, Lawrence
Got it, thanks! Would it also be wise to also reuse the Vecs r,x,lb,ub? Or is the performance difference going to be negligible for these? Thanks, Justin On Thu, Nov 12, 2015 at 4:29 AM, Lawrence Mitchell < lawrence.mitchell@imperial.ac.uk> wrote:
On 12/11/15 10:38, Justin Chang wrote:
Okay, I have no problems after reinstalling. I must have unintentionally done something taboo to my previous MacOSX firedrake installation.
OK, good.
I note two things about your code.
1. If VI is on, you never call VecDestroy on ub.
2. If you don't want to create and destroy a SNES and Mat for every dof, you can do the following, and noting that the FormJacobian and FormFunction calls are the same for every kernel (it's just the data you feed them that changes):
import ctypes
class CtxPointer(ctypes.Structure): _fields_ = [("snes", ctypes.c_voidp), ("J", ctypes.c_voidp), ("setup", ctypes.c_int)]
reaction_snes = PETSc.SNES().create(comm=PETSc.COMM_SELF) reaction_snes.setOptionsPrefix("reactions_") reaction_snes.setFromOptions() reaction_mat = PETSc.Mat().create(comm=PETSc.COMM_SELF) reaction_mat.setOptionsPrefix("reactions_") reaction_mat.setType(PETSc.Mat.Type.SEQDENSE) reaction_mat.setSizes((3, 3)) reaction_mat.setUp() reaction_mat.setFromOptions()
AppCtx = CtxPointer() AppCtx.snes = reaction_snes.handle AppCtx.J = reaction_mat.handle AppCtx.setup = 0
ctx = op2.Global(1, data=ctypes.addressof(AppCtx), dtype=np.uintp)
Add an extra argument to the reactions par_loop: op2.par_loop(reactions3_kernel,Q.dof_dset,u0_A.dat(op2.READ),u0_ modify your solve_kernel:B.dat(op2.READ), c_A.dat(op2.WRITE),c_B.dat(op2.WRITE), c_C.dat(op2.WRITE),k1.dat(op2.READ), ctx(op2.READ))
Abusing op2.Globals for this purpose is a bit icky. We should probably add the ability to pass in an AppCtx-like object, but I don't know exactly how this should interface with the python-level stuff.
Now you need to modify your solve kernel:
typedef struct { SNES snes; Mat J; int setup; } AppCtx;
void solve_kernel(const double *psiA, const double *psiB, double *cA, double *cB, double *cC, const double *k1, const uintptr_t *ctx_) { AppCtx *ctx = (AppCtx *)(ctx_[0]); SNES snes = ctx->snes; Mat J = ctx->J; Vec x,r,lb,ub; const int numC = 3; const int VI = %(VI)d; const PetscScalar *xx; PetscInt i; PetscErrorCode ierr; /* Initialize */ problem_type = %(problem)d; VecCreateSeq(PETSC_COMM_SELF, numC, &x); VecDuplicate(x,&r);
/* Initialize values */ VecSet(x,1); psiValues[0] = psiA[0]; psiValues[1] = psiB[0]; keq[0] = k1[0];
/* Solver */ if (!ctx->setup) { SNESSetFunction(snes,r,FormResidual,NULL); SNESSetJacobian(snes,J,J,FormJacobian,NULL); ctx->setup = 1; } if (VI) { VecDuplicate(x,&lb); VecDuplicate(x,&ub); VecSet(lb,0); VecSet(ub,2); SNESSetType(snes,SNESVINEWTONRSLS); SNESVISetVariableBounds(snes,lb,ub); } SNESSolve(snes,NULL,x);
/* Store solution */ VecGetArrayRead(x,&xx); cA[0] = xx[0]; cB[0] = xx[1]; cC[0] = xx[2]; VecRestoreArrayRead(x,&xx);
/* Destroy */ VecDestroy(&x); VecDestroy(&r); if (VI) { VecDestroy(&lb); VecDestroy(&ub); } }
"""
Now you reuse the same SNES and Mat every call, this way you don't build and destroy thousands of SNESes. And you can set up the solver in the usual way using the -reactions_ options prefix.
Cheers,
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Justin Chang
-
Lawrence Mitchell