Python PETSc Example BP1 by tylercollins5737 · Pull Request #1806 · CEED/libCEED

tylercollins5737 · 2025-04-19T09:04:56Z

Hello!

This is a Python implementation of bpsraw.c for BP1 and was initially discussed in #1774 . I'm aiming to recreate the benchmark behavior using libceed and petsc4py. I appreciate any feedback or indications I'm headed in the right direction!

It does run, but this is a work in progress. Remaining tasks:

jeremylt · 2025-04-21T14:54:51Z

+    matrix = PETSc.Mat().createAIJ([num_dofs, num_dofs], comm=PETSc.COMM_WORLD)
+    matrix.setFromOptions()
+    matrix.setUp()
+    for i in range(num_dofs):
+        e_vec = ceed.Vector(num_dofs)
+        r_vec = ceed.Vector(num_dofs)
+        e_vec.set_value(0.0)
+        r_vec.set_value(0.0)
+        with e_vec.array_write() as arr:
+            arr[i] = 1.0
+        op.apply(e_vec, r_vec)
+        with r_vec.array_read() as r_arr:
+            matrix.setValues(range(num_dofs), [i], r_arr)
+    matrix.assemble()


You'll want a MatShell, not a MatAIJ

Apologies for leaving this a bit stale. Addressing this led to a lot of other changes. I had trouble setting the matrix type to "shell" and applying context, kept getting segmentation faults. So, instead of using the "shell" type, I used the "python" type and set the context through a python class. Pulled that implementation from the petsc4py 2DPoisson example (about half way down) --> link. I'll tidy up what I have a submit an update soon.

Overall this is looking like the right approach. The big thing is to try to avoid copying around data during frequent operations, like the matrix application

Thank you for all the suggestions! Definitely, keeping the data in place would mean less overhead and better efficiency

valeriabarra · 2025-05-20T16:34:25Z

@tylercollins5737 , thank you for your work. Are you able to address the reviewer's comments?

tylercollins5737 · 2025-05-20T22:54:10Z

@valeriabarra I was able to get it working using the "python" matrix type, MatAIJ was not the best. I believe the "python" type and MatCreateShell can achieve the same thing.

jeremylt · 2025-05-21T16:57:06Z

+        with self.v.array_read() as va:
+            y.setValues(range(lo, hi), va[lo:hi])
+        y.assemble()


Ideally, you'd set v to use y's array before you apply the libCEED operator

That makes more sense, thank you!

jeremylt · 2025-05-21T16:58:04Z

+def CreateRestriction(ceed, mesh_elem, P, num_comp, ndofs, offsets):
+    return ceed.ElemRestriction(mesh_elem, P, num_comp, 1, ndofs, offsets, cmode=libceed.USE_POINTER)


If this just wraps another function call, I think its clearer to use the function call itself

Totally agree, I have another CreateRestriction function that reflects C but its a WIP

valeriabarra · 2025-10-08T21:08:57Z

Hey @tylercollins5737 , how are you? I was wondering if you intended to complete this PR or if there were any impeding issues that you encountered.

There are a couple of other potential contributors that are eyeing this project and I was wondering if there was a missing interface for which this BP could not be completed or if it was just for personal reasons. Thank you!

tylercollins5737 · 2025-10-09T22:45:19Z

Hey @tylercollins5737 , how are you? I was wondering if you intended to complete this PR or if there were any impeding issues that you encountered.

There are a couple of other potential contributors that are eyeing this project and I was wondering if there was a missing interface for which this BP could not be completed or if it was just for personal reasons. Thank you!

Hey @valeriabarra hope you are doing well! In short, progress didn't stop because of any major limitations in the libceed or petsc libraries. But the qfunction for building the RHS "SetupMassRHS" and "Error" I believe aren't available so a workaround would have to be implemented which could be very slow, or an addition to the gallery of qfunctions, or someway to use the qfunctions used by the BPs, or a way to define a custom one in python. If there are interested contributors they're more than welcome to work on it. If no one picks it up, I'll start on it again and get the development environment setup on a computer I actually own. Cheers!

jeremylt · 2025-10-10T14:45:35Z

Yeah, since there are QFunctions that are not available via the gallery, you'd have to use the technique shown in this folder for ex3_volume.py (and associated files) to use a QFunction written in C.

tylercollins5737 · 2025-10-10T16:39:26Z

@jeremylt Cool, yeah that should do the trick.

jeremylt · 2025-12-16T16:42:35Z

Hi @tylercollins5737, I'm hoping to have this ready in case next semester's students want to build on this for BP 2-6. This is good work so far and will be very helpful to future users. Is there anything I can do to help here? I'm happy to help out on any pieces that have you stuck so we can get this moving again.

tylercollins5737 · 2025-12-18T03:38:47Z

Hi @tylercollins5737, I'm hoping to have this ready in case next semester's students want to build on this for BP 2-6. This is good work so far and will be very helpful to future users. Is there anything I can do to help here? I'm happy to help out on any pieces that have you stuck so we can get this moving again.

Cool, thank you for the support. I have the development environment setup again with everything updated, and I'm making progress on the example! I will reach out if I get stumped.

zatkins-dev · 2026-07-08T21:26:47Z

@tylercollins5737 are you planning on coming back to this? I think there's a lot of good work here! If not, it may be time to close this PR.

tylercollins5737 · 2026-07-09T00:54:58Z

@tylercollins5737 are you planning on coming back to this? I think there's a lot of good work here! If not, it may be time to close this PR.

@zatkins-dev Yes! I should have have a pretty complete update soon, ideally by the end of the weekend.

zatkins-dev · 2026-07-09T02:29:38Z

Awesome! Happy to keep it open then, thanks!

The gallery has Mass3DBuild and MassApply, but nothing like SetupMassRhs or Error, which depend on the exact solution the benchmark solves for. These follow the ex1-ex3 pattern: C sources next to the Python examples, built into the QFunctions library by qfunctions.c.

Port examples/petsc/bpsraw.c to Python with petsc4py and the libCEED Python bindings, solving BP1. A python PETSc Mat takes the place of the MatShell; its mult() runs the libCEED operator. MatMult hands the PETSc Vec arrays straight to libCEED with CEED_USE_POINTER instead of copying them, on the host and on the device. petsc4py has no VecGetArrayAndMemType, so on the device the pointer comes from getCUDAHandle and is passed through the CUDA array interface.

The example picks its iteration limit from how long a first solve takes, so the tests pin it with ksp_max_it_clip to keep the iteration count the same every run, as the C test does.

tylercollins5737 · 2026-07-13T00:49:06Z

QFunctions are in the PR now (qfunctions/bps/bp1.h, common.h, plus qfunctions.c), using the ex3_volume.py approach you pointed me to. The gallery has Mass3DBuild/MassApply but nothing like SetupMassRhs or Error, so those have to be C QFunctions.
MatMult no longer copies data. The PETSc Vec arrays go to libCEED with CEED_USE_POINTER, on the host and on the device.
Added tests (BP1, degrees 1–3) and matched the CLI to bpsraw.c (-problem bp1, -ksp_max_it_clip, -write_solution).
Output matches bpsraw.c on CPU and GPU, across 1, 2 and 4 ranks and degrees 1–4. For example at 1 rank, degree 3: 16 iterations, max error 1.730043e-09, same as the C.

Two notes:

I copied bp1.h/common.h into examples/python/qfunctions/bps/, the way ex1–ex3 keep their own copies, so the example is self contained.
CeedVectorTakeArray isn't wrapped in the Python bindings, so the borrowed arrays get rebound on the next MatMult instead of detached. Works as written, but not the same as the C.

jeremylt · 2026-07-13T06:37:29Z

 import ex1_volume
 import ex2_surface
 import ex3_volume
+import bpsraw


since this file relies upon PETSc4Py, we'll want to test it separately from ex1-3, which only need libCEED

jeremylt · 2026-07-13T06:40:29Z

+        # Apply libCEED operator
+        self.op_apply.apply(self.x_ceed, self.y_ceed)
+
+        # Restore arrays; the C also detaches them with CeedVectorTakeArray, which the Python bindings do not


We really should just add CeedVectorTakeArarry to the Python bindings - @zatkins-dev or I can hit that, though dunno how long it'll take to publish that

jeremylt · 2026-07-13T06:43:22Z

+                         if isinstance(value, int)}
+
+
+def Split3(size, reverse=False):


we should use the same function naming, style throughout. Also a brief docstring for each would be handy (though BPsRaw probably omitted that? If so, then no worries here)

jeremylt · 2026-07-13T06:44:52Z

I have just a handful of comments, but this looks great!

jeremylt reviewed Apr 21, 2025

View reviewed changes

jeremylt reviewed May 21, 2025

View reviewed changes

Comment thread examples/python/bpsraw.py Outdated

jeremylt reviewed May 21, 2025

View reviewed changes

jeremylt added examples PETSc 0-WIP labels May 28, 2025

tylercollins5737 added 3 commits July 12, 2026 17:07

python - test BP1 at degrees 1 to 3

f833966

The example picks its iteration limit from how long a first solve takes, so the tests pin it with ksp_max_it_clip to keep the iteration count the same every run, as the C test does.

tylercollins5737 force-pushed the tyler/bp1_python branch from c387450 to f833966 Compare July 13, 2026 00:26

jeremylt reviewed Jul 13, 2026

View reviewed changes

		def CreateRestriction(ceed, mesh_elem, P, num_comp, ndofs, offsets):
		return ceed.ElemRestriction(mesh_elem, P, num_comp, 1, ndofs, offsets, cmode=libceed.USE_POINTER)

Uh oh!

Conversation

tylercollins5737 commented Apr 19, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

valeriabarra commented May 20, 2025

Uh oh!

tylercollins5737 commented May 20, 2025

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

valeriabarra commented Oct 8, 2025

Uh oh!

tylercollins5737 commented Oct 9, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jeremylt commented Oct 10, 2025

Uh oh!

tylercollins5737 commented Oct 10, 2025

Uh oh!

jeremylt commented Dec 16, 2025

Uh oh!

tylercollins5737 commented Dec 18, 2025

Uh oh!

zatkins-dev commented Jul 8, 2026

Uh oh!

tylercollins5737 commented Jul 9, 2026

Uh oh!

zatkins-dev commented Jul 9, 2026

Uh oh!

tylercollins5737 commented Jul 13, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jeremylt commented Jul 13, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

tylercollins5737 commented Apr 19, 2025 •

edited

Loading

tylercollins5737 commented Oct 9, 2025 •

edited

Loading