Kartiki p ex3 volume julia

[KartikiP]

This is an effort at contributing a Julia version of Example 3 - volume.

[jeremylt]

This is looking pretty good. There's a few comments that I'll have to make this a bit more idiomatic and tie into the test system, but this is 85% good to go. Thanks for the good work!

[jeremylt]

Note that the linting job has some formatting changes:
https://github.com/CEED/libCEED/actions/runs/19895474658/job/57042483042#step:4:1

[jeremylt]

Note that the linting job has some formatting changes: https://github.com/CEED/libCEED/actions/runs/19895474658/job/57042483042#step:4:1

This is the command that we use in GitHub Actions to run the formatter and check if any changes are needed:

libCEED/.github/workflows/julia-test-with-style.yml

Line 41 in ce3aeaf

julia --project=.style/ -e 'import Pkg; Pkg.instantiate()' && julia --project=.style/ .style/ceed_style.jl && git diff --exit-code src test examples

[KartikiP]

Thank you so much for your review. I will take a look, try and fix the formatting and remove the part mentioned asap.

[jeremylt]

When I run

julia --project=.style/ -e 'import Pkg; Pkg.instantiate()' && julia --project=.style/ .style/ceed_style.jl && git diff --exit-code src test examples 

I get this diff

diff --git a/julia/LibCEED.jl/examples/ex3-volume.jl b/julia/LibCEED.jl/examples/ex3-volume.jl
index 3dfd68f1..e1f11462 100644
--- a/julia/LibCEED.jl/examples/ex3-volume.jl
+++ b/julia/LibCEED.jl/examples/ex3-volume.jl
@@ -3,7 +3,7 @@ using LibCEED, Printf
 include("common.jl")
 
 function transform_mesh_coords!(dim, mesh_size, mesh_coords)
-     @witharray coords = mesh_coords begin
+    @witharray coords = mesh_coords begin
         if dim == 1
             for i = 1:mesh_size
                 # map [0,1] to [0,1] varying the mesh density
@@ -34,31 +34,38 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
     prob_size < 0 && (prob_size = 256*1024)
 
     ceed = Ceed(ceed_spec)
-    mesh_basis = 
-    create_tensor_h1_lagrange_basis(ceed, dim, ncompx, mesh_order + 1, num_qpts, GAUSS)
-    sol_basis = 
-    create_tensor_h1_lagrange_basis(ceed, dim, 1, sol_order + 1, num_qpts, GAUSS)
+    mesh_basis =
+        create_tensor_h1_lagrange_basis(ceed, dim, ncompx, mesh_order + 1, num_qpts, GAUSS)
+    sol_basis =
+        create_tensor_h1_lagrange_basis(ceed, dim, 1, sol_order + 1, num_qpts, GAUSS)
 
-    nxyz = get_cartesian_mesh_size(dim, sol_order, prob_size) 
-    println("Mesh size:",nxyz)
+    nxyz = get_cartesian_mesh_size(dim, sol_order, prob_size)
+    println("Mesh size:", nxyz)
 
     #Build CeedElemRestriction objects describing the mesh and solution discrete
     # mesh_rstr: for building (no qdata restriction needed)
-    mesh_size, mesh_rstr, _ = build_cartesian_restriction(ceed, dim, nxyz, mesh_order, ncompx, num_qpts)
+    mesh_size, mesh_rstr, _ =
+        build_cartesian_restriction(ceed, dim, nxyz, mesh_order, ncompx, num_qpts)
 
-   # sol_rstr + sol_rstr_i: for solving
+    # sol_rstr + sol_rstr_i: for solving
     sol_size, sol_rstr, sol_rstr_i = build_cartesian_restriction(
-        ceed, dim, nxyz, sol_order, 1, num_qpts, mode=RestrictionAndStrided
+        ceed,
+        dim,
+        nxyz,
+        sol_order,
+        1,
+        num_qpts,
+        mode=RestrictionAndStrided,
     )
 
-   # mesh_coords
+    # mesh_coords
     mesh_coords = CeedVector(ceed, mesh_size)
     set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords)
     exact_vol = transform_mesh_coords!(dim, mesh_size, mesh_coords)
 
     #Create the Q-function that builds the mass operator ( i.e it computes the quadrature data) and set its context data.
-    num_q_comp = 1 + div(dim*(dim+1), 2)
-    
+    num_q_comp = 1 + div(dim*(dim + 1), 2)
+
     @interior_qf build_qfunc = (
         ceed,
         dim=dim,
@@ -68,15 +75,15 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
         begin
             # Compute determinant
             det_J = det(dx)
-            
+
             # Store mass component
-            qdata[1] = weights * det_J
-            
+            qdata[1] = weights*det_J
+
             # Store diffusion components (J^T * J)
             idx = 2
             for i = 1:dim
                 for j = i:dim
-                    qdata[idx] = dx[:, i]' * dx[:, j]
+                    qdata[idx] = dx[:, i]'*dx[:, j]
                     idx += 1
                 end
             end
@@ -97,7 +104,7 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
     # Apply to get qdata
     elem_qpts = num_qpts^dim
     num_elem = prod(nxyz)
-    qdata = CeedVector(ceed, num_elem * elem_qpts * num_q_comp)
+    qdata = CeedVector(ceed, num_elem*elem_qpts*num_q_comp)
     print("Computing the quadrature data for the mass operator ...")
     flush(stdout)
     apply!(build_oper, mesh_coords, qdata)
@@ -114,12 +121,12 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
         (dv, :out, EVAL_GRAD, dim),
         begin
             # Apply mass: v = qdata[1] * u
-            v .= qdata[1] .* u
-            
+            v .= qdata[1].*u
+
             # Apply diffusion: dv = (qdata[2:end]) * du
             # The qdata contains the symmetric diffusion tensor (J^T*J)
             # dv_i = sum_j (J^T*J)_{i,j} * du_j
-            
+
             # For efficiency, rebuild the matrix from stored components
             idx = 2
             for i = 1:dim
@@ -127,42 +134,42 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
                 for j = 1:dim
                     # Reconstruct symmetric matrix element
                     if j >= i
-                        mat_idx = idx + div((j-1)*j, 2) + (i-1)
+                        mat_idx = idx + div((j - 1)*j, 2) + (i - 1)
                     else
-                        mat_idx = idx + div((i-1)*i, 2) + (j-1)
+                        mat_idx = idx + div((i - 1)*i, 2) + (j - 1)
                     end
-                    dv_i += qdata[mat_idx] * du[j]
+                    dv_i += qdata[mat_idx]*du[j]
                 end
                 dv[i] = dv_i
             end
         end,
     )
     apply_oper = Operator(
-    ceed,
-    qf=apply_qfunc,
-    fields=[
-        (:u, sol_rstr, sol_basis, CeedVectorActive()),
-        (:du, sol_rstr, sol_basis, CeedVectorActive()),
-        (:qdata, sol_rstr_i, BasisNone(), qdata),
-        (:v, sol_rstr, sol_basis, CeedVectorActive()),
-        (:dv, sol_rstr, sol_basis, CeedVectorActive()),
-    ],
-)
+        ceed,
+        qf=apply_qfunc,
+        fields=[
+            (:u, sol_rstr, sol_basis, CeedVectorActive()),
+            (:du, sol_rstr, sol_basis, CeedVectorActive()),
+            (:qdata, sol_rstr_i, BasisNone(), qdata),
+            (:v, sol_rstr, sol_basis, CeedVectorActive()),
+            (:dv, sol_rstr, sol_basis, CeedVectorActive()),
+        ],
+    )
 
-# # Compute the mesh volume using the massdiff operator
+    # # Compute the mesh volume using the massdiff operator
     print("Computing the mesh volume using the formula: vol = 1^T * (M + K) * 1...")
     flush(stdout)
-    
+
     u = CeedVector(ceed, sol_size)
     v = CeedVector(ceed, sol_size)
     u[] = 1.0
-    
+
     # Apply operator
     apply!(apply_oper, u, v)
-    
+
     # Compute volume
     vol = witharray_read(sum, v, MEM_HOST)
-    
+
     @printf("Exact mesh volume    : % .14g\n", exact_vol)
     @printf("Computed mesh volume : % .14g\n", vol)
     @printf("Volume error         : % .14g\n", vol - exact_vol)

[jeremylt]

If you're unfamiliar with git diff, here's a brief explainer on how it works: https://stackoverflow.com/questions/12320863/how-do-you-take-a-git-diff-file-and-apply-it-to-a-local-branch-that-is-a-copy-o

Basically, copy the diff from my message above into a file foo.diff and put it into the root of your libCEED dir, then run git apply foo.diff and commit the changes.

[jeremylt]

index 8f6a93c9..e17d7974 100644
--- a/julia/LibCEED.jl/examples/ex3-volume.jl
+++ b/julia/LibCEED.jl/examples/ex3-volume.jl
@@ -1,7 +1,7 @@
 using LibCEED, Printf
 
 include("common.jl")   # defines helper functions: create_tensor_h1_lagrange_basis,
-                       # build_cartesian_restriction, get_cartesian_mesh_size, etc.
+# build_cartesian_restriction, get_cartesian_mesh_size, etc.
 
 function transform_mesh_coords!(dim, mesh_size, mesh_coords)
     @witharray coords = mesh_coords begin
@@ -11,45 +11,59 @@ function transform_mesh_coords!(dim, mesh_size, mesh_coords)
             end
             exact_volume = 1.0
         else
-            num_nodes = mesh_size ÷ dim
+            num_nodes = mesh_size÷dim
             @inbounds @simd for i = 1:num_nodes
                 u = coords[i]            # ξ ∈ [0,1]
                 v = coords[i+num_nodes]  # η ∈ [0,1]
                 r = 1.0 + u
-                θ = (pi/2) * v
-                coords[i]           = r * cos(θ)
-                coords[i+num_nodes] = r * sin(θ)
+                θ = (pi/2)*v
+                coords[i] = r*cos(θ)
+                coords[i+num_nodes] = r*sin(θ)
             end
-            exact_volume = 3.0/4.0 * pi
+            exact_volume = 3.0/4.0*pi
         end
         return exact_volume
     end
 end
 
-function run_ex3(; ceed_spec="/cpu/self", dim=2, mesh_order=3, sol_order=3,
-                 num_qpts=4, prob_size=-1, gallery=false)
-
+function run_ex3(;
+    ceed_spec="/cpu/self",
+    dim=2,
+    mesh_order=3,
+    sol_order=3,
+    num_qpts=4,
+    prob_size=-1,
+    gallery=false,
+)
     prob_size < 0 && (prob_size = 256*1024)
     ncompx = dim
 
     ceed = Ceed(ceed_spec)
 
-    
     # Bases
-    
-    mesh_basis = create_tensor_h1_lagrange_basis(ceed, dim, ncompx, mesh_order + 1, num_qpts, GAUSS)
-    sol_basis  = create_tensor_h1_lagrange_basis(ceed, dim, 1,       sol_order + 1,  num_qpts, GAUSS)
 
-  
+    mesh_basis =
+        create_tensor_h1_lagrange_basis(ceed, dim, ncompx, mesh_order + 1, num_qpts, GAUSS)
+    sol_basis =
+        create_tensor_h1_lagrange_basis(ceed, dim, 1, sol_order + 1, num_qpts, GAUSS)
+
     # Mesh size & restrictions
-  
+
     nxyz = get_cartesian_mesh_size(dim, sol_order, prob_size)
     println("Mesh size: ", nxyz)
 
-    mesh_size, mesh_rstr, _ = build_cartesian_restriction(ceed, dim, nxyz, mesh_order, ncompx, num_qpts)
+    mesh_size, mesh_rstr, _ =
+        build_cartesian_restriction(ceed, dim, nxyz, mesh_order, ncompx, num_qpts)
 
     sol_size, sol_rstr, _ = build_cartesian_restriction(
-        ceed, dim, nxyz, sol_order, 1, num_qpts, mode=RestrictionAndStrided)
+        ceed,
+        dim,
+        nxyz,
+        sol_order,
+        1,
+        num_qpts,
+        mode=RestrictionAndStrided,
+    )
 
     # Physical mesh coordinates + exact volume
 
@@ -57,12 +71,11 @@ function run_ex3(; ceed_spec="/cpu/self", dim=2, mesh_order=3, sol_order=3,
     set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords)
     exact_vol = transform_mesh_coords!(dim, mesh_size, mesh_coords)
 
-  
     num_q_comp = 1 + div(dim*(dim + 1), 2)   # 1 for detJ (mass), rest for symmetric JᵀJ
-    num_elem   = prod(nxyz)
-    elem_qpts  = num_qpts^dim
+    num_elem = prod(nxyz)
+    elem_qpts = num_qpts^dim
 
-    qdata_size = num_elem * elem_qpts * num_q_comp
+    qdata_size = num_elem*elem_qpts*num_q_comp
     qdata = CeedVector(ceed, qdata_size)
     qdata[] = 0.0
 
@@ -73,29 +86,30 @@ function run_ex3(; ceed_spec="/cpu/self", dim=2, mesh_order=3, sol_order=3,
         elem_qpts,          # nodes per element = quadrature points
         num_q_comp,         # components per quadrature point
         qdata_size,
-        CeedInt[1, num_q_comp, num_q_comp * elem_qpts]  # strides: [comp, elem, total]
+        CeedInt[1, num_q_comp, num_q_comp*elem_qpts],  # strides: [comp, elem, total]
     )
 
-    
     @interior_qf build_qfunc = (
         ceed,
         dim=dim,
-        (dx,      :in,  EVAL_GRAD,    dim, dim),
-        (weights, :in,  EVAL_WEIGHT),
-        (qdata,   :out, EVAL_NONE,   num_q_comp),
+        (dx, :in, EVAL_GRAD, dim, dim),
+        (weights, :in, EVAL_WEIGHT),
+        (qdata, :out, EVAL_NONE, num_q_comp),
         begin
             # det(J)
             if dim == 1
-                detJ = dx[1,1]
+                detJ = dx[1, 1]
             elseif dim == 2
-                detJ = dx[1,1]*dx[2,2] - dx[1,2]*dx[2,1]
+                detJ = dx[1, 1]*dx[2, 2] - dx[1, 2]*dx[2, 1]
             else
-                detJ = (dx[1,1]*(dx[2,2]*dx[3,3] - dx[2,3]*dx[3,2]) -
-                        dx[1,2]*(dx[2,1]*dx[3,3] - dx[2,3]*dx[3,1]) +
-                        dx[1,3]*(dx[2,1]*dx[3,2] - dx[2,2]*dx[3,1]))
+                detJ = (
+                    dx[1, 1]*(dx[2, 2]*dx[3, 3] - dx[2, 3]*dx[3, 2]) -
+                    dx[1, 2]*(dx[2, 1]*dx[3, 3] - dx[2, 3]*dx[3, 1]) +
+                    dx[1, 3]*(dx[2, 1]*dx[3, 2] - dx[2, 2]*dx[3, 1])
+                )
             end
 
-            qdata[1] = weights * detJ   # mass contribution
+            qdata[1] = weights*detJ   # mass contribution
 
             # Upper triangle of JᵀJ (symmetric storage)
             k = 2
@@ -103,22 +117,22 @@ function run_ex3(; ceed_spec="/cpu/self", dim=2, mesh_order=3, sol_order=3,
                 for j = i:dim
                     s = zero(eltype(dx))
                     @inbounds @simd for m = 1:dim
-                        s += dx[m,i] * dx[m,j]
+                        s += dx[m, i]*dx[m, j]
                     end
                     qdata[k] = s
                     k += 1
                 end
             end
-        end
+        end,
     )
 
     build_oper = Operator(
         ceed,
         qf=build_qfunc,
         fields=[
-            (:dx,      mesh_rstr, mesh_basis, mesh_coords),
+            (:dx, mesh_rstr, mesh_basis, mesh_coords),
             (:weights, ElemRestrictionNone(), mesh_basis, CeedVectorNone()),
-            (:qdata,   qdata_rstr, BasisNone(), qdata),
+            (:qdata, qdata_rstr, BasisNone(), qdata),
         ],
     )
 
@@ -130,13 +144,13 @@ function run_ex3(; ceed_spec="/cpu/self", dim=2, mesh_order=3, sol_order=3,
     @interior_qf apply_qfunc = (
         ceed,
         dim=dim,
-        (u,     :in,  EVAL_INTERP),
-        (du,    :in,  EVAL_GRAD, dim),
-        (qdata, :in,  EVAL_NONE, num_q_comp),
-        (v,     :out, EVAL_INTERP),
-        (dv,    :out, EVAL_GRAD, dim),
+        (u, :in, EVAL_INTERP),
+        (du, :in, EVAL_GRAD, dim),
+        (qdata, :in, EVAL_NONE, num_q_comp),
+        (v, :out, EVAL_INTERP),
+        (dv, :out, EVAL_GRAD, dim),
         begin
-            v .= qdata[1] .* u                     # mass part
+            v .= qdata[1].*u                     # mass part
 
             # diffusion part: reconstruct symmetric matrix from upper triangle
             idx = 2
@@ -144,30 +158,29 @@ function run_ex3(; ceed_spec="/cpu/self", dim=2, mesh_order=3, sol_order=3,
                 acc = zero(eltype(dv))
                 for j = 1:dim
                     if j >= i
-                        mat_idx = idx + (j*(j-1))÷2 + (i-1)
+                        mat_idx = idx + (j*(j - 1))÷2 + (i - 1)
                     else
-                        mat_idx = idx + (i*(i-1))÷2 + (j-1)
+                        mat_idx = idx + (i*(i - 1))÷2 + (j - 1)
                     end
-                    acc += qdata[mat_idx] * du[j]
+                    acc += qdata[mat_idx]*du[j]
                 end
                 dv[i] = acc
             end
-        end
+        end,
     )
 
     apply_oper = Operator(
         ceed,
         qf=apply_qfunc,
         fields=[
-            (:u,     sol_rstr,  sol_basis,  CeedVectorActive()),
-            (:du,    sol_rstr,  sol_basis,  CeedVectorActive()),
+            (:u, sol_rstr, sol_basis, CeedVectorActive()),
+            (:du, sol_rstr, sol_basis, CeedVectorActive()),
             (:qdata, qdata_rstr, BasisNone(), qdata),
-            (:v,     sol_rstr,  sol_basis,  CeedVectorActive()),
-            (:dv,    sol_rstr,  sol_basis,  CeedVectorActive()),
+            (:v, sol_rstr, sol_basis, CeedVectorActive()),
+            (:dv, sol_rstr, sol_basis, CeedVectorActive()),
         ],
     )
 
-    
     print("Computing volume via 1ᵀ (M + K) 1 ... ")
     flush(stdout)
     u_vec = CeedVector(ceed, sol_size)
@@ -182,14 +195,13 @@ function run_ex3(; ceed_spec="/cpu/self", dim=2, mesh_order=3, sol_order=3,
     @printf("Volume error         : %.14g\n", computed_vol - exact_vol)
 end
 
-
 # Run
 
 run_ex3(
-    ceed_spec  = "/cpu/self",   # use "/gpu/cuda" or "/gpu/hip" if available
-    dim        = 2,
-    mesh_order = 3,
-    sol_order  = 3,
-    num_qpts   = 4,
-    prob_size  = -1,            # -1 → auto-size (~256K dofs)
+    ceed_spec="/cpu/self",   # use "/gpu/cuda" or "/gpu/hip" if available
+    dim=2,
+    mesh_order=3,
+    sol_order=3,
+    num_qpts=4,
+    prob_size=-1,            # -1 → auto-size (~256K dofs)
 )

[jeremylt]

WRT the comments you added - I think its best to mirror the amount and location of comments found in ex1-volume.jl and ex2-surface.jl. Consistency is important here

[jeremylt]

Suggested change

	function run_ex3(; ceed_spec="/cpu/self", dim=2, mesh_order=3, sol_order=3,
	num_qpts=4, prob_size=-1, gallery=false)
	function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, gallery)

This is the way the signature is given for ex1 and ex2

[jeremylt]

Honestly, I'd ditch 5c36feb - it introduced a lot of changes that push your code further from the other examples. I would instead apply the first diff I gave you.

[KartikiP]

I am getting this error when I try to apply the diff
(base) kartikipande@Kartikis-MacBook-Pro libCEED % git apply foo.diff
error: patch failed: julia/LibCEED.jl/examples/ex3-volume.jl:34
error: julia/LibCEED.jl/examples/ex3-volume.jl: patch does not apply

[jeremylt]

Did you remove 5c36feb ?

That patch will only work if your branch is on commit de961ac

[jeremylt]

See here for how to remove the last commit in a branch

https://stackoverflow.com/questions/1338728/how-do-i-delete-a-commit-from-a-branch

git reset --hard HEAD~1

(Note, this will require a force push, which usually is one of those "know you absolutely want to do this" situations, but destroying the most recent commit from a branch is one of those usages that's less drastic)

[jeremylt]

Huzzah, CI is happy. You accidentally committed foo.diff though, so that should be removed

[jeremylt]

Ok, we are close but you had some differences when you ran ex1 vs ex3 when the outputs should be very similar. I went through and identified all differences between your code and ex1 and modified them for the different qdata sizes and came up with this:

index e1f11462..68edf598 100644
--- a/julia/LibCEED.jl/examples/ex3-volume.jl
+++ b/julia/LibCEED.jl/examples/ex3-volume.jl
@@ -28,8 +28,7 @@ function transform_mesh_coords!(dim, mesh_size, mesh_coords)
     end
 end
 
-function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, gallery)
-    # Main implementation goes here
+function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
     ncompx = dim
     prob_size < 0 && (prob_size = 256*1024)
 
@@ -39,15 +38,24 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
     sol_basis =
         create_tensor_h1_lagrange_basis(ceed, dim, 1, sol_order + 1, num_qpts, GAUSS)
:...skipping...
diff --git a/julia/LibCEED.jl/examples/ex3-volume.jl b/julia/LibCEED.jl/examples/ex3-volume.jl
index e1f11462..68edf598 100644
--- a/julia/LibCEED.jl/examples/ex3-volume.jl
+++ b/julia/LibCEED.jl/examples/ex3-volume.jl
@@ -28,8 +28,7 @@ function transform_mesh_coords!(dim, mesh_size, mesh_coords)
     end
 end
 
-function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, gallery)
-    # Main implementation goes here
+function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
     ncompx = dim
     prob_size < 0 && (prob_size = 256*1024)
 
@@ -39,15 +38,24 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
     sol_basis =
         create_tensor_h1_lagrange_basis(ceed, dim, 1, sol_order + 1, num_qpts, GAUSS)
 
+    # Determine the mesh size based on the given approximate problem size.
     nxyz = get_cartesian_mesh_size(dim, sol_order, prob_size)
-    println("Mesh size:", nxyz)
+    println("Mesh size: ", nxyz)
 
-    #Build CeedElemRestriction objects describing the mesh and solution discrete
-    # mesh_rstr: for building (no qdata restriction needed)
+    # Build CeedElemRestriction objects describing the mesh and solution discrete
+    # representations.
     mesh_size, mesh_rstr, _ =
         build_cartesian_restriction(ceed, dim, nxyz, mesh_order, ncompx, num_qpts)
-
-    # sol_rstr + sol_rstr_i: for solving
+    num_q_comp = 1 + div(dim*(dim + 1), 2)
+    sol_size, _, qdata_rstr_i = build_cartesian_restriction(
+        ceed,
+        dim,
+        nxyz,
+        sol_order,
+        num_q_comp,
+        num_qpts,
+        mode=StridedOnly,
+    )
     sol_size, sol_rstr, sol_rstr_i = build_cartesian_restriction(
         ceed,
         dim,
@@ -57,15 +65,16 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
         num_qpts,
         mode=RestrictionAndStrided,
     )
+    println("Number of mesh nodes     : ", div(mesh_size, dim))
+    println("Number of solution nodes : ", sol_size)
 
-    # mesh_coords
+    # Create a CeedVector with the mesh coordinates.
     mesh_coords = CeedVector(ceed, mesh_size)
     set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords)
+    # Apply a transformation to the mesh.
     exact_vol = transform_mesh_coords!(dim, mesh_size, mesh_coords)
 
-    #Create the Q-function that builds the mass operator ( i.e it computes the quadrature data) and set its context data.
-    num_q_comp = 1 + div(dim*(dim + 1), 2)
-
+    #Create the Q-function that builds the mass+diffusion operator ( i.e it computes the quadrature data) and set its context data.
     @interior_qf build_qfunc = (
         ceed,
         dim=dim,
@@ -90,27 +99,27 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
         end,
     )
 
-    #Create the operator that builds the quadrature data for the mass operator
+    # Create the operator that builds the quadrature data for the mass+diffusion operator.
     build_oper = Operator(
         ceed,
         qf=build_qfunc,
         fields=[
             (:dx, mesh_rstr, mesh_basis, CeedVectorActive()),
             (:weights, ElemRestrictionNone(), mesh_basis, CeedVectorNone()),
-            (:qdata, sol_rstr_i, BasisNone(), CeedVectorActive()),
+            (:qdata, qdata_rstr_i, BasisNone(), CeedVectorActive()),
         ],
     )
 
-    # Apply to get qdata
+    # Compute the quadrature data for the mass+diff operator.
     elem_qpts = num_qpts^dim
     num_elem = prod(nxyz)
     qdata = CeedVector(ceed, num_elem*elem_qpts*num_q_comp)
-    print("Computing the quadrature data for the mass operator ...")
+    print("Computing the quadrature data for the mass+diffusion operator ...")
     flush(stdout)
     apply!(build_oper, mesh_coords, qdata)
     println(" done.")
 
-    #Create QFunction for applying the mass+diffusion operator
+    # Create the Q-function that defines the action of the mass+diffusion operator.
     @interior_qf apply_qfunc = (
         ceed,
         dim=dim,
@@ -144,32 +153,34 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size, g
             end
         end,
     )
-    apply_oper = Operator(
+
+    # Create the mass+diffusion operator.
+    oper = Operator(
         ceed,
         qf=apply_qfunc,
         fields=[
             (:u, sol_rstr, sol_basis, CeedVectorActive()),
             (:du, sol_rstr, sol_basis, CeedVectorActive()),
-            (:qdata, sol_rstr_i, BasisNone(), qdata),
+            (:qdata, qdata_rstr_i, BasisNone(), qdata),
             (:v, sol_rstr, sol_basis, CeedVectorActive()),
             (:dv, sol_rstr, sol_basis, CeedVectorActive()),
         ],
     )
 
-    # # Compute the mesh volume using the massdiff operator
+    # Compute the mesh volume using the mass+diffusion operator: vol = 1^T \cdot (M + K) \cdot 1
     print("Computing the mesh volume using the formula: vol = 1^T * (M + K) * 1...")
     flush(stdout)
-
+    # Create auxiliary solution-size vectors.
     u = CeedVector(ceed, sol_size)
     v = CeedVector(ceed, sol_size)
+    # Initialize 'u' with ones.
     u[] = 1.0
-
-    # Apply operator
-    apply!(apply_oper, u, v)
-
-    # Compute volume
+    # Apply the mass+diffusion operator: 'u' -> 'v'.
+    apply!(oper, u, v)
+    # Compute and print the sum of the entries of 'v' giving the mesh volume.
     vol = witharray_read(sum, v, MEM_HOST)
 
+    println(" done.")
     @printf("Exact mesh volume    : % .14g\n", exact_vol)
     @printf("Computed mesh volume : % .14g\n", vol)
     @printf("Volume error         : % .14g\n", vol - exact_vol)
@@ -178,10 +189,9 @@ end
 # Entry point
 run_ex3(
     ceed_spec="/cpu/self",
-    dim=2,
-    mesh_order=2,
-    sol_order=2,
-    num_qpts=3,
+    dim=3,
+    mesh_order=4,
+    sol_order=4,
+    num_qpts=4 + 2,
     prob_size=-1,
-    gallery=false,
 )

[jeremylt]

Ope, the CI job lists a couple of more formatting changes to make

[jeremylt]

That last commit looks like it is using a different formatting convention that CI. Did you use the formatting command I gave above?

[jeremylt]

Yeah, the CI job failed because you didn't use the formatting command that CI uses

[KartikiP]

I finally figured out how to run the formatter CI uses. I got this as output

diff --git a/julia/LibCEED.jl/examples/common.jl b/julia/LibCEED.jl/examples/common.jl
index 0e363639..1a93166f 100644
--- a/julia/LibCEED.jl/examples/common.jl
+++ b/julia/LibCEED.jl/examples/common.jl
@@ -112,7 +112,7 @@ function set_cartesian_mesh_coords!(dim, nxyz, mesh_order, mesh_coords)
     nodes = 0.5 .+ 0.5*nodes
 
     @witharray coords = mesh_coords begin
-        @inbounds @simd for gsnodes = 0:scalar_size-1
+        @inbounds @simd for gsnodes = 0:(scalar_size-1)
             rnodes = gsnodes
             for d = 1:dim
                 ndd = nd[d]
diff --git a/julia/LibCEED.jl/examples/ex3-volume.jl b/julia/LibCEED.jl/examples/ex3-volume.jl
index a5df6776..68edf598 100644
--- a/julia/LibCEED.jl/examples/ex3-volume.jl
+++ b/julia/LibCEED.jl/examples/ex3-volume.jl
@@ -54,7 +54,7 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
         sol_order,
         num_q_comp,
         num_qpts,
-        mode = StridedOnly,
+        mode=StridedOnly,
     )
     sol_size, sol_rstr, sol_rstr_i = build_cartesian_restriction(
         ceed,
@@ -63,7 +63,7 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
         sol_order,
         1,
         num_qpts,
-        mode = RestrictionAndStrided,
+        mode=RestrictionAndStrided,
     )
     println("Number of mesh nodes     : ", div(mesh_size, dim))
     println("Number of solution nodes : ", sol_size)
@@ -77,7 +77,7 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
     #Create the Q-function that builds the mass+diffusion operator ( i.e it computes the quadrature data) and set its context data.
     @interior_qf build_qfunc = (
         ceed,
-        dim = dim,
+        dim=dim,
         (dx, :in, EVAL_GRAD, dim, dim),      # ← THIS LINE: dx input
         (weights, :in, EVAL_WEIGHT),         # ← weights input
         (qdata, :out, EVAL_NONE, num_q_comp), # ← qdata output
@@ -102,8 +102,8 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
     # Create the operator that builds the quadrature data for the mass+diffusion operator.
     build_oper = Operator(
         ceed,
-        qf = build_qfunc,
-        fields = [
+        qf=build_qfunc,
+        fields=[
             (:dx, mesh_rstr, mesh_basis, CeedVectorActive()),
             (:weights, ElemRestrictionNone(), mesh_basis, CeedVectorNone()),
             (:qdata, qdata_rstr_i, BasisNone(), CeedVectorActive()),
@@ -122,7 +122,7 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
     # Create the Q-function that defines the action of the mass+diffusion operator.
     @interior_qf apply_qfunc = (
         ceed,
-        dim = dim,
+        dim=dim,
         (u, :in, EVAL_INTERP),
         (du, :in, EVAL_GRAD, dim),
         (qdata, :in, EVAL_NONE, num_q_comp),
@@ -130,7 +130,7 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
         (dv, :out, EVAL_GRAD, dim),
         begin
             # Apply mass: v = qdata[1] * u
-            v .= qdata[1] .* u
+            v .= qdata[1].*u
 
             # Apply diffusion: dv = (qdata[2:end]) * du
             # The qdata contains the symmetric diffusion tensor (J^T*J)
@@ -157,8 +157,8 @@ function run_ex3(; ceed_spec, dim, mesh_order, sol_order, num_qpts, prob_size)
     # Create the mass+diffusion operator.
     oper = Operator(
         ceed,
-        qf = apply_qfunc,
-        fields = [
+        qf=apply_qfunc,
+        fields=[
             (:u, sol_rstr, sol_basis, CeedVectorActive()),
             (:du, sol_rstr, sol_basis, CeedVectorActive()),
             (:qdata, qdata_rstr_i, BasisNone(), qdata),
@@ -188,10 +188,10 @@ end
 
 # Entry point
 run_ex3(
-    ceed_spec = "/cpu/self",
-    dim = 3,
-    mesh_order = 4,
-    sol_order = 4,
-    num_qpts = 4 + 2,
-    prob_size = -1,
+    ceed_spec="/cpu/self",
+    dim=3,
+    mesh_order=4,
+    sol_order=4,
+    num_qpts=4 + 2,
+    prob_size=-1,
 )
diff --git a/julia/LibCEED.jl/src/CeedVector.jl b/julia/LibCEED.jl/src/CeedVector.jl
index 6db67ca9..4487bbfa 100644
--- a/julia/LibCEED.jl/src/CeedVector.jl
+++ b/julia/LibCEED.jl/src/CeedVector.jl
@@ -166,7 +166,7 @@ function witharray_parse(assignment, args)
     mtype = MEM_HOST
     sz = :((length($(esc(v))),))
     body = args[end]
-    for i = 1:length(args)-1
+    for i = 1:(length(args)-1)
         a = args[i]
         if !Meta.isexpr(a, :(=))
             error("Incorrect call to @witharray or @witharray_read") # COV_EXCL_LINE
diff --git a/julia/LibCEED.jl/src/Cuda.jl b/julia/LibCEED.jl/src/Cuda.jl
index f6ce1f01..0d87056a 100644
--- a/julia/LibCEED.jl/src/Cuda.jl
+++ b/julia/LibCEED.jl/src/Cuda.jl
@@ -92,7 +92,7 @@ function generate_kernel(qf_name, kf, dims_in, dims_out)
 
     read_quads_in = [
         :(
-            for j = 1:$(input_sz[i])
+            for j = 1:($(input_sz[i]))
                 $(f_ins_j[i]) = unsafe_load(
                     reinterpret($device_ptr_type, fields.inputs[$i]),
                     q + (j - 1)*Q,
@@ -104,7 +104,7 @@ function generate_kernel(qf_name, kf, dims_in, dims_out)
 
     write_quads_out = [
         :(
-            for j = 1:$(output_sz[i])
+            for j = 1:($(output_sz[i]))
                 unsafe_store!(
                     reinterpret($device_ptr_type, fields.outputs[$i]),
                     $(f_outs[i])[j],
diff --git a/julia/LibCEED.jl/src/UserQFunction.jl b/julia/LibCEED.jl/src/UserQFunction.jl
index b057aa6d..ee3eb48d 100644
--- a/julia/LibCEED.jl/src/UserQFunction.jl
+++ b/julia/LibCEED.jl/src/UserQFunction.jl
@@ -78,7 +78,7 @@ function generate_user_qfunction(
                 $(const_assignments...)
                 $ctx_assignment
                 $(arrays...)
-                @inbounds @simd for $idx = 1:$Q
+                @inbounds @simd for $idx = 1:($Q)
                     $(array_views...)
                     $body
                 end
@@ -143,7 +143,7 @@ function meta_user_qfunction(ceed, def_module, qf, args)
     names_in = Symbol[]
     names_out = Symbol[]
 
-    for a ∈ args[1:end-1]
+    for a ∈ args[1:(end-1)]
         if Meta.isexpr(a, :(=))
             a1 = Meta.quot(a.args[1])
             a2 = esc(a.args[2])
diff --git a/julia/LibCEED.jl/test/rundevtests.jl b/julia/LibCEED.jl/test/rundevtests.jl
index 59d0e484..dcc506f6 100644
--- a/julia/LibCEED.jl/test/rundevtests.jl
+++ b/julia/LibCEED.jl/test/rundevtests.jl
@@ -21,7 +21,7 @@ end
         )
         b = create_tensor_h1_lagrange_basis(c, 3, 1, 3, 3, GAUSS_LOBATTO)
         n = getnumnodes(b)
-        offsets = Vector{CeedInt}(0:n-1)
+        offsets = Vector{CeedInt}(0:(n-1))
         r = create_elem_restriction(c, 1, n, 1, 1, n, offsets)
         op = Operator(
             c;
diff --git a/julia/LibCEED.jl/test/runtests.jl b/julia/LibCEED.jl/test/runtests.jl
index 724240d7..dc460117 100644
--- a/julia/LibCEED.jl/test/runtests.jl
+++ b/julia/LibCEED.jl/test/runtests.jl
@@ -229,7 +229,7 @@ else
             )
             b = create_tensor_h1_lagrange_basis(c, 3, 1, 3, 3, GAUSS_LOBATTO)
             n = getnumnodes(b)
-            offsets = Vector{CeedInt}(0:n-1)
+            offsets = Vector{CeedInt}(0:(n-1))
             r = create_elem_restriction(c, 1, n, 1, 1, n, offsets)
             op = Operator(
                 c;
@@ -277,7 +277,7 @@ else
         @testset "ElemRestriction" begin
             c = Ceed()
             n = 10
-            offsets = Vector{CeedInt}([0:n-1; n-1:2*n-2])
+            offsets = Vector{CeedInt}([0:(n-1); (n-1):(2*n-2)])
             lsize = 2*n - 1
             r = create_elem_restriction(c, 2, n, 1, lsize, lsize, offsets)
             @test getcompstride(r) == lsize
@@ -489,7 +489,7 @@ else
             )
 
             lv = Vector{CeedScalar}(undef, nelem + 1)
-            for i = 1:nelem+1
+            for i = 1:(nelem+1)
                 lv[i] = 10 + i - 1
             end

[jeremylt]

Getting close. The two different formatters sometimes have odd interactions, so I'd

1. roll back to 417cd67
2. Apply my last patch
3. Run the formatting command that CI uses
4. commit the result
5. push both (my patch and your formatting) as two commits

With the two, you should see ex3 give low error when you run it manually and CI should pass

[KartikiP]

First of all, Thank you so much for being soo kind and patient :)

Second,
By "my last patch" you mean this, right?

[jeremylt]

Yup, that patch.

And no worries - it can be a lot all at once to make a new contribution to a project, especially in a new language