Volume example-3 in Python

[ArshiaIlaty]

Hi libCEED team and @jeremylt

This PR is related to #1780 and is still a work in progress, but I’d really appreciate any early feedback you could provide.

While the general structure is in place, there are still areas I’m actively working on, especially around:

•	Vector dimensionality mismatch — I’m in the process of debugging and aligning these correctly.

              For 2D: Solution size (total)    : 263169
              Mesh coordinates vector size: 526338
              Element restriction expected size: 526338
              
              For 3D: Solution size (total)    : 274625
              Mesh coordinates vector size: 823875
              Element restriction expected size: 823875

•	Code optimization and cleanup — definitely room for improvement!

Before I invest more time into polishing and resolving all edge cases, I’d love a quick review or initial thoughts — particularly regarding design, integration approach, or anything that might be misaligned with libCEED’s principles or architecture.

Thanks so much in advance — looking forward to your feedback!

[jeremylt]

See

• https://github.com/CEED/libCEED/blob/main/python/tests/test-4-qfunction.py
• https://github.com/CEED/libCEED/blob/main/python/tests/setup-qfunctions.py
• https://github.com/CEED/libCEED/blob/main/python/tests/test-qfunctions.h
• https://github.com/CEED/libCEED/blob/main/python/tests/libceed-qfunctions.c

to see how creating a QFunction works

[ArshiaIlaty]

Thank you Jeremy, for guiding me. I have solved the issue, and the example should work in all dimensions. Please let me know if you have any suggestions. Does the below output look good to you?

Exact mesh volume    : 1
Computed mesh volume : 0.99999999999942
Volume error         : -5.7887028503956e-13

=== Dimension: 1, Degree: 3 ===

Exact mesh volume    : 1
Computed mesh volume : 0.99999999997779
Volume error         : -2.2206458893947e-11

=== Dimension: 1, Degree: 4 ===

Exact mesh volume    : 1
Computed mesh volume : 0.99999999997913
Volume error         : -2.0868640149274e-11

=== Dimension: 1, Degree: 5 ===

Exact mesh volume    : 1
Computed mesh volume : 0.99999999997596
Volume error         : -2.4043433910492e-11

=== Dimension: 2, Degree: 2 ===

Exact mesh volume    : 2.3561944901923
Computed mesh volume : 2.3561944901854
Volume error         : -6.9699801485967e-12

=== Dimension: 2, Degree: 3 ===

Exact mesh volume    : 2.3561944901923
Computed mesh volume : 2.356194490192
Volume error         : -3.1796787425264e-13

=== Dimension: 2, Degree: 4 ===

Exact mesh volume    : 2.3561944901923
Computed mesh volume : 2.356194490192
Volume error         : -2.9842794901924e-13

=== Dimension: 2, Degree: 5 ===

Exact mesh volume    : 2.3561944901923
Computed mesh volume : 2.3561944901917
Volume error         : -6.883382752676e-13

=== Dimension: 3, Degree: 2 ===

Exact mesh volume    : 2.3561944901923
Computed mesh volume : 2.3561944616941
Volume error         : -2.849826596929e-08

=== Dimension: 3, Degree: 3 ===

Exact mesh volume    : 2.3561944901923
Computed mesh volume : 2.3561944901644
Volume error         : -2.7991386986059e-11

=== Dimension: 3, Degree: 4 ===

Exact mesh volume    : 2.3561944901923
Computed mesh volume : 2.3561944901923
Volume error         : -5.0182080713057e-14

=== Dimension: 3, Degree: 5 ===

Exact mesh volume    : 2.3561944901923
Computed mesh volume : 2.3561944901922
Volume error         : -1.0036416142611e-13```

[jeremylt]

Suggested change

	main()
	main()

[jeremylt]

Cool, I forgot we added the ability to use native Python QFunction functions. I'll look sometime in the next few days. Hopefully I'll have an answer about how to test this file soon too.

[ArshiaIlaty]

Thank you so much for reviewing the code. Sure. Take your time

[jeremylt]

CI is giving some formatting problems - I would run make format-py to fix them and commit those fixes

[ArshiaIlaty]

CI is giving some formatting problems - I would run make format-py to fix them and commit those fixes

@jeremylt oh, sorry I should have done it myself. I truly appreciate it. Thanks for letting me know about it

[ArshiaIlaty]

I’ve used Black to reformat the code. Please let me know if you’d like me to commit and push the changes.

reformatted ex3-volume.py

All done! ✨ 🍰 ✨
1 file reformatted.```

[jeremylt]

You need to use the same formatter that our CI uses - make format-py runs autopep8. Without using autopep8, this file will fail our formatting checks

[ArshiaIlaty]

Got it. I made the changes. Please let me know if there are any other problems.

[jeremylt]

autopep8 still has changes:

diff --git a/examples/python/ex3_volume_qfunctions.py b/examples/python/ex3_volume_qfunctions.py
index 75e394f..c791f41 100644
--- a/examples/python/ex3_volume_qfunctions.py
+++ b/examples/python/ex3_volume_qfunctions.py
@@ -9,10 +9,10 @@
 
 def build_mass_diff(ctx, q, inputs, outputs):
     """Build quadrature data for a mass + diffusion operator.
-    
+
     At every quadrature point, compute w/det(J).adj(J).adj(J)^T and store
     the symmetric part of the result.
-    
+
     Args:
         ctx: QFunction context with 'dim' and 'space_dim' fields
         q: Number of quadrature points
@@ -25,14 +25,14 @@ def build_mass_diff(ctx, q, inputs, outputs):
     w = inputs[1]
     q_data = outputs[0]
     dim = ctx["dim"]
-    
+
     if dim == 1:
         for i in range(q):
             # Mass
             q_data[0][i] = w[i] * J[0][0][i]
             # Diffusion
             q_data[1][i] = w[i] / J[0][0][i]
-    
+
     elif dim == 2:
         for i in range(q):
             # J: 0 2   q_data: 0 2   adj(J):  J22 -J12
@@ -42,35 +42,35 @@ def build_mass_diff(ctx, q, inputs, outputs):
             J01 = J[1][0][i]
             J11 = J[1][1][i]
             qw = w[i] / (J00 * J11 - J10 * J01)
-            
+
             # Mass
             q_data[0][i] = w[i] * (J00 * J11 - J10 * J01)
             # Diffusion
             q_data[1][i] = qw * (J01 * J01 + J11 * J11)
             q_data[2][i] = qw * (J00 * J00 + J10 * J10)
             q_data[3][i] = -qw * (J00 * J01 + J10 * J11)
-    
+
     elif dim == 3:
         for i in range(q):
             # Compute the adjoint
             A = [[0 for _ in range(3)] for _ in range(3)]
-            
+
             for j in range(3):
                 for k in range(3):
                     # Equivalent code with J as a VLA and no mod operations:
                     # A[k][j] = J[j+1][k+1]*J[j+2][k+2] - J[j+1][k+2]*J[j+2][k+1]
                     A[k][j] = (
-                        J[(j+1)%3][(k+1)%3][i] * J[(j+2)%3][(k+2)%3][i] -
-                        J[(j+1)%3][(k+2)%3][i] * J[(j+2)%3][(k+1)%3][i]
+                        J[(j + 1) % 3][(k + 1) % 3][i] * J[(j + 2) % 3][(k + 2) % 3][i] -
+                        J[(j + 1) % 3][(k + 2) % 3][i] * J[(j + 2) % 3][(k + 1) % 3][i]
                     )
-            
+
             # Compute quadrature weight / det(J)
             det_j = J[0][0][i] * A[0][0] + J[0][1][i] * A[0][1] + J[0][2][i] * A[0][2]
             qw = w[i] / det_j
-            
+
             # Mass
             q_data[0][i] = w[i] * det_j
-            
+
             # Diffusion - stored in Voigt convention
             # 1 6 5
             # 6 2 4
@@ -85,11 +85,11 @@ def build_mass_diff(ctx, q, inputs, outputs):
 
 def build_mass_diff_single(ctx, q, inputs, outputs):
     """Build quadrature data for a mass + diffusion operator using a single component.
-    
+
     This QFunction is similar to build_mass_diff, but it works with a single component
     (x-coordinate) instead of all components. It assumes that the mesh is a Cartesian
     grid, so the Jacobians can be computed from the x-coordinates.
-    
+
     Args:
         ctx: QFunction context with 'dim' field
         q: Number of quadrature points
@@ -102,14 +102,14 @@ def build_mass_diff_single(ctx, q, inputs, outputs):
     w = inputs[1]
     q_data = outputs[0]
     dim = ctx["dim"]
-    
+
     if dim == 1:
         for i in range(q):
             # Mass
             q_data[0][i] = w[i] * x[0][i]
             # Diffusion
             q_data[1][i] = w[i] / x[0][i]
-    
+
     elif dim == 2:
         for i in range(q):
             # For a Cartesian grid, the Jacobian is diagonal
@@ -118,14 +118,14 @@ def build_mass_diff_single(ctx, q, inputs, outputs):
             J11 = x[0][i]
             det_j = J00 * J11
             qw = w[i] / det_j
-            
+
             # Mass
             q_data[0][i] = w[i] * det_j
             # Diffusion
             q_data[1][i] = qw * J11 * J11
             q_data[2][i] = qw * J00 * J00
             q_data[3][i] = 0.0
-    
+
     elif dim == 3:
         for i in range(q):
             # For a Cartesian grid, the Jacobian is diagonal
@@ -135,7 +135,7 @@ def build_mass_diff_single(ctx, q, inputs, outputs):
             J22 = x[0][i]
             det_j = J00 * J11 * J22
             qw = w[i] / det_j
-            
+
             # Mass
             q_data[0][i] = w[i] * det_j
             # Diffusion - stored in Voigt convention
@@ -152,9 +152,9 @@ def build_mass_diff_single(ctx, q, inputs, outputs):
 
 def apply_mass_diff(ctx, q, inputs, outputs):
     """Apply mass + diffusion operator.
-    
+
     Apply the action of the operator at quadrature points.
-    
+
     Args:
         ctx: QFunction context with 'dim' field
         q: Number of quadrature points
@@ -167,14 +167,14 @@ def apply_mass_diff(ctx, q, inputs, outputs):
     v = outputs[0]      # Result values
     vg = outputs[1]     # Result gradients
     dim = ctx["dim"]
-    
+
     if dim == 1:
         for i in range(q):
             # Mass
             v[0][i] = q_data[0][i] * u[0][i]
             # Diffusion
             vg[0][i] = q_data[1][i] * ug[0][i]
-    
+
     elif dim == 2:
         for i in range(q):
             # Mass
@@ -182,7 +182,7 @@ def apply_mass_diff(ctx, q, inputs, outputs):
             # Diffusion
             vg[0][i] = q_data[1][i] * ug[0][i] + q_data[3][i] * ug[1][i]
             vg[1][i] = q_data[3][i] * ug[0][i] + q_data[2][i] * ug[1][i]
-    
+
     elif dim == 3:
         for i in range(q):
             # Mass
@@ -190,4 +190,4 @@ def apply_mass_diff(ctx, q, inputs, outputs):
             # Diffusion
             vg[0][i] = q_data[1][i] * ug[0][i] + q_data[6][i] * ug[1][i] + q_data[5][i] * ug[2][i]
             vg[1][i] = q_data[6][i] * ug[0][i] + q_data[2][i] * ug[1][i] + q_data[4][i] * ug[2][i]
-            vg[2][i] = q_data[5][i] * ug[0][i] + q_data[4][i] * ug[1][i] + q_data[3][i] * ug[2][i]
\ No newline at end of file
+            vg[2][i] = q_data[5][i] * ug[0][i] + q_data[4][i] * ug[1][i] + q_data[3][i] * ug[2][i]

[ArshiaIlaty]

Since I didn’t use ex3_volume_qfunctions.py in the example file, I believe it’s not necessary to include it in the PR, and reviewing it may not be relevant.

I’m thinking of removing it from the PR to keep things clean—what do you think?

[ArshiaIlaty]

Hi @jeremylt , I hope you are doing well. I have tried to address your comments and made the changes. Please let me know what steps I need to take to improve this version. Thanks a lot. I appreciate your time.

└── examples/
    └── python/
        ├── ex3-volume.py          # Main example implementation
        ├── ex3-volume.h           # QFunction definitions
        ├── qfunctions.c           # QFunction wrapper
        ├── setup-qfunctions.py    # Build configuration
        └── test_ex3_volume.py     # Test suite```

[jeremylt]

I'll take a closer look tomorrow, but at first glance this looks pretty good

[jeremylt]

You have both Python and C formatting issues according to CI

[ArshiaIlaty]

You have both Python and C formatting issues according to CI

I used Make Format C and Python for all directories. I wonder where the issue comes from. Should I redo it?

[jeremylt]

Ok, we're close here. Open todos

1. Volume example-3 in Python #1805 (comment)
2. Volume example-3 in Python #1805 (comment)
3. remove extra files
4. style fixes

[ArshiaIlaty]

Ok, we're close here. Open todos

1. Volume example-3 in Python #1805 (comment)
2. remove extra files
3. style fixes

Fantastic. Would you please let me know which files should get removed? I guess the only one is test_ex3_volume.py. Am I right?

In case of style fixing, do you have any suggestions? I ran make format-c and py multiple times and did not see any modifications in the code.

[jeremylt]

Removal - we don't need ex3_volume_qfunctions.py or setup.py

Style - the changes are automatically inserted into the changed files. I usually do

$ make format
$ git add .
$ git commit -m "style fixes"

then push

Also, added one more comment: #1805 (comment)

[ArshiaIlaty]

Hi @jeremylt Hope you are doing well. I checked the code and reformatted it. Please let me know if you have any suggestions. I think it should pass the CI tests

[jeremylt]

Between the three C examples, this one had some especially tricky details. Nice work. I'll merge this to the staging branch. Thanks for you work on this - this PR is accepted

CC @valeriabarra

[ArshiaIlaty]

@jeremylt Thank you so much. I appreciate your kind words. Happy to see my contribution got accepted.

[jeremylt]

Manually merged into the staging branch