From 91dab136585410ea38c5818f1de8c7dda7ca9ead Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jesus=20Mart=C3=ADnez=20Frutos?= Date: Tue, 10 Nov 2026 17:51:17 +0100 Subject: [PATCH 1/3] =?UTF-8?q?Second=20derivatives=20(=E2=88=87=E2=88=87,?= =?UTF-8?q?=20=CE=94)=20on=20skeleton=20integrals?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit jump(∇∇(u)), mean(Δ(u)), ... of the unknown, test and trial functions are compiled by the facet kernels instead of sending the skeleton to the Gridap fallback. The Hessian of each side uses the Jt⁻¹ and the curvature of that side's cell map at its mapped points (side_geometry; hess_use of a skeleton class), in the values of the unknown (field_values_sides with one Φ per side and field), the physical weights of test/trial functions and the pull-back of the residual coefficients, so ⟦∇∇u⟧ is exact on non-affine quadrilaterals and hexahedra too. The matrix-free data of a skeleton class stores the curvature of each side (only with Hessians on non-affine cells); weight buffers sized with wcols. Tests (test_hessian.jl, replacing the fallback test): C0 interior penalty with {Δu}⟦∂ₙv⟧ and a nonlinear Hessian-jump term, and the third-medium regularisation of the gear problem (σβ/h_F ⟦∇u⟧:⟦∇v⟧ + c₂βh_F ⟦∇∇u⟧:⟦∇∇v⟧ on the skeleton of a subdomain) against Gridap on P2/P3 triangles, Q2 quadrilaterals and P2 tetrahedra, against hessian_nonaffine on non-affine quadrilaterals and hexahedra; zero Hessian jump of a quadratic field; multi-field, Hessian only in the Jacobian, update, matrix-free rebind; zero allocations, threads. --- CHANGELOG.md | 8 ++ README.md | 9 ++- src/kernels.jl | 11 ++- src/matrixfree.jl | 24 ++++-- src/skeleton.jl | 52 +++++++++---- test/README.md | 2 +- test/test_hessian.jl | 174 ++++++++++++++++++++++++++++++++++++++++--- 7 files changed, 241 insertions(+), 39 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 595c28b..c6e913a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,13 @@ # Changelog +## Unreleased + +- Second derivatives (`∇∇`, `Δ`) on skeleton integrals: `jump(∇∇(u))`, `mean(∇∇(u))`, `mean(Δ(u))`, … of + the unknown, test and trial functions are compiled by the facet kernels (C⁰ interior-penalty terms, + Hessian-jump regularisations ⟦∇∇u⟧:⟦∇∇v⟧), with the Hessian push-forward of each side (its Jt⁻¹ and the + curvature of its cell map: exact on non-affine quadrilaterals and hexahedra); residual, assembled and + matrix-free Jacobian. They no longer send the skeleton to the Gridap fallback. + ## 0.2.2 - p-multigrid preconditioner `PMultigrid` for `MatrixFreeSolver`: nested Lagrangian levels on the same mesh diff --git a/README.md b/README.md index 705ff4b..e3f3735 100644 --- a/README.md +++ b/README.md @@ -31,7 +31,8 @@ backend for its drivers (`backend = :compiled`). (Hessian regularisations `β ∫ ∇∇u ⊙ ∇∇v`, biharmonic and C⁰-interior-penalty terms) with Gridap's layout (`TensorValue{D,D}`, `ThirdOrderTensorValue{D,D,N}` with `H[a,b,k] = ∂²u_k/∂x_a∂x_b`), on simplices and Cartesian cells as Gridap's `∇∇` and, unlike it, also on non-affine - quadrilaterals and hexahedra (the curvature of the cell map is included). + quadrilaterals and hexahedra (the curvature of the cell map is included). Also on skeleton + integrals: `jump(∇∇(u))`, `mean(Δ(u))`, … with the Hessian of each side. - **Spaces:** Lagrangian spaces of any order and H(div)-conforming Raviart–Thomas spaces (contravariant Piola map, per-cell dof signs) on hexahedra, quadrilaterals, tetrahedra and triangles; single- and multi-field (e.g. Biot u–W, mixed Darcy); Dirichlet masks; subdomains. Lagrangian spaces also on @@ -318,9 +319,9 @@ run with both backends as part of the test suite. - On manifold triangulations (cells of lower dimension than the space): Lagrangian spaces in cell integrals. Boundary and skeleton integrals of those spaces and Raviart–Thomas fields on them go to the Gridap fallback; internal variables on them are not supported yet. -- Second derivatives of Lagrangian fields in cell and boundary integrals and in the matrix-free - Jacobian (dense weights; sum factorisation is used for the fields without Hessian terms). `∇∇` on - skeleton integrals and of Raviart–Thomas fields goes to the Gridap fallback; third derivatives too. +- Second derivatives of Lagrangian fields in cell, boundary and skeleton integrals and in the matrix-free + Jacobian (dense weights; sum factorisation is used for the fields without Hessian terms). `∇∇` of + Raviart–Thomas fields goes to the Gridap fallback; third derivatives too. - p-multigrid for single-field problems; multi-field (u, φ) saddle points need a block smoother and are not covered yet (use `LaggedAssembled` with HyperFEM's block preconditioners). - Shared-memory parallelism (threads). Distributed runs with GridapDistributed and GPU execution diff --git a/src/kernels.jl b/src/kernels.jl index c206b7a..b0515b5 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -263,9 +263,14 @@ node loops of every field are compiled in the caller's loop. `Φ` is the curvatu :(($(ex...),)) end -"Same with one `iJt` per field (skeleton: the inverse Jacobian of each side)." -@inline @generated function field_values_sides(tables::NTuple{NF,Any}, locs, use, xl, q, iJs) where {NF} - ex = [:(field_point_values(tables[$i], xl, locs[$i], q, iJs[$i], use[$i])) for i in 1:NF] +""" +Same with one `iJt` and one curvature `Φ` per field (skeleton: the inverse Jacobian and the curvature of +the cell map of each side; `Φs` a tuple of `nothing` when no field uses Hessians). +""" +@inline field_values_sides(tables::NTuple{NF,Any}, locs, use, xl, q, iJs) where {NF} = + field_values_sides(tables, locs, use, xl, q, iJs, ntuple(_ -> nothing, Val(NF))) +@inline @generated function field_values_sides(tables::NTuple{NF,Any}, locs, use, xl, q, iJs, Φs) where {NF} + ex = [:(field_point_values(tables[$i], xl, locs[$i], q, iJs[$i], nothing, Φs[$i], use[$i])) for i in 1:NF] :(($(ex...),)) end diff --git a/src/matrixfree.jl b/src/matrixfree.jl index 5a1f007..5be634e 100644 --- a/src/matrixfree.jl +++ b/src/matrixfree.jl @@ -8,7 +8,8 @@ # H(div) fields use the dense path with the stored Piola geometry (Jt and Jt⁻¹) and the dof signs # at gather/scatter. # Skeleton groups do the same per facet with the increments of both cells (dense weights, the -# Jt⁻¹ of each side; the T and dS arrays are shared by the two `GroupData` of a class). +# Jt⁻¹ of each side and, with Hessians on non-affine cells, its curvature; the T and dS arrays are +# shared by the two `GroupData` of a class). struct GroupData{D} T::Array{Float64,3} # nsv*nsu × np × ncells (column-major: T[sv + (su-1)*nsv, q, cell]) @@ -52,10 +53,12 @@ function MatrixFreeJacobian(op::KernelOperator) D = group_dim(g) nsv, nsu = nslots(g) if g isa SkeletonGroup - # one tangent and measure per facet of the class; Jt⁻¹ of the plus and of the minus cell + # one tangent and measure per facet of the class; Jt⁻¹ (and the curvature) of the plus and of the minus cell n = length(g.facets) T, dV = zeros(nsv * nsu, np, n), zeros(np, n) - (GroupData{D}(T, zeros(D * D, np, n), dV), GroupData{D}(T, zeros(D * D, np, n), dV)) + curved = hess_use(g) isa Val{true} && !(g.geo.affine && g.geom.affine) + side() = GroupData{D}(T, zeros(D * D, np, n), zeros(D * D, np, n), dV, zeros(curved ? D * D * D : 0, np, n)) + (side(), side()) else nc = length(g.cellmap) curved = hess_use(g) isa Val{true} && !g.geo.affine @@ -162,15 +165,16 @@ function prepare_cells!(gds::Tuple{GroupData{D},GroupData{D}}, g::SkeletonGroup{ ndl = op.ndl xl = zeros(2ndl) zsv = zero(SVector{NSV,Float64}); zsu = zero(SVector{NSU,Float64}) + hu = hess_use(g) + curved = size(gp.Phi, 1) > 0 for lf in lfacets f = g.facets[lf] gather_cell!(xl, op.fields, x, dvals, g.skel.cells[1, f]) gather_cell!(xl, op.fields, x, dvals, g.skel.cells[2, f], ndl) for q in 1:g.geo.np - dS, iJp, Jp = point_geometry(g.geo, q, lf) - _, iJm, Jm = point_geometry(g.geom, q, lf) - iJs = side_tuple(iJp, iJm, Val(NF2)) - FV = field_values_sides(tables, locs, g.use, xl, q, iJs) + dS, iJs, Φs, (Jp, Jm), (Φp, Φm) = side_geometry(g, q, lf, hu) + iJp, iJm = iJs[1], iJs[end] + FV = field_values_sides(tables, locs, g.use, xl, q, iJs, Φs) Aq = map(a -> @inbounds(a[q, lf]), auxd) T = if LF === Nothing ForwardDiff.jacobian(r -> ForwardDiff.gradient(s -> pf.core(pf.ops, pf.coefs, FV, Aq, s, r), zsv), zsu) @@ -186,6 +190,10 @@ function prepare_cells!(gds::Tuple{GroupData{D},GroupData{D}}, g::SkeletonGroup{ gp.Jt[i, q, lf] = Jp[i] gm.Jt[i, q, lf] = Jm[i] end + if curved + store_curvature!(gp, Φp, q, lf) + store_curvature!(gm, Φm, q, lf) + end gp.dV[q, lf] = dS end end @@ -275,7 +283,7 @@ function apply_buffers(g::SkeletonGroup{NSV,NSU,PF,TT,AUX,SV,SU,D}, tps, ndl::In np = g.geo.np (xl=zeros(2ndl), re=zeros(2ndl), su=zeros(NSU, np), sv=zeros(NSV, np), - Wb=map(t -> zeros(t.nn, D + 1), g.tables), + Wb=map(t -> zeros(t.nn, wcols(t)), g.tables), uc=map(t -> zeros(t.nn), g.tables), rc=map(t -> zeros(t.nn), g.tables), pts=map(t -> zeros(np, D + 1), g.tables), coef=map(t -> zeros(np, D + 1), g.tables), tmp=map(_ -> zeros(0, 0), g.tables)) diff --git a/src/skeleton.jl b/src/skeleton.jl index 0036eb1..6bb0cfc 100644 --- a/src/skeleton.jl +++ b/src/skeleton.jl @@ -11,6 +11,9 @@ # facet (`SkeletonData.nzpos`, the cross blocks are added to the sparsity pattern). The pointwise # integrand sees the unknown of both sides (`FV = (plus fields..., minus fields...)`) and test/trial # slots per side (`Leaf.side`); the measure is that of the plus face (`w · meas(Jf⁺ ⋅ Jt⁺)`). +# Second derivatives (∇∇, Δ) use the Hessian push-forward of the volume kernels on each side, with the +# Jt⁻¹ and the curvature Φ of that side's cell map at its mapped points (`point_curvature` of `geo` and +# `geom`; zero on affine cells), e.g. the C⁰ interior-penalty terms ⟦∇∇u⟧:⟦∇∇v⟧ or {Δu}⟦∇v⋅n⟧. "FE-space cell of the plus and minus side of each facet (2 × nfacets)." function skeleton_cells(trian::Triangulation, space_trian::Triangulation) @@ -122,7 +125,6 @@ function SkeletonGroup(terms::Vector{Term}, trian::Triangulation, Uspaces, locs, ls = vcat(pf.test.leaves, pf.trial.leaves, leaves(terms, :unknown)) isempty(ls) && error("HyperFEMKernels: integral without unknown, test or trial functions") any(l -> l.side === :cell, ls) && error("HyperFEMKernels: a function is used on the skeleton without a side (use plus/minus, jump or mean)") - any(l -> l.deriv == 2, ls) && error("HyperFEMKernels: second derivatives (∇∇) on skeleton integrals are not supported yet") w = uniform_value(terms[1].weights, "the quadrature weights") k = findfirst(g -> g isa SkeletonGroup && g.trian === trian && g.key == cls.key && g.geo.w == w && g.facets == cls.facets, previous) if k === nothing @@ -156,6 +158,32 @@ function SkeletonGroup(terms::Vector{Term}, trian::Triangulation, Uspaces, locs, objectid(terms[1].weights)) end +""" + hess_use(g::SkeletonGroup) -> Val{true/false} + +Whether any (side, field) of the class uses Hessians: the curvature of both cell maps is then +evaluated per point (see `hess_use(::Group)`). +""" +@generated function hess_use(::SkeletonGroup{NSV,NSU,PF,TT,AUX,SV,SU,D,NF2,LF,JF,UU}) where {NSV,NSU,PF,TT,AUX,SV,SU,D,NF2,LF,JF,UU} + hs(T) = any(t -> t.parameters[5] === true, T.parameters) # SlotSpec{NC,LV,HV,HG,HH,...} + hu(T) = any(t -> t.parameters[3] === true, T.parameters) # FieldUse{HV,HG,HH} + :(Val{$(hs(SV) || hs(SU) || hu(UU))}()) +end + +""" + side_geometry(g, q, lf, hu) -> (dS, iJs, Φs) + +Geometry of point `q` of facet `lf` of a class: the facet measure (plus face), and per (side, field) +the inverse Jacobian and the curvature (`nothing` without Hessians) of that side's cell map. +""" +@inline function side_geometry(g::SkeletonGroup{NSV,NSU,PF,TT,AUX,SV,SU,D,NF2}, q::Int, lf::Integer, hu) where {NSV,NSU,PF,TT,AUX,SV,SU,D,NF2} + dS, iJp, Jp = point_geometry(g.geo, q, lf) + _, iJm, Jm = point_geometry(g.geom, q, lf) + Φp = point_curvature(g.geo, q, lf, hu) + Φm = point_curvature(g.geom, q, lf, hu) + (dS, side_tuple(iJp, iJm, Val(NF2)), side_tuple(Φp, Φm, Val(NF2)), (Jp, Jm), (Φp, Φm)) +end + "The skeleton class `g` with its values read from `terms` (see `rebind(::Group, terms)`)." function rebind(g::SkeletonGroup, terms::Vector{Term}) pf = PointForm(terms, g.pf) @@ -214,18 +242,17 @@ function element_residual!(re, g::SkeletonGroup{NSV,NSU,PF,TT,AUX,SV,SU,D,NF2,LF fill!(re, 0.0) tables, locs, pf = g.tables, g.locs, g.pf zsv = zero(SVector{NSV,Float64}) + hu = hess_use(g) for q in 1:g.geo.np - dS, iJp = point_geometry(g.geo, q, lf) - _, iJm = point_geometry(g.geom, q, lf) - iJs = side_tuple(iJp, iJm, Val(NF2)) - FV = field_values_sides(tables, locs, g.use, xl, q, iJs) + dS, iJs, Φs = side_geometry(g, q, lf, hu) + FV = field_values_sides(tables, locs, g.use, xl, q, iJs, Φs) Aq = map(a -> @inbounds(a[q, lf]), auxd) R = if LF === Nothing ForwardDiff.gradient(s -> pf.core(pf.ops, pf.coefs, FV, Aq, s, EMPTY_SLOTS), zsv) else apply_maps(g.lin, g.lin.state(g.lin.ops, FV, Aq), Val(NSV)) end - map((t, l, sm, iJ) -> contract_residual!(re, t, l, sm, R, q, iJ, dS), tables, locs, g.smv, iJs) + map((t, l, sm, iJ, Φ) -> contract_residual!(re, t, l, sm, R, q, iJ, nothing, Φ, dS), tables, locs, g.smv, iJs, Φs) end re end @@ -236,12 +263,11 @@ function element_jacobian!(KeT, g::SkeletonGroup{NSV,NSU,PF,TT,AUX,SV,SU,D,NF2,L tables, locs, pf = g.tables, g.locs, g.pf zsv = zero(SVector{NSV,Float64}) zsu = zero(SVector{NSU,Float64}) + hu = hess_use(g) for q in 1:g.geo.np - dS, iJp = point_geometry(g.geo, q, lf) - _, iJm = point_geometry(g.geom, q, lf) - iJs = side_tuple(iJp, iJm, Val(NF2)) - map((t, W, iJ) -> point_weights!(W, t, q, iJ), tables, Wb, iJs) - FV = field_values_sides(tables, locs, g.use, xl, q, iJs) + dS, iJs, Φs = side_geometry(g, q, lf, hu) + map((t, W, iJ, Φ, sv, su) -> point_weights!(W, t, q, iJ, nothing, Φ, hess_weights(sv, su)), tables, Wb, iJs, Φs, g.smv, g.smu) + FV = field_values_sides(tables, locs, g.use, xl, q, iJs, Φs) Aq = map(a -> @inbounds(a[q, lf]), auxd) T = if LF === Nothing ForwardDiff.jacobian(r -> ForwardDiff.gradient(s -> pf.core(pf.ops, pf.coefs, FV, Aq, s, r), zsv), zsu) # T[test slot, trial slot] @@ -262,7 +288,7 @@ function cells!(::Val{K}, b, A, op::KernelOperator, g::SkeletonGroup{NSV,NSU,PF, ndl = op.ndl n2 = 2ndl xl = zeros(n2) - Wb = map(t -> zeros(t.nn, D + 1), g.tables) + Wb = map(t -> zeros(t.nn, wcols(t)), g.tables) skel, facets = g.skel, g.facets if K === :res re = zeros(n2) @@ -344,7 +370,7 @@ function _element_allocations(kind::Symbol, op::KernelOperator, g::SkeletonGroup f = g.facets[1] gather_cell!(xl, op.fields, x, dvals, g.skel.cells[1, f]) gather_cell!(xl, op.fields, x, dvals, g.skel.cells[2, f], ndl) - Wb = map(t -> zeros(t.nn, D + 1), g.tables) + Wb = map(t -> zeros(t.nn, wcols(t)), g.tables) if kind === :residual re = zeros(2ndl) element_residual!(re, g, xl, auxd, 1) diff --git a/test/README.md b/test/README.md index 8513cc3..d423880 100644 --- a/test/README.md +++ b/test/README.md @@ -17,7 +17,7 @@ tolerance is 1e-12 relative. | `test_boundary.jl` | boundary integrals in the engine: Neumann/Robin on quads and triangles (orders 1–3), whole boundary with several faces per cell (threads, colouring), normal vector and ∇u on Γ, 3D follower pressure load on hexahedra and tetrahedra, multi-field, volume fields on Γ as auxiliary data, matrix-free, `update`, zero allocations; every case against pure Gridap | | `test_hdiv.jl` | H(div) (Raviart–Thomas) fields in the engine: Piola tables against Gridap's basis (values, gradients, dof signs; 2D/3D, n-cubes/simplices, orders 0–1), mixed Darcy RT_k–P_k/Q_k (orders 0–2, non-affine meshes), Biot u–W (Lagrangian u + RT W: Biot coupling through the divergences, mass couplings, deformation-dependent permeability, ∇W term; linear maps on/off, threads, `update`, zero allocations), matrix-product vs loop paths, boundary integrals with w ⋅ n on curved faces, matrix-free Jacobian (dense weights with the Piola geometry and the dof signs), skeleton terms with RT sent to the fallback with a clear reason; every case against pure Gridap | | `test_cfd.jl` | weak form of the aorta CFD example (Oseen–BDF2 Navier–Stokes, P2–P1, SUPG/PSPG/LSIC, elementwise Laplacian of the previous step, backflow and Windkessel Robin terms on the outlet) on a small tetrahedral channel: backward Euler and BDF2, linear system with lifted Dirichlet values (A = J, b = −r(0)), data overwritten in place between assemblies, `update` after the first step; every case against pure Gridap | -| `test_hessian.jl` | second derivatives (∇∇, Δ) of the unknown, test and trial functions: exact Hessians of an interpolated quadratic vector field (Q2 affine and non-affine, P2), scalar and vector forms (St Venant–Kirchhoff + β∫∇∇u⊙∇∇v) against Gridap on simplices P1–P3 and Cartesian Q1–Q2 in 2D/3D (linear maps on/off, matrix-free, threads), non-affine quadrilaterals/hexahedra against `hessian_nonaffine` (HyperFEM's push-forward, reproduced in `common.jl`), mixed/multi-field use, `update` without recompilation, ∇∇ of FE functions read from their nodal values, fallback of ∇∇ on skeletons, zero allocations | +| `test_hessian.jl` | second derivatives (∇∇, Δ) of the unknown, test and trial functions: exact Hessians of an interpolated quadratic vector field (Q2 affine and non-affine, P2), scalar and vector forms (St Venant–Kirchhoff + β∫∇∇u⊙∇∇v) against Gridap on simplices P1–P3 and Cartesian Q1–Q2 in 2D/3D (linear maps on/off, matrix-free, threads), non-affine quadrilaterals/hexahedra against `hessian_nonaffine` (HyperFEM's push-forward, reproduced in `common.jl`), mixed/multi-field use, `update` without recompilation, ∇∇ of FE functions read from their nodal values, ∇∇ on skeleton integrals (C⁰ interior penalty with `mean(Δ)` and a nonlinear Hessian-jump term; the third-medium regularisation σβ/h_F ⟦∇u⟧:⟦∇v⟧ + c₂βh_F ⟦∇∇u⟧:⟦∇∇v⟧ on the skeleton of a subdomain; triangles P2–P3, quadrilaterals, tetrahedra against Gridap, non-affine quadrilaterals and hexahedra against `hessian_nonaffine`; zero Hessian jump of a quadratic field; multi-field, Hessian only in the Jacobian, `update`, matrix-free), zero allocations | | `test_update.jl` | cheap `update`: forms of the same structure rebind the compiled groups (shared cores, tables and auxiliary buffers) with residual and Jacobian equal to a fresh operator at 1e-14 — numbers by value, by `Ref`, FE functions changed in place; volume, boundary and skeleton terms with ∇∇, multi-field, `:assembled` and `:matrixfree`; structure changes (term added, other triangulation, other operation, other constant type) take the full path; linear-map values follow the parameters; `HFK_CHECK_REBIND`; allocations independent of the mesh size; `KernelFunctional` | | `test_skeleton.jl` | skeleton integrals in the engine: interior penalty on quads/triangles (orders 1–3), sparsity pattern with the cross blocks and facet positions, facet colouring and threads, `mean`/`jump` of vector fields with the normal on hexahedra and tetrahedra, the weighted jump of ∇u of third-medium contact (nonlinear weight on both sides, facet size as data, non-symmetric tangent, linear maps on/off), skeleton of a subdomain with the space on a sub-triangulation, multi-field, `update` (parameters, reuse of tables and positions, pattern recomputed when a skeleton appears), skeleton term only in the residual, matrix-free Jacobian with facets (scalar/vector, 2D/3D, multi-field, threads), errors; every case against pure Gridap, zero allocations | | `test_subdomains.jl` | configuration of the hMRE + vacuum example: FE spaces on a sub-triangulation with integrals on its parts and different laws, the boundary of a subdomain as internal interface, an auxiliary FE function living on that interface, Dirichlet data on the interface changed in place; every case against pure Gridap | diff --git a/test/test_hessian.jl b/test/test_hessian.jl index 58fac9a..ccbb612 100644 --- a/test/test_hessian.jl +++ b/test/test_hessian.jl @@ -3,7 +3,7 @@ # defined) at 1e-13 relative, and against the push-forward of `hessian_nonaffine` (HyperFEM's # `push_hessian`, reproduced in common.jl) on non-affine quadrilaterals and hexahedra; exact # Hessians of an interpolated quadratic field; FE functions (not the unknown) with ∇∇; the -# matrix-free Jacobian; the fallback of ∇∇ on skeletons; zero allocations. +# matrix-free Jacobian; ∇∇ on skeleton integrals (C⁰ interior penalty, Hessian jumps); zero allocations. Random.seed!(3) const HTOL = 1e-13 @@ -221,17 +221,171 @@ end end end -@testset "∇∇ on a skeleton: clean fallback of that triangulation only" begin +# --- ∇∇ on skeleton integrals -------------------------------------------------------------- +# The Hessian of each side uses the Jt⁻¹ and the curvature of that side's cell map at its mapped points. +# Reference forms: Gridap's ∇∇/Δ on affine cells, `hessian_nonaffine` (and its trace) on mapped n-cubes. + +sk_noalloc(op, x) = all(a -> a.bytes == 0, allocations(op, x)) +Hn_tr(a) = tr(hessian_nonaffine(a)) + +""" +C⁰ interior penalty (biharmonic-like) scalar problem on the skeleton of the whole mesh: jump of the +normal derivative, consistency terms {Δ}⟦∂ₙ⟧, and a nonlinear Hessian-jump term +(1 + {Δu}²) ⟦∇∇u⟧:⟦∇∇v⟧ (dual numbers). `H`, `L` are the Hessian and Laplacian used in the forms +(∇∇ and Δ for the kernels; `hessian_nonaffine` and its trace for the reference on mapped cells). +""" +function c0ip_problem(model; order=2, degree=2order, H=∇∇, L=Δ) + Ω = Triangulation(model); dΩ = Measure(Ω, degree) + Λ = SkeletonTriangulation(model); dΛ = Measure(Λ, degree) + n = get_normal_vector(Λ) + V = TestFESpace(model, ReferenceFE(lagrangian, Float64, order), dirichlet_tags="boundary") + U = TrialFESpace(V, x -> x[1]^2 - 0.5 * x[2]) + hF = CellField(get_array(∫(1)dΛ), Λ) + γ = 3.0 + dn(a) = jump(∇(a) ⋅ n) + res(u, v) = ∫(H(v) ⊙ H(u) + v * (u * u))dΩ + + ∫((γ / hF) * (dn(v) * dn(u)) - mean(L(u)) * dn(v) - mean(L(v)) * dn(u) + + (0.2 * hF) * ((jump(H(u)) ⊙ jump(H(v))) * (1.0 + mean(L(u)) * mean(L(u)))))dΛ + jac(u, du, v) = ∫(H(v) ⊙ H(du) + 2.0 * (v * (u * du)))dΩ + + ∫((γ / hF) * (dn(v) * dn(du)) - mean(L(du)) * dn(v) - mean(L(v)) * dn(du) + + (0.2 * hF) * ((jump(H(du)) ⊙ jump(H(v))) * (1.0 + mean(L(u)) * mean(L(u))) + + (jump(H(u)) ⊙ jump(H(v))) * (2.0 * mean(L(u)) * mean(L(du)))))dΛ + (; model, Ω, dΩ, Λ, dΛ, U, V, res, jac) +end + +""" +Vector problem as in the third-medium regularisation of TMC_gears (G1): SVK on the whole mesh, β ∫∇∇u⊙∇∇v on +a subdomain Ωc and, on the skeleton of Ωc (a composite triangulation), σβ/h_F ⟦∇u⟧:⟦∇v⟧ + c₂βh_F ⟦∇∇u⟧:⟦∇∇v⟧. +""" +function g1_problem(model, D; order=2, H=∇∇) + Ω = Triangulation(model); dΩ = Measure(Ω, 2order) + cc = [sum(c) / length(c) for c in get_cell_coordinates(model)] + Ωc = Triangulation(model, [x[1] > 0.4 for x in cc]); dΩc = Measure(Ωc, 2order) + Λ = SkeletonTriangulation(Ωc); dΛ = Measure(Λ, 2order) + hF = CellField(get_array(∫(1)dΛ), Λ) + β = 5e-3 + cJ, cJ2 = β / hF, (2.0 * β) * hF + V = TestFESpace(model, ReferenceFE(lagrangian, VectorValue{D,Float64}, order), dirichlet_tags="boundary") + U = TrialFESpace(V, x -> VectorValue(ntuple(i -> 0.02 * x[i] * x[mod1(i + 1, D)], D))) + qJ(a, b) = jump(∇(a)) ⊙ jump(∇(b)) + qJ2(a, b) = jump(H(a)) ⊙ jump(H(b)) + res(u, v) = ∫(∇(v)' ⊙ (P_svk ∘ (Fk ∘ (∇(u)'))))dΩ + β * ∫(H(u) ⊙ H(v))dΩc + ∫(cJ * qJ(u, v) + cJ2 * qJ2(u, v))dΛ + jac(u, du, v) = ∫(∇(v)' ⊙ (dP_svk ∘ (Fk ∘ (∇(u)'), ∇(du)')))dΩ + β * ∫(H(du) ⊙ H(v))dΩc + ∫(cJ * qJ(du, v) + cJ2 * qJ2(du, v))dΛ + (; model, Ω, dΩ, Ωc, Λ, dΛ, U, V, res, jac) +end + +@testset "∇∇ on skeleton integrals: C⁰ interior penalty and G1 regularisation against Gridap (affine cells)" begin + # 3D: the scalar problem on tetrahedra here; the vector one (54 slots per side pair) on the hexahedra of the next testset + for (D, simplex, order) in ((2, true, 2), (2, true, 3), (2, false, 2), (3, true, 2)) + model = cartesian_model(D, D == 2 ? (4, 3) : (2, 2, 2); simplex) + probs = D == 2 ? ((c0ip_problem(model; order), "C0-IP"), (g1_problem(model, D; order), "G1")) : ((c0ip_problem(model; order), "C0-IP"),) + for (prob, name) in probs + for linear_maps in (true, false) + c = compare(prob, prob.res, prob.jac; linear_maps) + @test c.res < HTOL && c.jac < HTOL && c.mf < 1e-11 + @test any(g -> g isa HK.SkeletonGroup && HK.hess_use(g) isa Val{true}, c.op.jgroups) + end + op = KernelOperator(prob.res, prob.jac, prob.U, prob.V) + x = randx(prob.U) + @test sk_noalloc(op, x) + op1 = KernelOperator(prob.res, prob.jac, prob.U, prob.V; nthreads=1) + @test reldiff(Algebra.residual(op1, x), Algebra.residual(op, x)) < HTOL + @test reldiff(Algebra.jacobian(op1, x), Algebra.jacobian(op, x)) < HTOL + # simplices: no curvature stored for the matrix-free skeleton groups (the Q1 reference Hessians of + # quadrilaterals do not vanish, so Cartesian cells store a curvature that is zero at run time) + mf = MatrixFreeJacobian(op) + for (g, gd) in zip(op.jgroups, mf.data) + g isa HK.SkeletonGroup && @test (size(gd[1].Phi, 1) == 0 && size(gd[2].Phi, 1) == 0) == simplex + end + end + end +end + +@testset "∇∇ on skeleton integrals: non-affine quadrilaterals and hexahedra against hessian_nonaffine" begin + for D in (2, 3) + model = cartesian_model(D, D == 2 ? (3, 3) : (2, 2, 2); mapped=true) + # 3D: Gridap takes the faces of a mapped Cartesian model as flat (normals, see test_hdiv.jl); not those of an unstructured one + D == 3 && (model = Gridap.Geometry.UnstructuredDiscreteModel(model)) + # 3D: the scalar problem only (the vector one, 54 × 54 slots per facet, takes minutes to compile; its + # curvature path is the same, checked in 2D and by the scalar problem here) + probs = ((c0ip_problem(model), c0ip_problem(model; H=hessian_nonaffine, L=Hn_tr)),) + D == 2 && (probs = (probs..., (g1_problem(model, D), g1_problem(model, D; H=hessian_nonaffine)))) + for (kprob, gprob) in probs + for linear_maps in (true, false) + c = compare(kprob, gprob.res, gprob.jac; linear_maps) + @test c.res < 1e-12 && c.jac < 1e-12 && c.mf < 1e-11 + end + op = KernelOperator(kprob.res, kprob.jac, kprob.U, kprob.V) + @test sk_noalloc(op, randx(kprob.U)) + sk = [g for g in op.jgroups if g isa HK.SkeletonGroup] + @test !isempty(sk) && all(g -> !g.geo.affine && !g.geom.affine, sk) + mf = MatrixFreeJacobian(op) + for (g, gd) in zip(op.jgroups, mf.data) + g isa HK.SkeletonGroup && @test size(gd[1].Phi, 1) == D^3 && size(gd[2].Phi, 1) == D^3 + end + end + end +end + +@testset "∇∇ on skeleton integrals: a field with a constant Hessian has no Hessian jump (Q2 non-affine, P2)" begin + # uq is quadratic: interpolated exactly, its Hessian is the same constant on both sides of every facet + for (mapped, simplex) in ((false, false), (true, false), (false, true)) + model = cartesian_model(2, (4, 3); mapped, simplex) + Λ = SkeletonTriangulation(model); dΛ = Measure(Λ, 4) + V = TestFESpace(model, ReferenceFE(lagrangian, VectorValue{2,Float64}, 2)) + res(u, v) = ∫(jump(∇∇(u)) ⊙ jump(∇∇(v)) + (mean(∇∇(u)) ⊙ mean(∇∇(u)) - Hq ⊙ Hq) * mean(tr(∇(v))))dΛ + jac(u, du, v) = ∫(jump(∇∇(du)) ⊙ jump(∇∇(v)) + (2.0 * (mean(∇∇(u)) ⊙ mean(∇∇(du)))) * mean(tr(∇(v))))dΛ + op = KernelOperator(res, jac, V, V) + @test isempty(op.fallback.rtrians) && isempty(op.fallback.jtrians) + x = get_free_dof_values(interpolate(uq, V)) + b = Algebra.residual(op, x) + @test norm(b, Inf) < 1e-10 + @test norm(Algebra.residual(op, x .+ 1e-3 .* randn(length(x))), Inf) > 1e-6 + end +end + +@testset "∇∇ on skeleton integrals: Hessian only in the Jacobian, multi-field, update, matrix-free rebind" begin model = cartesian_model(2, (3, 3); simplex=true) Ω = Triangulation(model); dΩ = Measure(Ω, 4) Λ = SkeletonTriangulation(model); dΛ = Measure(Λ, 4) - V = TestFESpace(model, ReferenceFE(lagrangian, Float64, 2), dirichlet_tags="boundary") - U = TrialFESpace(V, x -> x[1]^2) - res(u, v) = ∫(∇∇(v) ⊙ ∇∇(u) + v * (u * u))dΩ + ∫(mean(∇∇(u)) ⊙ mean(∇∇(v)) + jump(u) * jump(v))dΛ - jac(u, du, v) = ∫(∇∇(v) ⊙ ∇∇(du) + 2.0 * (v * (u * du)))dΩ + ∫(mean(∇∇(du)) ⊙ mean(∇∇(v)) + jump(du) * jump(v))dΛ - op = @test_logs (:warn, r"second derivatives \(∇∇\) on skeleton integrals are not supported yet") match_mode = :any KernelOperator(res, jac, U, V) - @test op.fallback.rtrians == [Λ] && op.fallback.jtrians == [Λ] - @test length(op.rgroups) == 1 && op.rgroups[1].trian === Ω - c = check(res, jac, U, V, randx(U); op) + n = get_normal_vector(Λ) + Vu = TestFESpace(model, ReferenceFE(lagrangian, VectorValue{2,Float64}, 2), dirichlet_tags="boundary") + Vφ = TestFESpace(model, ReferenceFE(lagrangian, Float64, 2), dirichlet_tags="boundary") + V = MultiFieldFESpace([Vu, Vφ]) + U = MultiFieldFESpace([TrialFESpace(Vu, x -> VectorValue(0.01x[2], 0.0)), TrialFESpace(Vφ, x -> x[1])]) + # coupling on the skeleton: mean(Δφ) with jump(∇v), jump(∇∇u) with mean(φ) + mk(γ) = (((u, φ), (v, vφ)) -> ∫(∇(v)' ⊙ (P_svk ∘ (Fk ∘ (∇(u)'))) + ∇(vφ) ⋅ ∇(φ))dΩ + + ∫(γ * (jump(∇∇(u)) ⊙ jump(∇∇(v))) * (1.0 + mean(φ) * mean(φ)) + mean(Δ(φ)) * (jump(∇(v)) ⊙ jump(∇(u))) + + jump(∇(vφ) ⋅ n) * jump(∇(φ) ⋅ n) + mean(Δ(vφ)) * mean(Δ(φ)))dΛ, + ((u, φ), (du, dφ), (v, vφ)) -> ∫(∇(v)' ⊙ (dP_svk ∘ (Fk ∘ (∇(u)'), ∇(du)')) + ∇(vφ) ⋅ ∇(dφ))dΩ + + ∫(γ * (jump(∇∇(du)) ⊙ jump(∇∇(v))) * (1.0 + mean(φ) * mean(φ)) + + γ * (jump(∇∇(u)) ⊙ jump(∇∇(v))) * (2.0 * mean(φ) * mean(dφ)) + + mean(Δ(dφ)) * (jump(∇(v)) ⊙ jump(∇(u))) + mean(Δ(φ)) * (jump(∇(v)) ⊙ jump(∇(du))) + + jump(∇(vφ) ⋅ n) * jump(∇(dφ) ⋅ n) + mean(Δ(vφ)) * mean(Δ(dφ)))dΛ) + r1, j1 = mk(0.3); r2, j2 = mk(1.7) + for linear_maps in (true, false) + c = compare((; U, V, res=r1, jac=j1), r1, j1; linear_maps) + @test c.res < HTOL && c.jac < HTOL && c.mf < 1e-11 + end + # Hessian only in the Jacobian (a residual without skeleton Hessians) + r0((u, φ), (v, vφ)) = ∫(∇(v)' ⊙ (P_svk ∘ (Fk ∘ (∇(u)'))) + ∇(vφ) ⋅ ∇(φ))dΩ + ∫(jump(∇(vφ) ⋅ n) * jump(∇(φ) ⋅ n))dΛ + j0((u, φ), (du, dφ), (v, vφ)) = ∫(∇(v)' ⊙ (dP_svk ∘ (Fk ∘ (∇(u)'), ∇(du)')) + ∇(vφ) ⋅ ∇(dφ))dΩ + + ∫(jump(∇(vφ) ⋅ n) * jump(∇(dφ) ⋅ n) + 0.5 * (jump(∇∇(dφ)) ⊙ jump(∇∇(vφ))))dΛ + c0 = compare((; U, V, res=r0, jac=j0), r0, j0) + @test c0.res < HTOL && c0.jac < HTOL && c0.mf < 1e-11 + @test all(g -> !(g isa HK.SkeletonGroup) || HK.hess_use(g) isa Val{false}, c0.op.rgroups) + @test any(g -> g isa HK.SkeletonGroup && HK.hess_use(g) isa Val{true}, c0.op.jgroups) + # update with another γ: same compiled kernels and tables; matrix-free rebound + op1 = KernelOperator(r1, j1, U, V) + op2 = update(op1, r2, j2) + for (g1, g2) in zip(vcat(op1.rgroups, op1.jgroups), vcat(op2.rgroups, op2.jgroups)) + @test typeof(g2) === typeof(g1) && typeof(g2.pf.core) === typeof(g1.pf.core) && g2.tables === g1.tables + end + x = randx(U); dx = randn(num_free_dofs(U)) + c = check(r2, j2, U, V, x; op=op2) @test c.res < HTOL && c.jac < HTOL + mf = MatrixFreeJacobian(op1) + jacobian!(mf, op2, x) + @test reldiff(mf * dx, Algebra.jacobian(op2, x) * dx) < 1e-11 + @test sk_noalloc(op2, x) end From 3ccfecd10e4e9449a2e3048b2c6fd4a7a46c4b4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jesus=20Mart=C3=ADnez=20Frutos?= Date: Tue, 10 Nov 2026 18:03:57 +0100 Subject: [PATCH 2/3] KernelFunctional: skeleton integrals evaluated by the facet kernels MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Energies with jump penalties (½∫c⟦∇u⟧:⟦∇u⟧, ½∫c⟦∇∇u⟧:⟦∇∇u⟧ on a skeleton, the Armijo energy of the third-medium regularisation) were summed by the Gridap fallback: compiling the skeleton failed because the facet colouring needs the test rows of an operator. compile_groups(...; facet_colors = false) builds the skeleton data with one colour (a reduction has no races); functional_cells for a SkeletonGroup gathers both cells and evaluates the integrand at the facet points with the values, gradients and Hessians of both sides (:sum with the facet measure, :min/:max over the points). The functional keeps its skeleton data for update; aux_data of a skeleton class serves functionals too. Tests (test_functionals.jl): the energy of the regularisation on the skeleton of a subdomain plus a mean term, against Gridap on triangles, Cartesian and non-affine quadrilaterals; update reuses groups, tables and skeleton data; threads; gradient equals the operator's residual; :min/:max over the facet points; no allocations in the facet loop. --- CHANGELOG.md | 3 ++ src/assembly.jl | 4 +- src/functionals.jl | 52 +++++++++++++++++++++----- src/skeleton.jl | 13 +++++-- test/solvers/test_functionals.jl | 63 ++++++++++++++++++++++++++++++++ 5 files changed, 120 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c6e913a..44e85e0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,9 @@ Hessian-jump regularisations ⟦∇∇u⟧:⟦∇∇v⟧), with the Hessian push-forward of each side (its Jt⁻¹ and the curvature of its cell map: exact on non-affine quadrilaterals and hexahedra); residual, assembled and matrix-free Jacobian. They no longer send the skeleton to the Gridap fallback. +- `KernelFunctional` evaluates skeleton integrals with the facet kernels (energies of jump penalties such + as `½ ∫ c ⟦∇u⟧:⟦∇u⟧ dΛ` and `½ ∫ c ⟦∇∇u⟧:⟦∇∇u⟧ dΛ`, used by the Armijo line search; `:min`/`:max` over the + facet points). Before, those terms were summed by the Gridap fallback. ## 0.2.2 diff --git a/src/assembly.jl b/src/assembly.jl index 512fa8c..dc08021 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -288,7 +288,7 @@ quadrature or operation) are added to `fallback` with the reason and are assembl function compile_groups(forms::Vector{Form}, U, locs, D::Int, colors, previous=Any[]; linear_maps::Bool=true, fallback=Tuple{Any,String}[], skel::Vector{SkeletonData}=SkeletonData[], previous_skel::Vector{SkeletonData}=SkeletonData[], - fields::Vector{FieldLayout}=FieldLayout[], nrows::Int=0, ndl::Int=0) + fields::Vector{FieldLayout}=FieldLayout[], nrows::Int=0, ndl::Int=0, facet_colors::Bool=true) groups = Any[] Us = field_spaces(U) strian = get_triangulation(Us[1]) @@ -302,7 +302,7 @@ function compile_groups(forms::Vector{Form}, U, locs, D::Int, colors, previous=A k = findfirst(s -> s.trian === form.trian, skel) if k === nothing kp = findfirst(s -> s.trian === form.trian, previous_skel) - sd = kp === nothing ? skeleton_data(form.trian, strian, fields, nrows) : + sd = kp === nothing ? skeleton_data(form.trian, strian, fields, nrows; colored=facet_colors) : SkeletonData(form.trian, previous_skel[kp].cells, previous_skel[kp].colors, previous_skel[kp].nzpos) push!(skel, sd) else diff --git a/src/functionals.jl b/src/functionals.jl index 3f002ca..187d606 100644 --- a/src/functionals.jl +++ b/src/functionals.jl @@ -3,7 +3,9 @@ # a line search, `min_q α(∇u, ∇du)`). The functional is captured like a residual without test # function: its terms are compiled into a pointwise function of the field values and the auxiliary # data (FE functions in the form read their nodal values at every evaluation, so an increment -# stored in an FE function can change between calls without recompiling). +# stored in an FE function can change between calls without recompiling). Skeleton integrals +# (`jump`, `mean`, `plus`/`minus`, with gradients and Hessians of both sides) are evaluated per facet +# with the facet kernels of the operators. """ KernelFunctional(f, U; reduce = :sum, nthreads = Threads.nthreads()) @@ -18,6 +20,9 @@ quantity of the state and an increment). The functional is called on the free va - `reduce = :min` / `:max`: the extremum of the integrand over the quadrature points of the measure, without weights (the injectivity step `min_q α(∇u_q, ∇du_q)` of `InjectivityLS`). +Cell, boundary and skeleton integrals (e.g. the energy `½ ∫ c ⟦∇u⟧:⟦∇u⟧ dΛ` of a jump penalty) are +evaluated by the kernels. + Contributions the kernels cannot read are summed by Gridap (`:sum` only; `:min`/`:max` need every contribution in the kernels). [`update`](@ref) gives a functional with another `f` on the same space reusing the tables. @@ -32,6 +37,7 @@ struct KernelFunctional{F,US} nthreads::Int cache_aux::Bool fallback::Vector{Any} # triangulations whose contributions Gridap sums (`:sum`) + skel::Vector{SkeletonData} # skeletons of the form (cells per facet, one colour), reused by `update` signed::Bool sig::Any # structural signature of the capture (`signature`; cheap `update`) end @@ -57,15 +63,16 @@ function KernelFunctional(f, U::FESpace; reduce::Symbol=:sum, nthreads::Int=Thre end forms, fb = capture(f, probes(U, :unknown)) sig = signature(forms, fb) - colors = [Int32.(1:ncells)] # a reduction has no races: one colour - groups = compile_groups(forms, U, Tuple(fl.loc for fl in fields), D, colors; linear_maps=false, fallback=fb, fields, ndl=loc) + colors = [Int32.(1:ncells)] # a reduction has no races: one colour (cells and facets) + skel = SkeletonData[] + groups = compile_groups(forms, U, Tuple(fl.loc for fl in fields), D, colors; linear_maps=false, fallback=fb, skel, fields, ndl=loc, + facet_colors=false) _check_functional(groups, fb, reduce) - KernelFunctional(f, U, fields, loc, groups, reduce, nthreads, false, Any[t for (t, _) in fb], any(has_signs, fields), sig) + KernelFunctional(f, U, fields, loc, groups, reduce, nthreads, false, Any[t for (t, _) in fb], skel, any(has_signs, fields), sig) end function _check_functional(groups, fb, reduce) for g in groups - g isa Group || error("HyperFEMKernels: KernelFunctional does not support skeleton integrals yet") nslots(g) == (0, 0) || error("HyperFEMKernels: a KernelFunctional has no test or trial functions (got a form with them)") end if !isempty(fb) @@ -93,12 +100,15 @@ function update(fun::KernelFunctional, f) if sig == fun.sig groups = rebind_groups(fun.groups, forms) check_rebind() && verify_rebind(groups, compile_groups(forms, U, locs, D, colors, fun.groups; linear_maps=false, fallback=copy(fb), - fields=fun.fields, ndl=fun.ndl), "functional") - return KernelFunctional(f, U, fun.fields, fun.ndl, groups, fun.reduce, fun.nthreads, fun.cache_aux, fun.fallback, fun.signed, sig) + skel=SkeletonData[], previous_skel=fun.skel, fields=fun.fields, ndl=fun.ndl, + facet_colors=false), "functional") + return KernelFunctional(f, U, fun.fields, fun.ndl, groups, fun.reduce, fun.nthreads, fun.cache_aux, fun.fallback, fun.skel, fun.signed, sig) end - groups = compile_groups(forms, U, locs, D, colors, fun.groups; linear_maps=false, fallback=fb, fields=fun.fields, ndl=fun.ndl) + skel = SkeletonData[] + groups = compile_groups(forms, U, locs, D, colors, fun.groups; linear_maps=false, fallback=fb, skel, previous_skel=fun.skel, + fields=fun.fields, ndl=fun.ndl, facet_colors=false) _check_functional(groups, fb, fun.reduce) - KernelFunctional(f, U, fun.fields, fun.ndl, groups, fun.reduce, fun.nthreads, fun.cache_aux, Any[t for (t, _) in fb], fun.signed, sig) + KernelFunctional(f, U, fun.fields, fun.ndl, groups, fun.reduce, fun.nthreads, fun.cache_aux, Any[t for (t, _) in fb], skel, fun.signed, sig) end _init(::Val{:sum}) = 0.0 @@ -167,3 +177,27 @@ function functional_cells(R::Val, fun::KernelFunctional, g::Group{NSV,NSU,PF,TT, end acc end + +"Facets `lfacets` (local indices of the class) of a skeleton class: both cells gathered, the integrand at the facet points." +function functional_cells(R::Val, fun::KernelFunctional, g::SkeletonGroup{NSV,NSU,PF,TT,AUX,SV,SU,D}, x, dvals, auxd, + lfacets) where {NSV,NSU,PF,TT,AUX,SV,SU,D} + ndl = fun.ndl + xl = zeros(2ndl) + tables, locs, pf = g.tables, g.locs, g.pf + skel, facets = g.skel, g.facets + hu = hess_use(g) + acc = _init(R) + for lf in lfacets + f = facets[lf] + gather_cell!(xl, fun.fields, x, dvals, skel.cells[1, f]) + gather_cell!(xl, fun.fields, x, dvals, skel.cells[2, f], ndl) + for q in 1:g.geo.np + dS, iJs, Φs = side_geometry(g, q, lf, hu) + FV = field_values_sides(tables, locs, g.use, xl, q, iJs, Φs) + Aq = map(a -> @inbounds(a[q, lf]), auxd) + v = pf.core(pf.ops, pf.coefs, FV, Aq, EMPTY_SLOTS, EMPTY_SLOTS)::Float64 + acc = _combine(R, acc, R isa Val{:sum} ? v * dS : v) + end + end + acc +end diff --git a/src/skeleton.jl b/src/skeleton.jl index 6bb0cfc..5402ed6 100644 --- a/src/skeleton.jl +++ b/src/skeleton.jl @@ -63,10 +63,15 @@ function color_facets(fields::Vector{FieldLayout}, cells::Matrix{Int32}, nrows:: groups end -"Data of a skeleton for the FE spaces on `space_trian`: cells per facet and facet colouring (positions come with the pattern)." -function skeleton_data(trian::Triangulation, space_trian::Triangulation, fields::Vector{FieldLayout}, nrows::Int) +""" +Data of a skeleton for the FE spaces on `space_trian`: cells per facet and facet colouring (positions come +with the pattern). `colored = false` (a reduction, `KernelFunctional`: no scatter, no races) puts every +facet in one colour. +""" +function skeleton_data(trian::Triangulation, space_trian::Triangulation, fields::Vector{FieldLayout}, nrows::Int; colored::Bool=true) cells = skeleton_cells(trian, space_trian) - SkeletonData(trian, cells, color_facets(fields, cells, nrows), Matrix{Int32}(undef, 0, 0)) + colors = colored ? color_facets(fields, cells, nrows) : [Int32.(1:size(cells, 2))] + SkeletonData(trian, cells, colors, Matrix{Int32}(undef, 0, 0)) end """ @@ -222,7 +227,7 @@ function _materialize_facets!(M::Matrix, a, facets, range, np::Int) M end -function aux_data(op::KernelOperator, g::SkeletonGroup{NSV,NSU,PF,TT,AUX}) where {NSV,NSU,PF,TT,AUX} +function aux_data(op, g::SkeletonGroup{NSV,NSU,PF,TT,AUX}) where {NSV,NSU,PF,TT,AUX} # op: KernelOperator or KernelFunctional if g.auxcache[] === nothing g.auxcache[] = materialize_aux(g, op.nthreads) elseif !op.cache_aux diff --git a/test/solvers/test_functionals.jl b/test/solvers/test_functionals.jl index 049fe9e..9ee098b 100644 --- a/test/solvers/test_functionals.jl +++ b/test/solvers/test_functionals.jl @@ -86,6 +86,69 @@ end @test fc(xc) ≈ lazy_extremum(g ∘ (∇(uh[1])', ∇(duc[1])') + uh[2] * uh[2], dΩc, min) rtol = 1e-13 end +# Skeleton integrals: the energy of the third-medium regularisation (jump penalties of ∇u and ∇∇u on the skeleton of a +# subdomain, facet size as data) and a mean term on the whole skeleton; non-affine quadrilaterals against hessian_nonaffine. +function skeleton_energy_problem(; simplex=true, mapped=false) + model = cartesian_model(2, (4, 3); simplex, mapped) + Ω = Triangulation(model); dΩ = Measure(Ω, 4) + cc = [sum(c) / length(c) for c in get_cell_coordinates(model)] + Ωc = Triangulation(model, [x[1] > 0.4 for x in cc]); dΩc = Measure(Ωc, 4) + Λ = SkeletonTriangulation(Ωc); dΛ = Measure(Λ, 4) + Λa = SkeletonTriangulation(model); dΛa = Measure(Λa, 4) + hF = CellField(get_array(∫(1)dΛ), Λ) + V = TestFESpace(model, ReferenceFE(lagrangian, VectorValue{2,Float64}, 2), dirichlet_tags="boundary") + U = TrialFESpace(V, x -> VectorValue(0.01x[2], -0.02x[1] * x[2])) + mk(H) = β -> u -> ∫(Ψ_svk ∘ (Fk ∘ (∇(u)')))dΩ + (0.5β) * ∫(H(u) ⊙ H(u))dΩc + + ∫((0.5β / hF) * (jump(∇(u)) ⊙ jump(∇(u))) + (β * hF) * (jump(H(u)) ⊙ jump(H(u))))dΛ + ∫(0.1 * (mean(u) ⋅ mean(u)))dΛa + # the residual of the energy (gradient), for the finite-difference check + res(β) = (u, v) -> ∫(∇(v)' ⊙ (P_svk ∘ (Fk ∘ (∇(u)'))))dΩ + β * ∫(∇∇(u) ⊙ ∇∇(v))dΩc + + ∫((β / hF) * (jump(∇(u)) ⊙ jump(∇(v))) + (2β * hF) * (jump(∇∇(u)) ⊙ jump(∇∇(v))))dΛ + ∫(0.2 * (mean(u) ⋅ mean(v)))dΛa + jac(β) = (u, du, v) -> ∫(∇(v)' ⊙ (dP_svk ∘ (Fk ∘ (∇(u)'), ∇(du)')))dΩ + β * ∫(∇∇(du) ⊙ ∇∇(v))dΩc + + ∫((β / hF) * (jump(∇(du)) ⊙ jump(∇(v))) + (2β * hF) * (jump(∇∇(du)) ⊙ jump(∇∇(v))))dΛ + ∫(0.2 * (mean(du) ⋅ mean(v)))dΛa + (; model, Λ, dΛ, U, V, en=mk(∇∇), enref=mk(mapped ? hessian_nonaffine : ∇∇), res, jac) +end + +@testset "skeleton integrals (jumps of ∇u and ∇∇u, mean): energy, update, threads, gradient, extrema" begin + for (simplex, mapped) in ((true, false), (false, false), (false, true)) + s = skeleton_energy_problem(; simplex, mapped) + fun = KernelFunctional(s.en(0.05), s.U) + @test isempty(fun.fallback) && count(g -> g isa HK.SkeletonGroup, fun.groups) >= 2 + @test length(fun.skel) == 2 && all(sd -> length(sd.colors) == 1, fun.skel) # a reduction: one colour per skeleton + x = randx(s.U; scale=0.1) + @test fun(x) ≈ lazy_sum(s.enref(0.05), s.U, x) rtol = 1e-12 + # update with another β: groups rebound, tables and skeleton data reused + fun2 = update(fun, s.en(0.3)) + @test fun2(x) ≈ lazy_sum(s.enref(0.3), s.U, x) rtol = 1e-12 + @test fun2.skel === fun.skel && all(((g1, g2),) -> g1.tables === g2.tables, zip(fun.groups, fun2.groups)) + @test KernelFunctional(s.en(0.05), s.U; nthreads=1)(x) ≈ fun(x) rtol = 1e-13 + if !mapped + # the gradient of the energy is the residual (finite differences along a random direction) + op = KernelOperator(s.res(0.05), s.jac(0.05), s.U, s.V) + dx = randn(num_free_dofs(s.U)); h = 1e-6 + @test (fun(x + h * dx) - fun(x - h * dx)) / 2h ≈ dot(residual(op, x), dx) rtol = 1e-6 + end + end + # pointwise extrema over the facet points + s = skeleton_energy_problem() + g(a, b) = a / (1.0 + b * b) + f(u) = ∫(g ∘ (jump(∇(u)) ⊙ jump(∇(u)), mean(tr ∘ ∇(u))))s.dΛ + x = randx(s.U; scale=0.1); uh = FEFunction(s.U, x) + cf = g ∘ (jump(∇(uh)) ⊙ jump(∇(uh)), mean(tr ∘ ∇(uh))) + @test KernelFunctional(f, s.U; reduce=:min)(x) ≈ lazy_extremum(cf, s.dΛ, min) rtol = 1e-13 + @test KernelFunctional(f, s.U; reduce=:max)(x) ≈ lazy_extremum(cf, s.dΛ, max) rtol = 1e-13 + # the facet loop does not allocate (same bytes for one facet and for all of them) + fun = KernelFunctional(s.en(0.05), s.U) + fun(x) + gs = first(gg for gg in fun.groups if gg isa HK.SkeletonGroup) + dvals = [get_dirichlet_dof_values(sp) for sp in HK.field_spaces(s.U)] + auxd = HK.aux_data(fun, gs) + one_facet, all_facets = Int32[1], gs.colors[1] + HK.functional_cells(Val(:sum), fun, gs, x, dvals, auxd, one_facet) + a1 = @allocated HK.functional_cells(Val(:sum), fun, gs, x, dvals, auxd, one_facet) + a2 = @allocated HK.functional_cells(Val(:sum), fun, gs, x, dvals, auxd, all_facets) + @test a1 == a2 +end + @testset "errors and allocations" begin p = poisson_2d(; order=2); dΩp = p.dΩ @test_throws ArgumentError KernelFunctional(u -> ∫(u * u)dΩp, p.U; reduce=:mean) From a914fae3d3d978d8403ef2477776ebb9b033d4d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jesus=20Mart=C3=ADnez=20Frutos?= Date: Wed, 11 Nov 2026 08:09:30 +0100 Subject: [PATCH 3/3] CI: Hessian tests in their own job (HFK_TEST_GROUP=hessian) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The skeleton ∇∇ tests pushed the core job past the 90-minute limit on ubuntu and on macOS with Julia 1.12 (29 → 82 min before the solver tests on ubuntu). test_hessian.jl now runs as a third group, as test_hdiv.jl does; the CI check requires the new job. --- .github/workflows/CI.yml | 36 +++++++++++++++++++++++++++++------- test/README.md | 2 +- test/runtests.jl | 11 +++++++---- 3 files changed, 37 insertions(+), 12 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 925ebf3..f9b76e0 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -30,9 +30,9 @@ jobs: - '!**/*.md' - '!docs/**' - '!LICENSE' - # The engine tests run in two groups (HFK_TEST_GROUP, see test/runtests.jl) so that each job stays - # well below the time limit: `test` runs everything but the H(div) tests (the job names are the - # required checks of main), `test-hdiv` runs the H(div) tests. + # The engine tests run in three groups (HFK_TEST_GROUP, see test/runtests.jl) so that each job stays + # well below the time limit: `test` runs everything but the H(div) and Hessian tests, `test-hdiv` the + # H(div) tests, `test-hessian` the second-derivative (∇∇) tests. test: name: Julia ${{ matrix.version }} - ${{ matrix.os }} needs: changes @@ -77,6 +77,28 @@ jobs: env: JULIA_NUM_THREADS: 2 HFK_TEST_GROUP: hdiv + test-hessian: + name: Hessian tests - Julia ${{ matrix.version }} - ${{ matrix.os }} + needs: changes + if: needs.changes.outputs.code == 'true' + runs-on: ${{ matrix.os }} + timeout-minutes: 90 + strategy: + fail-fast: false + matrix: + version: ['1.12', '1'] + os: [ubuntu-latest, macos-latest] + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + env: + JULIA_NUM_THREADS: 2 + HFK_TEST_GROUP: hessian hyperfem: name: HyperFEM weak forms - Julia ${{ matrix.version }} needs: changes @@ -109,16 +131,16 @@ jobs: ci: name: CI if: always() - needs: [changes, test, test-hdiv, hyperfem] + needs: [changes, test, test-hdiv, test-hessian, hyperfem] runs-on: ubuntu-latest steps: - name: Check results run: | echo "code changed: ${{ needs.changes.outputs.code }}" - echo "test: ${{ needs.test.result }}, test-hdiv: ${{ needs.test-hdiv.result }}, hyperfem: ${{ needs.hyperfem.result }}" + echo "test: ${{ needs.test.result }}, test-hdiv: ${{ needs.test-hdiv.result }}, test-hessian: ${{ needs.test-hessian.result }}, hyperfem: ${{ needs.hyperfem.result }}" [ "${{ needs.changes.result }}" = success ] || exit 1 if [ "${{ needs.changes.outputs.code }}" = true ]; then - [ "${{ needs.test.result }}" = success ] && [ "${{ needs.test-hdiv.result }}" = success ] && [ "${{ needs.hyperfem.result }}" = success ] || exit 1 + [ "${{ needs.test.result }}" = success ] && [ "${{ needs.test-hdiv.result }}" = success ] && [ "${{ needs.test-hessian.result }}" = success ] && [ "${{ needs.hyperfem.result }}" = success ] || exit 1 else - [ "${{ needs.test.result }}" = skipped ] && [ "${{ needs.test-hdiv.result }}" = skipped ] && [ "${{ needs.hyperfem.result }}" = skipped ] || exit 1 + [ "${{ needs.test.result }}" = skipped ] && [ "${{ needs.test-hdiv.result }}" = skipped ] && [ "${{ needs.test-hessian.result }}" = skipped ] && [ "${{ needs.hyperfem.result }}" = skipped ] || exit 1 fi diff --git a/test/README.md b/test/README.md index d423880..e0bd8b2 100644 --- a/test/README.md +++ b/test/README.md @@ -1,6 +1,6 @@ # Tests -`julia --project -e 'using Pkg; Pkg.test()'`. `HFK_TEST_GROUP=hdiv` runs only `test_hdiv.jl`, `HFK_TEST_GROUP=core` everything else (default `all`); CI runs the two groups as separate jobs. Every test compares against Gridap or against explicit sums; the +`julia --project -e 'using Pkg; Pkg.test()'`. `HFK_TEST_GROUP=hdiv` runs only `test_hdiv.jl`, `HFK_TEST_GROUP=hessian` only `test_hessian.jl`, `HFK_TEST_GROUP=core` everything else (default `all`); CI runs the three groups as separate jobs. Every test compares against Gridap or against explicit sums; the tolerance is 1e-12 relative. | File | What it checks | diff --git a/test/runtests.jl b/test/runtests.jl index 2d1e885..f5364c9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -10,10 +10,11 @@ const TOL = 1e-12 include("common.jl") -# Test groups (`HFK_TEST_GROUP`): "all" (default), "hdiv" (H(div) fields only) or "core" (everything -# else). CI runs "core" and "hdiv" as separate jobs so that each stays well below the time limit. +# Test groups (`HFK_TEST_GROUP`): "all" (default), "hdiv" (H(div) fields only), "hessian" (second +# derivatives only) or "core" (everything else). CI runs "core", "hdiv" and "hessian" as separate jobs so +# that each stays well below the time limit. const TEST_GROUP = get(ENV, "HFK_TEST_GROUP", "all") -TEST_GROUP in ("all", "core", "hdiv") || error("HFK_TEST_GROUP must be all, core or hdiv (got $TEST_GROUP)") +TEST_GROUP in ("all", "core", "hdiv", "hessian") || error("HFK_TEST_GROUP must be all, core, hdiv or hessian (got $TEST_GROUP)") const TEST_FILES = [ ("probes", "test_probes.jl"), ("tree", "test_tree.jl"), ("codegen", "test_codegen.jl"), ("tables", "test_tables.jl"), ("kernels", "test_kernels.jl"), ("linear maps", "test_linearmaps.jl"), @@ -22,7 +23,9 @@ const TEST_FILES = [ ("CFD (Oseen step)", "test_cfd.jl"), ("subdomains, interface", "test_subdomains.jl"), ("skeleton", "test_skeleton.jl"), ("Hessians (∇∇)", "test_hessian.jl"), ("cheap update", "test_update.jl"), ("solvers: operator mode", "solvers/test_operator.jl"), ("solvers: matrix-free solver", "solvers/test_matrixfree_solver.jl"), ("solvers: functionals", "solvers/test_functionals.jl"), ("solvers: p-multigrid transfers", "solvers/test_transfers.jl"), ("solvers: p-multigrid levels", "solvers/test_pmultigrid_levels.jl"), ("solvers: Chebyshev smoother", "solvers/test_chebyshev.jl"), ("solvers: p-multigrid", "solvers/test_pmultigrid.jl"), ("manifold spaces", "test_manifold.jl"), ("FE functions in the form", "test_fefunctions.jl"), ("states", "test_states.jl"), ("Gridap equality", "test_gridap_equality.jl"), ] -in_group(file) = TEST_GROUP == "all" || (TEST_GROUP == "hdiv") == (file == "test_hdiv.jl") +const GROUP_FILES = Dict("hdiv" => "test_hdiv.jl", "hessian" => "test_hessian.jl") +file_group(file) = something(findfirst(==(file), GROUP_FILES), "core") +in_group(file) = TEST_GROUP == "all" || TEST_GROUP == file_group(file) @testset "HyperFEMKernels ($TEST_GROUP)" begin for (name, file) in TEST_FILES