Pass the library to the pointwise macros

CHANGELOG.md

@@ -6,6 +6,7 @@
 
 - `local_arrays` and `restore_local_arrays!` warn: use `with_local_arrays!`. An exception between the two calls left both Vecs checked out.
 - `local_coordinate_array(dm)` warns: use `with_local_coordinate_array(f, dm)`. It now returns a copy. It used to leave the coordinate vector checked out until the garbage collector freed the array, and if the DM was destroyed first, that late hand-back crashed Julia.
+- `@simple_fn`, `@residual_fn`, `@jacobian_fn` and `@bd_fn` take the library first, `@residual_fn(petsclib, f0, Nf)`, and use its `PetscInt`, `PetscReal` and `PetscScalar` for the C signature. The library value or its type both work; it is evaluated when the macro expands, so it must be a global defined before the call. Without the library, the macros take `PetscInt`, `PetscScalar` and `PetscReal` from the calling module's globals as before, and warn. That form needed those three consts and allowed only one library per module.
 
 ### Added
 

docs/src/man/dmplex.md

@@ -98,8 +98,11 @@ PETSc.set_name!(petsclib, dm, "stokes")
 ## Pointwise Callbacks
 
 PETSc's FEM assembly calls user-supplied C function pointers at each
-quadrature point.  Three macros generate the required `@cfunction` wrappers
-from plain Julia functions.
+quadrature point.  Four macros (`@simple_fn`, `@residual_fn`, `@jacobian_fn`,
+`@bd_fn`) generate the required `@cfunction` wrappers from plain Julia
+functions.  Their first argument is the library (or its type), which sets the
+C types; it is evaluated when the macro expands, so it must be a global
+defined before the macro call.
 
 ### `@simple_fn` — boundary / exact-solution functions
 
@@ -111,7 +114,7 @@ function bg_vel(t, x, u, ctx)
     u[1] = x[1] * 0.1    # Vx = exx * x
     u[2] = x[2] * (-0.1) # Vy = eyy * y
 end
-const bg_vel_ptr = PETSc.@simple_fn(bg_vel)
+const bg_vel_ptr = PETSc.@simple_fn(petsclib, bg_vel)
 ```
 
 ### `@residual_fn` — f0 / f1 residual and post-processing functions
@@ -127,7 +130,7 @@ function f0_p(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     f0[1] = 0.0
     for d in 0:dim_-1; f0[1] -= u_x[d*dim_+d+1]; end
 end
-const f0_p_ptr = PETSc.@residual_fn(f0_p, 1)
+const f0_p_ptr = PETSc.@residual_fn(petsclib, f0_p, 1)
 ```
 
 ### `@jacobian_fn` — g0/g1/g2/g3 Jacobian functions
@@ -140,7 +143,7 @@ function g1_pu(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, utShift, x, nC, cst, g1)
     for d in 0:dim_-1; g1[d*dim_+d+1] = -1.0; end
 end
-const g1_pu_ptr = PETSc.@jacobian_fn(g1_pu, dim_*dim_)
+const g1_pu_ptr = PETSc.@jacobian_fn(petsclib, g1_pu, dim_*dim_)
 ```
 
 !!! note "Automatic differentiation"

examples/dmplex/ex12.jl

@@ -130,8 +130,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -148,13 +146,13 @@ let _b = get(NamedTuple(pairs(opts)), :bc,    nothing); global bc    = _b === no
 function quadratic_u(t, x, u, ctx)
     u[1] = sum(xi^2 for xi in x)
 end
-const quadratic_u_ptr = PETSc.@simple_fn(quadratic_u)
+const quadratic_u_ptr = PETSc.@simple_fn(petsclib, quadratic_u)
 
 # Zero initial guess for the interior DOFs.
 function zero_u(t, x, u, ctx)
     fill!(u, 0)
 end
-const zero_u_ptr = PETSc.@simple_fn(zero_u)
+const zero_u_ptr = PETSc.@simple_fn(petsclib, zero_u)
 
 # ── COEFF_FIELD:  κ(x) = x₂ + 1 ─────────────────────────────────────────────
 #
@@ -166,22 +164,22 @@ function f0_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0)
     f0[1] = (2*dim_+2)*x[2] + 2*dim_
 end
-const f0_field_ptr = PETSc.@residual_fn(f0_field, Nf)
+const f0_field_ptr = PETSc.@residual_fn(petsclib, f0_field, Nf)
 
 function f1_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f1)
     kap = x[2] + 1.0
     @inbounds for d in 1:dim_; f1[d] = kap * u_x[d]; end
 end
-const f1_field_ptr = PETSc.@residual_fn(f1_field, dim_)
+const f1_field_ptr = PETSc.@residual_fn(petsclib, f1_field, dim_)
 
 function g3_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                   aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g3)
     kap = x[2] + 1.0
     fill!(g3, 0)
     @inbounds for d in 1:dim_; g3[(d-1)*dim_ + d] = kap; end
 end
-const g3_field_ptr = PETSc.@jacobian_fn(g3_field, dim_*dim_)
+const g3_field_ptr = PETSc.@jacobian_fn(petsclib, g3_field, dim_*dim_)
 
 # Neumann flux for COEFF_FIELD:  h = κ(x)·∇u·n = (x₂+1)·2(x·n)
 # Convention: f0_bd = −h  (PETSc boundary residual sign)
@@ -192,7 +190,7 @@ function f0_bd_field(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     @inbounds for d in 1:dim_; xdotn += x[d]*n[d]; end
     f0[1] = -kap * 2.0 * xdotn
 end
-const f0_bd_field_ptr = PETSc.@bd_fn(f0_bd_field, Nf)
+const f0_bd_field_ptr = PETSc.@bd_fn(petsclib, f0_bd_field, Nf)
 
 # ── COEFF_NONLINEAR:  κ(u) = u + 1 ──────────────────────────────────────────
 #
@@ -214,28 +212,28 @@ function f0_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     u_val = u[Int(uOff[1])+1]
     f0[1] = (2*dim_+4)*u_val + 2*dim_
 end
-const f0_nl_ptr = PETSc.@residual_fn(f0_nl, Nf)
+const f0_nl_ptr = PETSc.@residual_fn(petsclib, f0_nl, Nf)
 
 function f1_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f1)
     kap = u[Int(uOff[1])+1] + 1.0
     @inbounds for d in 1:dim_; f1[d] = kap * u_x[d]; end
 end
-const f1_nl_ptr = PETSc.@residual_fn(f1_nl, dim_)
+const f1_nl_ptr = PETSc.@residual_fn(petsclib, f1_nl, dim_)
 
 # g0: ∂f0/∂u = 2·dim+4
 function g0_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g0)
     g0[1] = 2*dim_ + 4
 end
-const g0_nl_ptr = PETSc.@jacobian_fn(g0_nl, Nf*Nf)
+const g0_nl_ptr = PETSc.@jacobian_fn(petsclib, g0_nl, Nf*Nf)
 
 # g2: ∂f1[d]/∂u = u_x[d]   (κ′(u)=1, so ∂((u+1)·u_x[d])/∂u = u_x[d])
 function g2_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g2)
     @inbounds for d in 1:dim_; g2[d] = u_x[d]; end
 end
-const g2_nl_ptr = PETSc.@jacobian_fn(g2_nl, Nf*Nf*dim_)
+const g2_nl_ptr = PETSc.@jacobian_fn(petsclib, g2_nl, Nf*Nf*dim_)
 
 # g3: ∂f1[d]/∂(∂u/∂xd′) = (u+1)·δ_{dd′}
 function g3_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -244,7 +242,7 @@ function g3_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     fill!(g3, 0)
     @inbounds for d in 1:dim_; g3[(d-1)*dim_ + d] = kap; end
 end
-const g3_nl_ptr = PETSc.@jacobian_fn(g3_nl, dim_*dim_)
+const g3_nl_ptr = PETSc.@jacobian_fn(petsclib, g3_nl, dim_*dim_)
 
 # Neumann flux for COEFF_NONLINEAR:  h = (u+1)·2(x·n)
 # f0_bd = −h
@@ -255,7 +253,7 @@ function f0_bd_nl(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     @inbounds for d in 1:dim_; xdotn += x[d]*n[d]; end
     f0[1] = -kap * 2.0 * xdotn
 end
-const f0_bd_nl_ptr = PETSc.@bd_fn(f0_bd_nl, Nf)
+const f0_bd_nl_ptr = PETSc.@bd_fn(petsclib, f0_bd_nl, Nf)
 
 # ── Build the mesh ────────────────────────────────────────────────────────────
 # 3-D simplex (TetGen) is broken in the current PETSc_jll; use hexahedra for 3D.

examples/dmplex/ex17.jl

@@ -123,8 +123,6 @@ MPI.Initialized() || MPI.Init()
 PETSc.initialize(petsclib)
 
 const PetscInt    = petsclib.PetscInt
-const PetscScalar = Float64
-const PetscReal   = Float64
 
 comm = MPI.COMM_WORLD
 
@@ -154,39 +152,39 @@ VALID_SOL = ("vlap_quad", "elas_quad", "vlap_trig", "elas_trig",
 function zero_u(t, x, u, ctx)
     fill!(u, 0)
 end
-const zero_u_ptr = PETSc.@simple_fn(zero_u)
+const zero_u_ptr = PETSc.@simple_fn(petsclib, zero_u)
 
 function ge_shift_u(t, x, u, ctx)
     u[1] = 0.1
     for d in 2:length(u); u[d] = 0.0; end
 end
-const ge_shift_u_ptr = PETSc.@simple_fn(ge_shift_u)
+const ge_shift_u_ptr = PETSc.@simple_fn(petsclib, ge_shift_u)
 
 function quadratic_2d_u(t, x, u, ctx)
     u[1] = x[1]^2
     u[2] = x[2]^2 - 2*x[1]*x[2]
 end
-const quadratic_2d_u_ptr = PETSc.@simple_fn(quadratic_2d_u)
+const quadratic_2d_u_ptr = PETSc.@simple_fn(petsclib, quadratic_2d_u)
 
 function quadratic_3d_u(t, x, u, ctx)
     u[1] = x[1]^2
     u[2] = x[2]^2 - 2*x[1]*x[2]
     u[3] = x[3]^2 - 2*x[2]*x[3]
 end
-const quadratic_3d_u_ptr = PETSc.@simple_fn(quadratic_3d_u)
+const quadratic_3d_u_ptr = PETSc.@simple_fn(petsclib, quadratic_3d_u)
 
 function trig_2d_u(t, x, u, ctx)
     u[1] = sin(2π * x[1])
     u[2] = sin(2π * x[2]) - 2*x[1]*x[2]
 end
-const trig_2d_u_ptr = PETSc.@simple_fn(trig_2d_u)
+const trig_2d_u_ptr = PETSc.@simple_fn(petsclib, trig_2d_u)
 
 function trig_3d_u(t, x, u, ctx)
     u[1] = sin(2π * x[1])
     u[2] = sin(2π * x[2]) - 2*x[1]*x[2]
     u[3] = sin(2π * x[3]) - 2*x[2]*x[3]
 end
-const trig_3d_u_ptr = PETSc.@simple_fn(trig_3d_u)
+const trig_3d_u_ptr = PETSc.@simple_fn(petsclib, trig_3d_u)
 
 function axial_disp_u(t, x, u, ctx)
     mu = mu_val; lambda = lambda_val; N = N_val
@@ -195,15 +193,15 @@ function axial_disp_u(t, x, u, ctx)
     u[2] = 0.25*lambda / (mu*(lambda+mu)) * N * x[2]
     for d in 3:length(u); u[d] = 0.0; end
 end
-const axial_disp_u_ptr = PETSc.@simple_fn(axial_disp_u)
+const axial_disp_u_ptr = PETSc.@simple_fn(petsclib, axial_disp_u)
 
 function uniform_strain_u(t, x, u, ctx)
     eps_xx = 0.1; eps_xy = 0.3; eps_yy = 0.25
     u[1] = eps_xx * x[1] + eps_xy * x[2]
     u[2] = eps_xy * x[1] + eps_yy * x[2]
     for d in 3:length(u); u[d] = 0.0; end
 end
-const uniform_strain_u_ptr = PETSc.@simple_fn(uniform_strain_u)
+const uniform_strain_u_ptr = PETSc.@simple_fn(petsclib, uniform_strain_u)
 
 # ── Residual functions f0 (zeroth-order) ───────────────────────────────────────
 # f0 output size = Nc = dim (one entry per displacement component).
@@ -214,7 +212,7 @@ function f0_vlap_quadratic(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                            aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0)
     for d in 1:dim_; f0[d] += 2.0; end
 end
-const f0_vlap_quadratic_ptr = PETSc.@residual_fn(f0_vlap_quadratic, dim_)
+const f0_vlap_quadratic_ptr = PETSc.@residual_fn(petsclib, f0_vlap_quadratic, dim_)
 
 # Linear elasticity, quadratic solution: body force from −∇·σ.
 function f0_elas_quadratic(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -223,14 +221,14 @@ function f0_elas_quadratic(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     for d in 1:dim_-1; f0[d] += 2.0*mu; end
     f0[dim_] += 2.0*lambda + 4.0*mu
 end
-const f0_elas_quadratic_ptr = PETSc.@residual_fn(f0_elas_quadratic, dim_)
+const f0_elas_quadratic_ptr = PETSc.@residual_fn(petsclib, f0_elas_quadratic, dim_)
 
 # Vector Laplacian, trig solution: body force = 4π² sin(2πx_d) per component.
 function f0_vlap_trig(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                       aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f0)
     for d in 1:dim_; f0[d] += -4.0*π^2 * sin(2π * x[d]); end
 end
-const f0_vlap_trig_ptr = PETSc.@residual_fn(f0_vlap_trig, dim_)
+const f0_vlap_trig_ptr = PETSc.@residual_fn(petsclib, f0_vlap_trig, dim_)
 
 # Linear elasticity, trig solution: body force from −∇·σ.
 function f0_elas_trig(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
@@ -242,7 +240,7 @@ function f0_elas_trig(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                  (d < dim_ ? 2.0*(mu + lambda) : 0.0)
     end
 end
-const f0_elas_trig_ptr = PETSc.@residual_fn(f0_elas_trig, dim_)
+const f0_elas_trig_ptr = PETSc.@residual_fn(petsclib, f0_elas_trig, dim_)
 
 # ── Residual functions f1 (first-order, diffusion flux) ────────────────────────
 # f1 output size = Nc * dim = dim * dim (gradient of Nc-component field).
@@ -254,7 +252,7 @@ function f1_vlap(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_ + d + 1] += u_x[c*dim_ + d + 1]
     end
 end
-const f1_vlap_ptr = PETSc.@residual_fn(f1_vlap, dim_*dim_)
+const f1_vlap_ptr = PETSc.@residual_fn(petsclib, f1_vlap, dim_*dim_)
 
 function f1_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                  aOff, aOff_x, a, a_t, a_x, t, x, numConstants, constants, f1)
@@ -264,7 +262,7 @@ function f1_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         f1[c*dim_ + c + 1] += lambda * u_x[d*dim_ + d + 1]
     end
 end
-const f1_elas_ptr = PETSc.@residual_fn(f1_elas, dim_*dim_)
+const f1_elas_ptr = PETSc.@residual_fn(petsclib, f1_elas, dim_*dim_)
 
 # ── Jacobian functions ─────────────────────────────────────────────────────────
 # g3 output size = Nc*Nc*dim*dim = dim^4 (full stiffness tensor).
@@ -277,7 +275,7 @@ function g3_vlap(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((c*dim_ + c)*dim_ + d)*dim_ + d + 1] = 1.0
     end
 end
-const g3_vlap_ptr = PETSc.@jacobian_fn(g3_vlap, dim_*dim_*dim_*dim_)
+const g3_vlap_ptr = PETSc.@jacobian_fn(petsclib, g3_vlap, dim_*dim_*dim_*dim_)
 
 function g3_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                  aOff, aOff_x, a, a_t, a_x, t, u_tShift, x, numConstants, constants, g3)
@@ -289,7 +287,7 @@ function g3_elas(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
         g3[((c*dim_ + d)*dim_ + c)*dim_ + d + 1] += lambda
     end
 end
-const g3_elas_ptr = PETSc.@jacobian_fn(g3_elas, dim_*dim_*dim_*dim_)
+const g3_elas_ptr = PETSc.@jacobian_fn(petsclib, g3_elas, dim_*dim_*dim_*dim_)
 
 # ── VTK post-processing: displacement (3-vector) and stress (3×3 tensor) ──────
 #
@@ -307,7 +305,7 @@ function copy_displacement_3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     for c in 1:dim_; out[c] = u[c]; end
     for c in dim_+1:3; out[c] = 0.0; end
 end
-const copy_displacement_3_ptr = PETSc.@residual_fn(copy_displacement_3, 3)
+const copy_displacement_3_ptr = PETSc.@residual_fn(petsclib, copy_displacement_3, 3)
 
 # Compute Cauchy stress σ_{cd} = μ(∂u_c/∂x_d + ∂u_d/∂x_c) + λ tr(ε) δ_{cd}
 # and place it in a 3×3 output (9 components, row-major).
@@ -323,7 +321,7 @@ function compute_stress_3x3(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
     end
     for c in 0:dim_-1; out[c*3+c+1] += lambda * div_u; end
 end
-const compute_stress_3x3_ptr = PETSc.@residual_fn(compute_stress_3x3, 9)
+const compute_stress_3x3_ptr = PETSc.@residual_fn(petsclib, compute_stress_3x3, 9)
 
 # ── Neumann BC for elas_axial_disp: traction N on right wall ──────────────────
 # f0_bd output size = Nc = dim.
@@ -332,7 +330,7 @@ function f0_elas_axial_disp_bd(dim_, Nf, NfAux, uOff, uOff_x, u, u_t, u_x,
                                 numConstants, constants, f0)
     f0[1] = real(constants[3])   # N — tension force, constants[3] = N (1-indexed)
 end
-const f0_elas_axial_disp_bd_ptr = PETSc.@bd_fn(f0_elas_axial_disp_bd, dim_)
+const f0_elas_axial_disp_bd_ptr = PETSc.@bd_fn(petsclib, f0_elas_axial_disp_bd, dim_)
 
 # ── Near-null space constructor for AMG ────────────────────────────────────────
 # Called by PETSc for each level when setting up the preconditioner.