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

Changed the method for `SkeletonCellFieldPair` for AD to work with operations inside `mean` and `jump`

NEWS.md

@@ -10,6 +10,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - 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)
 
+### Fixed
+- Fixed the behavior of `gradient` for functionals involving operations of `CellFields` inside `mean` and `jump` terms of Skeleton Integration terms. Since PR [#800](https://github.com/gridap/Gridap.jl/pull/800)
+- Fixed the behavior `SkeletonCellFieldPair` at the Boundary integration terms. Since PR [#800](https://github.com/gridap/Gridap.jl/pull/800)
+
 ## [0.17.13] - 2022-05-31
 
 ### Added

src/CellData/SkeletonCellFieldPair.jl

@@ -68,24 +68,56 @@ function get_triangulation(a::SkeletonCellFieldPair)
 end
 
 #=
-The below change_domain has been overloaded for SkeletonCellFieldPair
-to work without errors for integrations over other triangulation - Ω
-(BodyFittedTriangulation) and Γ (BoundaryTriangulation).
-
-We arbitrarily choose plus side of SkeletonCellFieldPair for evaluations, as
-we don't use the DomainContributions associated with Ω and Γ while performing
-the sensitivities for SkeletonTriangulation (Λ) DomainContribution. This
-fix is for error free evaluation of the functional when a SkeletonCellFieldPair
-is passed into it. Ideally, if we could parse and extract only the Skeleton
-integration terms from the functional's julia function form, this fix is not
-required, but this is not trivial to do. On the positive side, since the
-evaluations are all lazy and not used, this doesn't put any noticable memory
-or computational overhead.
+We arbitrarily choose plus side of SkeletonCellFieldPair for evaluations, as we
+don't use the DomainContributions associated with Ω and Γ while performing the
+sensitivities for SkeletonTriangulation (Λ) DomainContribution. This fix is for
+error free evaluation of the functional when a SkeletonCellFieldPair is passed
+into it. Ideally, if we could parse and extract only the Skeleton integration
+terms from the functional's julia function form, this fix is not required, but
+this is not trivial to do. On the positive side, since the evaluations are all
+lazy and not used, this doesn't put any noticable memory or computational
+overhead. Ofcourse, it is made sure that the such plus side pick doesn't happen
+when the integration over the SkeletonTriangulation
 =#
-function change_domain(a::SkeletonCellFieldPair,target_trian::Triangulation,target_domain::DomainStyle)
-  change_domain(a.plus,DomainStyle(a),target_trian,target_domain)
+# If SkeletonCellFieldPair is evaluated we just pick the plus side parent
+function evaluate!(cache,f::SkeletonCellFieldPair,x::CellPoint)
+  _f, _x = _to_common_domain(f.plus.parent,x)
+  cell_field = get_data(_f)
+  cell_point = get_data(_x)
+  lazy_map(evaluate,cell_field,cell_point)
 end
 
+#=
+Fix for CellFieldAt{T}(parent::OperationCellField) for OperationCellField
+involving args which are not FEFunctions. It creates a problem with
+SkeletonCellFieldPair giving the wrong output or errors as it chooses the
+SkeletonCellFieldPair directly as the parent rather the parent CellField of the
+side of choice. In code terms it is the following:
+
+CellFieldAt{T}(OperationCellField) results in calling CellFieldAt{T}(operands).
+Without this override, CellFieldAt{T}(a::SkeletonCellField) results in the
+CellFieldAt structure with parent as SkeletonCellField and not a.T side
+CellField, which is not correct! This resulting in errors and is not the
+intended behaviour to have the parent as SkeletonCellFieldPair and is not
+consistent with our getproperty rules.
+=#
+function CellFieldAt{T}(parent::SkeletonCellFieldPair) where T
+  getproperty(parent,T)
+end
+
+#=
+To handle the evaluation of SkeletonCellFieldPair at BoundaryTriangulation
+making it consistent with plus side choice of direct evaluation of SCFP in
+general this is the case when trian of CellField is not the same as that of the
+quadrature.
+But we are protected from inconsistent behaviour at SkeletonTriangulation as it
+fails due to similar ambiguity as the CellField at the SkeletonTrian for direct
+evaluations, so get_data doesn't create any side effects.
+In case of BodyFittedTriangulation there is no hit to get_data directly and
+evaluate! handles it as intended
+=#
+get_data(a::SkeletonCellFieldPair) = get_data(a.plus)
+
 function Base.propertynames(a::SkeletonCellFieldPair, private::Bool=false)
   (fieldnames(typeof(a))...,:⁺,:plus,:⁻,:minus)
 end

test/FESpacesTests/FEAutodiffTests.jl

@@ -1,5 +1,6 @@
 module FEAutodiffTests
 
+using Test
 using LinearAlgebra
 using Gridap.Algebra
 using Gridap.FESpaces
@@ -120,6 +121,13 @@ a_Λ(u) = ∫( - jump(u*n_Λ)⊙mean(∇(u))
             - mean(∇(u))⊙jump(u*n_Λ)
             + jump(u*n_Λ)⊙jump(u*n_Λ) )dΛ
 
+# functionals having mean and jump of products of FEFunctions/CellFields
+g_Λ(uh) = ∫(mean(uh*uh))*dΛ
+h_Λ(uh) = ∫(jump(uh*uh*n_Λ)⋅jump(uh*uh*n_Λ))*dΛ
+j_Λ(uh) = ∫(mean(∇(uh)*uh)⋅jump((∇(uh)⋅∇(uh))*n_Λ))*dΛ
+
+@test sum(∫(mean(uh*uh))*dΛ) == 0.5*sum(∫(uh.plus*uh.plus + uh.minus*uh.minus)*dΛ)
+
 function f_uh_free_dofs(f,uh,θ)
   dir = similar(uh.dirichlet_values,eltype(θ))
   uh = FEFunction(U,θ,dir)
@@ -138,4 +146,20 @@ gridapgrada = assemble_vector(gradient(a_Λ,uh),U)
 fdgrada = ForwardDiff.gradient(a_Λ_,θ)
 test_array(gridapgrada,fdgrada,≈)
 
+g_Λ_(θ) = f_uh_free_dofs(g_Λ,uh,θ)
+h_Λ_(θ) = f_uh_free_dofs(h_Λ,uh,θ)
+j_Λ_(θ) = f_uh_free_dofs(j_Λ,uh,θ)
+
+gridapgradg = assemble_vector(gradient(g_Λ,uh),U)
+fdgradg = ForwardDiff.gradient(g_Λ_,θ)
+test_array(gridapgradg,fdgradg,≈)
+
+gridapgradh = assemble_vector(gradient(h_Λ,uh),U)
+fdgradh = ForwardDiff.gradient(h_Λ_,θ)
+test_array(gridapgradh,fdgradh,≈)
+
+gridapgradj = assemble_vector(gradient(j_Λ,uh),U)
+fdgradj = ForwardDiff.gradient(j_Λ_,θ)
+test_array(gridapgradj,fdgradj,≈)
+
 end # module

test/MultiFieldTests/MultiFieldFEAutodiffTests.jl

@@ -156,7 +156,7 @@ 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Λ
+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Λ
@@ -175,13 +175,13 @@ end
 # can also do the above by constructing a MultiFieldCellField
 
 g_Λ_(θ) = f_uh_free_dofs(g_Λ,xh,θ)
-# f_Λ_(θ) = f_uh_free_dofs(f_Λ,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(f_Λ(xh)) == f_Λ_(θ)
 @test sum(a_Λ(xh)) == a_Λ_(θ)
 @test sum(g_Λ(xh)) == g_Λ_(θ)
 
@@ -191,6 +191,10 @@ gridapgradf_Ω = assemble_vector(gradient(f_Ω,xh),X)
 fdgradf_Ω = ForwardDiff.gradient(f_Ω_,θ)
 test_array(gridapgradf_Ω,fdgradf_Ω,≈)
 
+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,≈)