From 163a42d928973202ce03d423958d3bf5401af971 Mon Sep 17 00:00:00 2001 From: Kishore Nori Date: Fri, 17 Jun 2022 13:22:18 +1000 Subject: [PATCH 1/5] Added methods for Skeleton AD with operation inside mean and jump * This involved removing the change_domain method SkeletonCellFieldPair * Adding an evaluate! method for SkeletonCellFieldPair * Adding tests for the feature in FEAutodiff and also MultiField case --- src/CellData/SkeletonCellFieldPair.jl | 33 +++++++++++++++++-- test/FESpacesTests/FEAutodiffTests.jl | 24 ++++++++++++++ .../MultiFieldFEAutodiffTests.jl | 10 ++++-- 3 files changed, 61 insertions(+), 6 deletions(-) diff --git a/src/CellData/SkeletonCellFieldPair.jl b/src/CellData/SkeletonCellFieldPair.jl index e8c0f2f03..a5d6c2d5b 100644 --- a/src/CellData/SkeletonCellFieldPair.jl +++ b/src/CellData/SkeletonCellFieldPair.jl @@ -80,12 +80,39 @@ 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. +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) +# function change_domain(a::SkeletonCellFieldPair,target_trian::Triangulation,target_domain::DomainStyle) +# change_domain(a.plus,DomainStyle(a),target_trian,target_domain) +# end +#= +The above is not good because it converts a SkeletonCellFieldPair involved +in an operation inside mean or jump, for example the one resulting from the +expression mean(uh*uh), into its plus side, as the change_domain acts while +building the internal operation tree, in the example for * operator and this +gives wrong result as it chooses plus even for the case of Skeleton integration +=# + +# 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) giving right output +# or else it chooses the SkeletonCellFieldPair directly as the parent +# resulting in errors and it 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 +get_data(a::SkeletonCellFieldPair) = get_data(a.plus) + function Base.propertynames(a::SkeletonCellFieldPair, private::Bool=false) (fieldnames(typeof(a))...,:⁺,:plus,:⁻,:minus) end diff --git a/test/FESpacesTests/FEAutodiffTests.jl b/test/FESpacesTests/FEAutodiffTests.jl index 42dce637f..d6590605d 100644 --- a/test/FESpacesTests/FEAutodiffTests.jl +++ b/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 diff --git a/test/MultiFieldTests/MultiFieldFEAutodiffTests.jl b/test/MultiFieldTests/MultiFieldFEAutodiffTests.jl index 9ecc3e002..2e9dcdc86 100644 --- a/test/MultiFieldTests/MultiFieldFEAutodiffTests.jl +++ b/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,≈) From b1a266a1882039a1e356780c5fbbb407084be77e Mon Sep 17 00:00:00 2001 From: Kishore Nori Date: Fri, 17 Jun 2022 13:26:37 +1000 Subject: [PATCH 2/5] removed an old comment and commented code --- src/CellData/SkeletonCellFieldPair.jl | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/src/CellData/SkeletonCellFieldPair.jl b/src/CellData/SkeletonCellFieldPair.jl index a5d6c2d5b..002eb7895 100644 --- a/src/CellData/SkeletonCellFieldPair.jl +++ b/src/CellData/SkeletonCellFieldPair.jl @@ -68,10 +68,6 @@ 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 @@ -83,17 +79,6 @@ 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) -# end -#= -The above is not good because it converts a SkeletonCellFieldPair involved -in an operation inside mean or jump, for example the one resulting from the -expression mean(uh*uh), into its plus side, as the change_domain acts while -building the internal operation tree, in the example for * operator and this -gives wrong result as it chooses plus even for the case of Skeleton integration -=# - # 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) From 9ede7c4ec173a302e6dfcf91e49897178acd305c Mon Sep 17 00:00:00 2001 From: Kishore Nori Date: Fri, 17 Jun 2022 13:47:09 +1000 Subject: [PATCH 3/5] Added the fixes to NEWS.md linked with PR --- NEWS.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/NEWS.md b/NEWS.md index cea2cd20a..28db23962 100644 --- a/NEWS.md +++ b/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 From 553fcc0a69b41598c2717a0c9bba934f1eef84dc Mon Sep 17 00:00:00 2001 From: Kishore Nori Date: Fri, 17 Jun 2022 14:59:21 +1000 Subject: [PATCH 4/5] Improved comments for the choice of overloads with SkeletonCellFieldPair --- src/CellData/SkeletonCellFieldPair.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/CellData/SkeletonCellFieldPair.jl b/src/CellData/SkeletonCellFieldPair.jl index 002eb7895..45e641378 100644 --- a/src/CellData/SkeletonCellFieldPair.jl +++ b/src/CellData/SkeletonCellFieldPair.jl @@ -87,15 +87,23 @@ function evaluate!(cache,f::SkeletonCellFieldPair,x::CellPoint) lazy_map(evaluate,cell_field,cell_point) end -# fix for CellFieldAt{T}(parent::OperationCellField) giving right output -# or else it chooses the SkeletonCellFieldPair directly as the parent -# resulting in errors and it not the intended behaviour to have the parent -# as SkeletonCellFieldPair and is not consistent with our getproperty rules +#= +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 the ambiguity as the CellField at the SkeletonTrian +# for direct evaluations get_data(a::SkeletonCellFieldPair) = get_data(a.plus) function Base.propertynames(a::SkeletonCellFieldPair, private::Bool=false) From d125d519a104a8a037de14bf1f23decd667660d4 Mon Sep 17 00:00:00 2001 From: Kishore Nori Date: Fri, 17 Jun 2022 15:04:29 +1000 Subject: [PATCH 5/5] fixed and added more details in comments --- src/CellData/SkeletonCellFieldPair.jl | 50 +++++++++++++++++---------- 1 file changed, 31 insertions(+), 19 deletions(-) diff --git a/src/CellData/SkeletonCellFieldPair.jl b/src/CellData/SkeletonCellFieldPair.jl index 45e641378..0e0c3a3c8 100644 --- a/src/CellData/SkeletonCellFieldPair.jl +++ b/src/CellData/SkeletonCellFieldPair.jl @@ -68,16 +68,16 @@ function get_triangulation(a::SkeletonCellFieldPair) end #= -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 +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 =# # If SkeletonCellFieldPair is evaluated we just pick the plus side parent function evaluate!(cache,f::SkeletonCellFieldPair,x::CellPoint) @@ -88,22 +88,34 @@ function evaluate!(cache,f::SkeletonCellFieldPair,x::CellPoint) 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: +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. +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 the ambiguity as the CellField at the SkeletonTrian -# for direct evaluations +#= +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)