Merge pull request #799 from kishore-nori/AD-for-SkeletonIntegration

Extended `gradient` for `MultiField` functionals involving Skeleton integration

Author: amartinhuertas

NEWS.md

@@ -7,7 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## Unreleased
 
 ### Added
-- Added functionality to take gradient of functional involving integration (`DomainContribution`) over Skeleton faces, with respect to the degrees of freedom `FEFunction`.  The interface remains the same - `gradient(f,uh)`. Since PR [#797](https://github.com/gridap/Gridap.jl/pull/797)
+- Added functionality to take gradient of functional involving integration (`DomainContribution`) over Skeleton faces, with respect to the degrees-of-freedom `FEFunction`.  The interface remains the same - `gradient(f,uh)`. Since PR [#797](https://github.com/gridap/Gridap.jl/pull/797)
+- Extended the `MultiField` functional gradient (with respect to degrees-of-freedom of `MultiFieldFEFunction`) to functionals involving Skeleton integration. The interface remains the same `gradient(f,xh)`. Since PR [#799](https://github.com/gridap/Gridap.jl/pull/799)
 
 ## [0.17.13] - 2022-05-31
 

src/MultiField/MultiField.jl

@@ -18,6 +18,8 @@ using Gridap.Fields
 using Gridap.FESpaces: FEBasis, TestBasis, TrialBasis, get_cell_dof_values
 using Gridap.FESpaces: SingleFieldFEBasis, TestBasis, TrialBasis
 using Gridap.CellData: CellFieldAt
+using Gridap.CellData: SkeletonCellFieldPair
+
 import Gridap.Fields: gradient, DIV, ∇∇
 
 using ForwardDiff

src/MultiField/MultiFieldFEAutodiff.jl

@@ -15,6 +15,22 @@ function _get_cell_dofs_field_offsets(uh::MultiFieldFEFunction)
   dofs_field_offsets
 end
 
+function _restructure_cell_grad!(
+  cell_grad, uh::MultiFieldFEFunction, trian)
+
+  monolithic_result=cell_grad
+  blocks = [] # TO-DO type unstable. How can I infer the type of its entries?
+  nfields = length(uh.fe_space.spaces)
+  cell_dofs_field_offsets=_get_cell_dofs_field_offsets(uh)
+  for i in 1:nfields
+    view_range=cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
+    block=lazy_map(x->view(x,view_range),monolithic_result)
+    append!(blocks,[block])
+  end
+  cell_grad=lazy_map(BlockMap(nfields,collect(1:nfields)),blocks...)
+  cell_grad
+end
+
 function FESpaces._gradient(f,uh::MultiFieldFEFunction,fuh::DomainContribution)
   terms = DomainContribution()
   U = get_fe_space(uh)
@@ -23,16 +39,7 @@ function FESpaces._gradient(f,uh::MultiFieldFEFunction,fuh::DomainContribution)
     cell_u = lazy_map(DensifyInnerMostBlockLevelMap(),get_cell_dof_values(uh))
     cell_id = FESpaces._compute_cell_ids(uh,trian)
     cell_grad = autodiff_array_gradient(g,cell_u,cell_id)
-    monolithic_result=cell_grad
-    blocks = [] # TO-DO type unstable. How can I infer the type of its entries?
-    nfields = length(U.spaces)
-    cell_dofs_field_offsets=_get_cell_dofs_field_offsets(uh)
-    for i in 1:nfields
-      view_range=cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
-      block=lazy_map(x->view(x,view_range),monolithic_result)
-      append!(blocks,[block])
-    end
-    cell_grad=lazy_map(BlockMap(nfields,collect(1:nfields)),blocks...)
+    cell_grad = _restructure_cell_grad!(cell_grad, uh, trian)
     add_contribution!(terms,trian,cell_grad)
   end
   terms
@@ -147,3 +154,62 @@ function FESpaces._hessian(f,uh::MultiFieldFEFunction,fuh::DomainContribution)
   end
   terms
 end
+
+# overloads for AD of SkeletonTriangulation DomainContribution with MultiFields
+
+function FESpaces._change_argument(
+  op::typeof(gradient),f,trian::SkeletonTriangulation,uh::MultiFieldFEFunction)
+
+  U = get_fe_space(uh)
+  function g(cell_u)
+    single_fields_plus = SkeletonCellFieldPair[]
+    single_fields_minus = SkeletonCellFieldPair[]
+    nfields = length(U.spaces)
+    cell_dofs_field_offsets=_get_cell_dofs_field_offsets(uh)
+    for i in 1:nfields
+      view_range=cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
+      cell_values_field = lazy_map(a->view(a,view_range),cell_u)
+      cf_dual = CellField(U.spaces[i],cell_values_field)
+      scfp_plus = SkeletonCellFieldPair(cf_dual, uh[i])
+      scfp_minus = SkeletonCellFieldPair(uh[i], cf_dual)
+      push!(single_fields_plus,scfp_plus)
+      push!(single_fields_minus,scfp_minus)
+    end
+    xh_plus = MultiFieldCellField(single_fields_plus)
+    xh_minus = MultiFieldCellField(single_fields_minus)
+    cell_grad_plus = f(xh_plus)
+    cell_grad_minus = f(xh_minus)
+    get_contribution(cell_grad_plus,trian), get_contribution(cell_grad_minus,trian)
+  end
+  g
+end
+
+function FESpaces._change_argument(
+  op::typeof(jacobian),f,trian::SkeletonTriangulation,uh::MultiFieldFEFunction)
+
+  @notimplemented
+end
+
+function FESpaces._change_argument(
+  op::typeof(hessian),f,trian::SkeletonTriangulation,uh::MultiFieldFEFunction)
+
+  @notimplemented
+end
+
+function _restructure_cell_grad!(
+  cell_grad, uh::MultiFieldFEFunction, trian::SkeletonTriangulation)
+
+  monolithic_result=cell_grad
+  blocks = [] # TO-DO type unstable. How can I infer the type of its entries?
+  nfields = length(uh.fe_space.spaces)
+  cell_dofs_field_offsets=_get_cell_dofs_field_offsets(uh)
+  for i in 1:nfields
+    view_range=cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
+    block_plus = lazy_map(x->view(x[1],view_range),monolithic_result)
+    block_minus = lazy_map(x->view(x[2],view_range),monolithic_result)
+    block = lazy_map(BlockMap(2,[1,2]),block_plus,block_minus)
+    push!(blocks,block)
+  end
+  cell_grad=lazy_map(BlockMap(nfields,collect(1:nfields)),blocks...)
+  cell_grad
+end

test/CellDataTests/SkeletonCellFieldPairTests.jl

@@ -9,6 +9,7 @@ using Gridap.Fields
 using Gridap.ReferenceFEs
 using Gridap.Geometry
 using Gridap.CellData
+using Gridap.MultiField
 using ForwardDiff
 
 model = CartesianDiscreteModel((0.,1.,0.,1.),(3,3))
@@ -88,4 +89,88 @@ scfp_dualplus = SkeletonCellFieldPair(uh_dual,uh)
 
 # AD of Integrals over Skeleton faces tested in FESpaces -> FEAutodiffTests.jl
 
+# testing the MultiField SkeletonCellFieldPair (SCFP)
+# this uses some low-level API to construct MultiFieldCellFields with SCFP
+
+Q = TestFESpace(model,reffe,conformity=:L2)
+P = TrialFESpace(Q)
+
+Y = MultiFieldFESpace([V, Q])
+X = MultiFieldFESpace([U, P])
+
+xh = FEFunction(X,rand(num_free_dofs(X)))
+uh,ph = xh
+
+cellu = lazy_map(DensifyInnerMostBlockLevelMap(),get_cell_dof_values(xh))
+cellu_dual = lazy_map(DualizeMap(ForwardDiff.gradient),cellu)
+single_fields_plus = SkeletonCellFieldPair[]
+single_fields = SkeletonCellFieldPair[]
+single_fields_dual = GenericCellField[]
+U = xh.fe_space
+nfields = length(U.spaces)
+cell_dofs_field_offsets = MultiField._get_cell_dofs_field_offsets(xh)
+
+for i in 1:nfields
+  view_range = cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
+  cell_values_field = lazy_map(a->view(a,view_range),cellu_dual)
+  cf_dual = CellField(U.spaces[i],cell_values_field)
+  scfp_plus = SkeletonCellFieldPair(cf_dual, xh[i])
+  scfp = SkeletonCellFieldPair(xh[i], xh[i])
+  push!(single_fields_dual,cf_dual)
+  push!(single_fields_plus,scfp_plus)
+  push!(single_fields,scfp)
+end
+
+xh_scfp = MultiFieldCellField(single_fields)
+uh_scfp, ph_scfp = xh_scfp
+
+xh_scfp_dual_plus = MultiFieldCellField(single_fields_plus)
+uh_scfp_dual_plus, ph_scfp_dual_plus = xh_scfp_dual_plus
+
+xh_dual = MultiFieldCellField(single_fields_dual)
+uh_dual, ph_dual = xh_dual
+
+g_Λ((uh,ph)) = ∫( mean(uh) + mean(ph) + mean(uh)*mean(ph) )dΛ
+a_Λ((uh,ph)) = ∫( - jump(uh*n_Λ)⊙mean(∇(ph))
+                  - mean(∇(uh))⊙jump(ph*n_Λ)
+                  + jump(uh*n_Λ)⊙jump(ph*n_Λ) )dΛ
+
+# few basic tests if integrations on trians are done fine
+@test sum(∫(uh_scfp*ph_scfp)*dΓ) == sum(∫(uh*ph)*dΓ)
+@test sum(∫(uh_scfp - ph_scfp)*dΓ) == sum(∫(uh - ph)*dΓ)
+@test sum(∫(uh_scfp*ph_scfp)*dΩ) == sum(∫(uh*ph)*dΩ)
+@test sum(∫(uh_scfp - ph_scfp)*dΩ) == sum(∫(uh - ph)*dΩ)
+
+# integrations involving SkeletonCellFieldPair derivative terms
+@test sum(∫(ph_scfp*∇(uh_scfp))*dΩ) == sum(∫(ph*∇(uh))*dΩ)
+@test sum(∫(uh_scfp*∇(ph_scfp))*dΩ) == sum(∫(uh*∇(ph))*dΩ)
+@test sum(∫(uh_scfp*∇(ph_scfp)⋅n_Γ)*dΓ) == sum(∫(uh*∇(ph)⋅n_Γ)*dΓ)
+@test sum(∫(ph_scfp*∇(uh_scfp)⋅n_Γ)*dΓ) == sum(∫(ph*∇(uh)⋅n_Γ)*dΓ)
+
+# dualized tests
+
+# few basic tests if integrations on trians are done fine
+@test sum(∫(uh_scfp_dual_plus*ph_scfp_dual_plus)*dΓ).value == sum(∫(uh*ph)*dΓ)
+@test sum(∫(uh_scfp_dual_plus*ph_scfp_dual_plus)*dΓ).partials == sum(∫(uh_dual*ph_dual)*dΓ).partials
+@test sum(∫(uh_scfp_dual_plus - ph_scfp_dual_plus)*dΓ).value == sum(∫(uh - ph)*dΓ)
+@test sum(∫(uh_scfp_dual_plus - ph_scfp_dual_plus)*dΓ).partials == sum(∫(uh_dual - ph_dual)*dΓ).partials
+@test sum(∫(uh_scfp_dual_plus*ph_scfp_dual_plus)*dΩ).value == sum(∫(uh*ph)*dΩ)
+@test sum(∫(uh_scfp_dual_plus*ph_scfp_dual_plus)*dΩ).partials == sum(∫(uh_dual*ph_dual)*dΩ).partials
+@test sum(∫(uh_scfp_dual_plus - ph_scfp_dual_plus)*dΩ).value == sum(∫(uh - ph)*dΩ)
+@test sum(∫(uh_scfp_dual_plus - ph_scfp_dual_plus)*dΩ).partials == sum(∫(uh_dual - ph_dual)*dΩ).partials
+
+# integrations involving SkeletonCellFieldPair derivative terms
+@test sum(∫(ph_scfp_dual_plus*∇(uh_scfp_dual_plus))*dΩ)[1].value == sum(∫(ph*∇(uh))*dΩ)[1]
+@test sum(∫(ph_scfp_dual_plus*∇(uh_scfp_dual_plus))*dΩ)[1].partials == sum(∫(ph_dual*∇(uh_dual))*dΩ)[1].partials
+@test sum(∫(uh_scfp_dual_plus*∇(ph_scfp_dual_plus))*dΩ)[1].value == sum(∫(uh*∇(ph))*dΩ)[1]
+@test sum(∫(uh_scfp_dual_plus*∇(ph_scfp_dual_plus))*dΩ)[1].partials == sum(∫(uh_dual*∇(ph_dual))*dΩ)[1].partials
+@test sum(∫(uh_scfp_dual_plus*∇(ph_scfp_dual_plus)⋅n_Γ)*dΓ).value == sum(∫(uh*∇(ph)⋅n_Γ)*dΓ)
+@test sum(∫(uh_scfp_dual_plus*∇(ph_scfp_dual_plus)⋅n_Γ)*dΓ).partials == sum(∫(uh_dual*∇(ph_dual)⋅n_Γ)*dΓ).partials
+@test sum(∫(ph_scfp_dual_plus*∇(uh_scfp_dual_plus)⋅n_Γ)*dΓ).value == sum(∫(ph*∇(uh)⋅n_Γ)*dΓ)
+@test sum(∫(ph_scfp_dual_plus*∇(uh_scfp_dual_plus)⋅n_Γ)*dΓ).partials == sum(∫(ph_dual*∇(uh_dual)⋅n_Γ)*dΓ).partials
+
+# test for SkeletonTriangulation
+@test sum(g_Λ((uh_scfp_dual_plus,ph_scfp_dual_plus))).value == sum(g_Λ((uh,ph)))
+@test sum(a_Λ((uh_scfp_dual_plus,ph_scfp_dual_plus))).value == sum(a_Λ((uh,ph)))
+
 end # module

test/MultiFieldTests/MultiFieldFEAutodiffTests.jl

@@ -3,12 +3,15 @@ module MultiFieldFEAutodiffTests
 using Test
 
 using Gridap.Algebra
+using Gridap.Arrays
 using Gridap.Geometry
 using Gridap.CellData
 using Gridap.ReferenceFEs
 using Gridap.FESpaces
 using Gridap.MultiField
 
+using ForwardDiff
+
 domain = (0,1,0,1)
 partition = (2,2)
 model = CartesianDiscreteModel(domain,partition)
@@ -132,4 +135,68 @@ c=jup_uh_ph[Ω]
 @test all(a .== map(x->x.array[1,1],c))
 @test all(b .== map(x->x.array[2,2],c))
 
+# AD tests for Multifiled functionals involving SkeletonTriangulation
+# compared with direct ForwardDiff results
+
+reffeV = ReferenceFE(lagrangian,Float64,2)
+reffeQ = ReferenceFE(lagrangian,Float64,1)
+V = TestFESpace(model,reffeV,conformity=:L2)
+Q = TestFESpace(model,reffeQ,conformity=:L2)
+
+U = TrialFESpace(V)
+P = TrialFESpace(Q)
+
+Y = MultiFieldFESpace([V, Q])
+X = MultiFieldFESpace([U, P])
+
+xh = FEFunction(X,rand(num_free_dofs(X)))
+
+Λ = SkeletonTriangulation(model)
+dΛ = Measure(Λ,2)
+n_Λ = get_normal_vector(Λ)
+
+g_Λ((uh,ph)) = ∫( mean(uh) + mean(ph) + mean(uh)*mean(ph) )dΛ
+# f_Λ((uh,ph)) = ∫( mean(uh*uh) + mean(uh*ph) + mean(ph*ph) )dΛ
+a_Λ((uh,ph)) = ∫( - jump(uh*n_Λ)⊙mean(∇(ph))
+                  - mean(∇(uh))⊙jump(ph*n_Λ)
+                  + jump(uh*n_Λ)⊙jump(ph*n_Λ) )dΛ
+
+function f_uh_free_dofs(f,xh,θ)
+  dir_u = similar(xh[1].dirichlet_values,eltype(θ))
+  dir_p = similar(xh[2].dirichlet_values,eltype(θ))
+  X = xh.fe_space
+  θ_uh = restrict_to_field(X,θ,1)
+  θ_ph = restrict_to_field(X,θ,2)
+  uh = FEFunction(X[1],θ_uh,dir_u)
+  ph = FEFunction(X[2],θ_ph,dir_p)
+  xh = MultiFieldFEFunction(θ,xh.fe_space,[uh,ph])
+  sum(f(xh))
+end
+# can also do the above by constructing a MultiFieldCellField
+
+g_Λ_(θ) = f_uh_free_dofs(g_Λ,xh,θ)
+# f_Λ_(θ) = f_uh_free_dofs(f_Λ,xh,θ)
+a_Λ_(θ) = f_uh_free_dofs(a_Λ,xh,θ)
+
+θ = get_free_dof_values(xh)
+
+# check if the evaluations are working
+# @test sum(f_Λ(xh)) == f_Λ_(θ)
+@test sum(a_Λ(xh)) == a_Λ_(θ)
+@test sum(g_Λ(xh)) == g_Λ_(θ)
+
+f_Ω((uh,ph)) = ∫(uh*uh + ph*uh + uh*ph)*dΩ
+f_Ω_(θ) = f_uh_free_dofs(f_Ω,xh,θ)
+gridapgradf_Ω = assemble_vector(gradient(f_Ω,xh),X)
+fdgradf_Ω = ForwardDiff.gradient(f_Ω_,θ)
+test_array(gridapgradf_Ω,fdgradf_Ω,≈)
+
+gridapgradg = assemble_vector(gradient(g_Λ,xh),X)
+fdgradg = ForwardDiff.gradient(g_Λ_,θ)
+test_array(gridapgradg,fdgradg,≈)
+
+gridapgrada = assemble_vector(gradient(a_Λ,xh),X)
+fdgrada = ForwardDiff.gradient(a_Λ_,θ)
+test_array(gridapgrada,fdgrada,≈)
+
 end # module