Second derivatives (∇∇, Δ) on skeleton integrals; skeleton integrals in KernelFunctional

.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

CHANGELOG.md

@@ -1,5 +1,16 @@
 # 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.
+- `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
 
 - p-multigrid preconditioner `PMultigrid` for `MatrixFreeSolver`: nested Lagrangian levels on the same mesh

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

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

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

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
 

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))