From a0c9d840cee08ca28e8605165492dc48c067f53f Mon Sep 17 00:00:00 2001 From: Kishore Nori Date: Tue, 21 Jun 2022 19:37:19 +1000 Subject: [PATCH 1/3] Extended the `jacobian` to functionals involving Skeleton terms * Added corresponding tests compared with ForwardDiff * The high-level api remains the same --- src/FESpaces/FEAutodiff.jl | 59 +++++++++++++++++++---- test/FESpacesTests/FEAutodiffTests.jl | 68 ++++++++++++++++++++++++++- 2 files changed, 117 insertions(+), 10 deletions(-) diff --git a/src/FESpaces/FEAutodiff.jl b/src/FESpaces/FEAutodiff.jl index 132dd50db..141c668e5 100644 --- a/src/FESpaces/FEAutodiff.jl +++ b/src/FESpaces/FEAutodiff.jl @@ -165,9 +165,7 @@ end Gridap, where we can use or import required functionalities from other modules without any circular dependency =# -function autodiff_array_gradient( - a, i_to_x, j_to_i::SkeletonPair) - +function autodiff_array_gradient(a, i_to_x, j_to_i::SkeletonPair) i_to_xdual = lazy_map(DualizeMap(ForwardDiff.gradient),i_to_x) # dual output of both sides at once @@ -191,12 +189,57 @@ function autodiff_array_gradient( lazy_map(BlockMap(2,[1,2]),j_to_result_plus,j_to_result_minus) end -function autodiff_array_jacobian( - a,i_to_x,j_to_i::SkeletonPair) - @notimplemented +function autodiff_array_jacobian(a,i_to_x,j_to_i::SkeletonPair) + i_to_xdual = lazy_map(DualizeMap(ForwardDiff.jacobian),i_to_x) + j_to_ydual_plus, j_to_ydual_minus = a(i_to_xdual) + + # bilinear form tested with test basis functions dv results in + # DomainContribution containing vector of vector blocks of 2 vectors + # coming from plus side and minus side of dv + # so we densify each of them into plain vectors and construct + # construct back the jacobian contributions into a MatrixBlock + + densify = DensifyInnerMostBlockLevelMap() + j_to_ydual_plus_dense = lazy_map(densify, j_to_ydual_plus) + j_to_ydual_minus_dense = lazy_map(densify, j_to_ydual_minus) + + # Work for plus side + j_to_x_plus = lazy_map(Reindex(i_to_x),j_to_i.plus) + j_to_cfg_plus = lazy_map(ConfigMap(ForwardDiff.jacobian),j_to_x_plus) + j_to_result_plus_dense = lazy_map(AutoDiffMap(ForwardDiff.jacobian), + j_to_ydual_plus_dense,j_to_x_plus,j_to_cfg_plus) + + # Work for minus side + j_to_x_minus = lazy_map(Reindex(i_to_x),j_to_i.minus) + j_to_cfg_minus = lazy_map(ConfigMap(ForwardDiff.jacobian),j_to_x_minus) + j_to_result_minus_dense = lazy_map(AutoDiffMap(ForwardDiff.jacobian), + j_to_ydual_minus_dense,j_to_x_minus,j_to_cfg_minus) + + # j_to_result_plus_dense/j_to_result_minus_dense can be (and must be) + # laid out into 2x2 block matrices + num_dofs_x_cell = lazy_map(length,i_to_x) + num_dofs_x_face_and_cell_plus = lazy_map(Reindex(num_dofs_x_cell),j_to_i.plus) + num_dofs_x_face_and_cell_minus = lazy_map(Reindex(num_dofs_x_cell),j_to_i.minus) + + J_11 = lazy_map((x,b)->view(x, 1:b,:), + j_to_result_plus_dense,num_dofs_x_face_and_cell_plus) + + J_21 = lazy_map((x,b)->view(x, b+1:size(x,1),:), + j_to_result_plus_dense,num_dofs_x_face_and_cell_minus) + + J_12 = lazy_map((x,b)->view(x, 1:b,:), + j_to_result_minus_dense,num_dofs_x_face_and_cell_plus) + + J_22 = lazy_map((x,b)->view(x, b+1:size(x,1),:), + j_to_result_minus_dense,num_dofs_x_face_and_cell_minus) + + # Assembly on SkeletonTriangulation expects an array of interior of facets + # where each entry is a 2x2-block MatrixBlock with the blocks of the Jacobian + # matrix + bm_jacobian = BlockMap((2,2),[(1,1),(2,1),(1,2),(2,2)]) + lazy_map(bm_jacobian, J_11, J_21, J_12, J_22) end -function autodiff_array_hessian( - a,i_to_x,i_to_j::SkeletonPair) +function autodiff_array_hessian(a,i_to_x,i_to_j::SkeletonPair) @notimplemented end diff --git a/test/FESpacesTests/FEAutodiffTests.jl b/test/FESpacesTests/FEAutodiffTests.jl index d6590605d..727fb197e 100644 --- a/test/FESpacesTests/FEAutodiffTests.jl +++ b/test/FESpacesTests/FEAutodiffTests.jl @@ -11,6 +11,7 @@ using Gridap.TensorValues using Gridap.CellData using Gridap.ReferenceFEs using ForwardDiff +using SparseArrays domain = (0,1,0,1) partition = (2,2) @@ -94,7 +95,7 @@ cell_j_auto = get_array(jacobian(u->res(u,dv),uh)) test_array(cell_j_auto,cell_j,≈) -# comparing AD of integration over Skeleton faces with ForwardDiff results +## comparing AD of integration over Skeleton faces with ForwardDiff results ## model = CartesianDiscreteModel((0.,1.,0.,1.),(3,3)) Ω = Triangulation(model) @@ -118,7 +119,7 @@ uh = FEFunction(U,rand(num_free_dofs(U))) f_Λ(uh) = ∫(mean(uh))*dΛ a_Λ(u) = ∫( - jump(u*n_Λ)⊙mean(∇(u)) - - mean(∇(u))⊙jump(u*n_Λ) + - mean(Δ(u)) + jump(u*n_Λ)⊙jump(u*n_Λ) )dΛ # functionals having mean and jump of products of FEFunctions/CellFields @@ -138,6 +139,8 @@ f_Λ_(θ) = f_uh_free_dofs(f_Λ,uh,θ) a_Λ_(θ) = f_uh_free_dofs(a_Λ,uh,θ) θ = get_free_dof_values(uh) +# testing gradients of functionals involving Skeleton integration terms + gridapgradf = assemble_vector(gradient(f_Λ,uh),U) fdgradf = ForwardDiff.gradient(f_Λ_,θ) test_array(gridapgradf,fdgradf,≈) @@ -162,4 +165,65 @@ gridapgradj = assemble_vector(gradient(j_Λ,uh),U) fdgradj = ForwardDiff.gradient(j_Λ_,θ) test_array(gridapgradj,fdgradj,≈) +# testing jacobians of functionals involving Skeleton inetgration terms +uh = FEFunction(U,rand(num_free_dofs(U))) +θ = get_free_dof_values(uh) +dv = get_fe_basis(V) +du = get_trial_fe_basis(U) + +a(uh,vh) = ∫(mean(uh)*mean(vh))*dΛ +b(uh,vh) = ∫(-jump(uh*n_Λ)⊙mean(∇(vh)) + + jump(vh*n_Λ)⊙mean(∇(uh)) + + mean(Δ(uh))*mean(Δ(vh)) + + jump(uh*n_Λ)⊙jump(vh*n_Λ))*dΛ +g(uh,vh) = ∫(mean(uh*vh))*dΛ +h(uh,vh) = ∫(jump(uh*vh*n_Λ)⊙jump(vh*vh*n_Λ))*dΛ +j(uh,vh) = ∫(mean(∇(vh)*uh)⊙jump((∇(vh)⊙∇(uh))*n_Λ))*dΛ + +function _change_input(f,θ,uh) + dir = similar(uh.dirichlet_values,eltype(θ)) + U = uh.fe_space + uh = FEFunction(U,θ,dir) + T = eltype(θ) + matrix_type = SparseMatrixCSC{T,Int} + vector_type = Vector{T} + assem = SparseMatrixAssembler(matrix_type,vector_type,U,U) + assemble_vector(f(uh),assem,U) +end + +f(uh) = a(uh,dv) +jac_gridap_a = assemble_matrix(jacobian(f,uh),U,V) +collect(get_array(jacobian(f,uh))) # just to check the working +f_(θ) = _change_input(f,θ,uh) +jac_forwdiff_a = ForwardDiff.jacobian(f_,θ) +test_array(jac_gridap_a,jac_forwdiff_a,≈) + +f(uh) = b(uh,dv) +jac_gridap_b = assemble_matrix(jacobian(f,uh),U,V) +collect(get_array(jacobian(f,uh))) # just to check the working +f_(θ) = _change_input(f,θ,uh) +jac_forwdiff_b = ForwardDiff.jacobian(f_,θ) +test_array(jac_gridap_b,jac_forwdiff_b,≈) + +f(uh) = g(uh,dv) +jac_gridap_g = assemble_matrix(jacobian(f,uh),U,V) +collect(get_array(jacobian(f,uh))) # just to check the working +f_(θ) = _change_input(f,θ,uh) +jac_forwdiff_g = ForwardDiff.jacobian(f_,θ) +test_array(jac_gridap_g,jac_forwdiff_g,≈) + +f(uh) = h(uh,dv) +jac_gridap_h = assemble_matrix(jacobian(f,uh),U,V) +collect(get_array(jacobian(f,uh))) # just to check the working +f_(θ) = _change_input(f,θ,uh) +jac_forwdiff_h = ForwardDiff.jacobian(f_,θ) +test_array(jac_gridap_h,jac_forwdiff_h,≈) + +f(uh) = j(uh,dv) +jac_gridap_j = assemble_matrix(jacobian(f,uh),U,V) +collect(get_array(jacobian(f,uh))) # just to check the working +f_(θ) = _change_input(f,θ,uh) +jac_forwdiff_j = ForwardDiff.jacobian(f_,θ) +test_array(jac_gridap_j,jac_forwdiff_j,≈) + end # module From 6a2966072c39a81f70ff23d17179e3fa61c3a9d2 Mon Sep 17 00:00:00 2001 From: Kishore Nori Date: Tue, 21 Jun 2022 23:36:30 +1000 Subject: [PATCH 2/3] Corrected earlier comments and added some more details --- NEWS.md | 3 ++- src/FESpaces/FEAutodiff.jl | 15 +++++++-------- test/FESpacesTests/FEAutodiffTests.jl | 10 ++++++++++ 3 files changed, 19 insertions(+), 9 deletions(-) diff --git a/NEWS.md b/NEWS.md index 28db23962..1bad95512 100644 --- a/NEWS.md +++ b/NEWS.md @@ -7,8 +7,9 @@ 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 of `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) +- Added functionality to take jacobian of functional involving integration (`DomainContribution`) over Skeleton faces (obtained from testing bilinear form with the whole set of test `fe_basis`), with respect to the degrees-of-freedom of `FEFunction`. The interface remains the same - `jacobian(f,uh)`. Since PR [#803](https://github.com/gridap/Gridap.jl/pull/803) ### 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) diff --git a/src/FESpaces/FEAutodiff.jl b/src/FESpaces/FEAutodiff.jl index 141c668e5..caff6c448 100644 --- a/src/FESpaces/FEAutodiff.jl +++ b/src/FESpaces/FEAutodiff.jl @@ -193,11 +193,11 @@ function autodiff_array_jacobian(a,i_to_x,j_to_i::SkeletonPair) i_to_xdual = lazy_map(DualizeMap(ForwardDiff.jacobian),i_to_x) j_to_ydual_plus, j_to_ydual_minus = a(i_to_xdual) - # bilinear form tested with test basis functions dv results in - # DomainContribution containing vector of vector blocks of 2 vectors - # coming from plus side and minus side of dv - # so we densify each of them into plain vectors and construct - # construct back the jacobian contributions into a MatrixBlock + # Bilinear form when tested with test basis functions dv results in + # DomainContribution containing vector of 2-block `VectorBlock{Vector}` + # each block coming from plus side and minus side of dv. + # So we densify each of the `VectorBlock`` into plain vectors and construct + # back the jacobian contributions into a MatrixBlock densify = DensifyInnerMostBlockLevelMap() j_to_ydual_plus_dense = lazy_map(densify, j_to_ydual_plus) @@ -233,9 +233,8 @@ function autodiff_array_jacobian(a,i_to_x,j_to_i::SkeletonPair) J_22 = lazy_map((x,b)->view(x, b+1:size(x,1),:), j_to_result_minus_dense,num_dofs_x_face_and_cell_minus) - # Assembly on SkeletonTriangulation expects an array of interior of facets - # where each entry is a 2x2-block MatrixBlock with the blocks of the Jacobian - # matrix + # Assembly on SkeletonTriangulation expects an array of facets where each + # entry is a 2x2-block MatrixBlock with the blocks of the Jacobian matrix bm_jacobian = BlockMap((2,2),[(1,1),(2,1),(1,2),(2,2)]) lazy_map(bm_jacobian, J_11, J_21, J_12, J_22) end diff --git a/test/FESpacesTests/FEAutodiffTests.jl b/test/FESpacesTests/FEAutodiffTests.jl index 727fb197e..03e678a67 100644 --- a/test/FESpacesTests/FEAutodiffTests.jl +++ b/test/FESpacesTests/FEAutodiffTests.jl @@ -180,6 +180,16 @@ g(uh,vh) = ∫(mean(uh*vh))*dΛ h(uh,vh) = ∫(jump(uh*vh*n_Λ)⊙jump(vh*vh*n_Λ))*dΛ j(uh,vh) = ∫(mean(∇(vh)*uh)⊙jump((∇(vh)⊙∇(uh))*n_Λ))*dΛ +# We use the strongly-typed lowel-level interface of SparseMatrixCSC{T,Int} here +# as ForwardDiff doesn't work directly through `assemble_vector``, this is due +# to the typing of SparseMatrixCSC{T,Int} and Vector{T} inside high-level +# `assemble_vector` API (via the `SparseMatrixAssembler`) using the type T of +# the eltype of FESpace U, which remains Float even though Dual numbers are +# passed into θ. +# So as to mitigate this problem, low-level interface of assemble_vector is +# being used with SparseMatrixCSC{T,Int} and Vector{T} constructed by hand with +# types of θ which can be dual numbers when using ForwardDiff or even # +# ReverseDiff types, when using ReverseDiff! This is just for testing purposes function _change_input(f,θ,uh) dir = similar(uh.dirichlet_values,eltype(θ)) U = uh.fe_space From 8b7d3bacad32e4f5618bde57686d29aaf4653e68 Mon Sep 17 00:00:00 2001 From: Kishore Nori Date: Wed, 22 Jun 2022 13:03:32 +1000 Subject: [PATCH 3/3] Made changes to skeleton `FEAutodiffTests.jl` to reduce compute time * removed higher order derivative tests for jacobian with skeleton functionals * reduced the dimension of problem to 2x2 and order = 1 for skeleton tests * wrapped repeated test code into a function to reduce compile time --- test/FESpacesTests/FEAutodiffTests.jl | 58 +++++++++++++-------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/test/FESpacesTests/FEAutodiffTests.jl b/test/FESpacesTests/FEAutodiffTests.jl index 03e678a67..2566ac4b3 100644 --- a/test/FESpacesTests/FEAutodiffTests.jl +++ b/test/FESpacesTests/FEAutodiffTests.jl @@ -97,7 +97,7 @@ test_array(cell_j_auto,cell_j,≈) ## comparing AD of integration over Skeleton faces with ForwardDiff results ## -model = CartesianDiscreteModel((0.,1.,0.,1.),(3,3)) +model = CartesianDiscreteModel((0.,1.,0.,1.),(2,2)) Ω = Triangulation(model) Γ = BoundaryTriangulation(model) Λ = SkeletonTriangulation(model) @@ -109,7 +109,7 @@ dΛ = Measure(Λ,2) n_Γ = get_normal_vector(Γ) n_Λ = get_normal_vector(Λ) -reffe = ReferenceFE(lagrangian,Float64,2) +reffe = ReferenceFE(lagrangian,Float64,1) V = TestFESpace(model,reffe,conformity=:L2) u(x) = sin(norm(x)) @@ -172,13 +172,13 @@ dv = get_fe_basis(V) du = get_trial_fe_basis(U) a(uh,vh) = ∫(mean(uh)*mean(vh))*dΛ -b(uh,vh) = ∫(-jump(uh*n_Λ)⊙mean(∇(vh)) + - jump(vh*n_Λ)⊙mean(∇(uh)) + - mean(Δ(uh))*mean(Δ(vh)) + - jump(uh*n_Λ)⊙jump(vh*n_Λ))*dΛ -g(uh,vh) = ∫(mean(uh*vh))*dΛ +# b(uh,vh) = ∫(-jump(uh*n_Λ)⊙mean(∇(vh)) + +# jump(vh*n_Λ)⊙mean(∇(uh)) + +# mean(Δ(uh))*mean(Δ(vh)) + +# jump(uh*n_Λ)⊙jump(vh*n_Λ))*dΛ +# g(uh,vh) = ∫(mean(uh*vh))*dΛ h(uh,vh) = ∫(jump(uh*vh*n_Λ)⊙jump(vh*vh*n_Λ))*dΛ -j(uh,vh) = ∫(mean(∇(vh)*uh)⊙jump((∇(vh)⊙∇(uh))*n_Λ))*dΛ +# j(uh,vh) = ∫(mean(∇(vh)*uh)⊙jump((∇(vh)⊙∇(uh))*n_Λ))*dΛ # We use the strongly-typed lowel-level interface of SparseMatrixCSC{T,Int} here # as ForwardDiff doesn't work directly through `assemble_vector``, this is due @@ -201,39 +201,39 @@ function _change_input(f,θ,uh) assemble_vector(f(uh),assem,U) end +function _assemble_jacobian(f,uh) + assemble_matrix(jacobian(f,uh),U,V) +end + f(uh) = a(uh,dv) -jac_gridap_a = assemble_matrix(jacobian(f,uh),U,V) +jac_gridap_a = _assemble_jacobian(f,uh) collect(get_array(jacobian(f,uh))) # just to check the working f_(θ) = _change_input(f,θ,uh) jac_forwdiff_a = ForwardDiff.jacobian(f_,θ) test_array(jac_gridap_a,jac_forwdiff_a,≈) -f(uh) = b(uh,dv) -jac_gridap_b = assemble_matrix(jacobian(f,uh),U,V) -collect(get_array(jacobian(f,uh))) # just to check the working -f_(θ) = _change_input(f,θ,uh) -jac_forwdiff_b = ForwardDiff.jacobian(f_,θ) -test_array(jac_gridap_b,jac_forwdiff_b,≈) +# f(uh) = b(uh,dv) +# jac_gridap_b = _assemble_jacobian(f,uh) +# f_(θ) = _change_input(f,θ,uh) +# jac_forwdiff_b = ForwardDiff.jacobian(f_,θ) +# test_array(jac_gridap_b,jac_forwdiff_b,≈) -f(uh) = g(uh,dv) -jac_gridap_g = assemble_matrix(jacobian(f,uh),U,V) -collect(get_array(jacobian(f,uh))) # just to check the working -f_(θ) = _change_input(f,θ,uh) -jac_forwdiff_g = ForwardDiff.jacobian(f_,θ) -test_array(jac_gridap_g,jac_forwdiff_g,≈) +# f(uh) = g(uh,dv) +# jac_gridap_g = _assemble_jacobian(f,uh) +# f_(θ) = _change_input(f,θ,uh) +# jac_forwdiff_g = ForwardDiff.jacobian(f_,θ) +# test_array(jac_gridap_g,jac_forwdiff_g,≈) f(uh) = h(uh,dv) -jac_gridap_h = assemble_matrix(jacobian(f,uh),U,V) -collect(get_array(jacobian(f,uh))) # just to check the working +jac_gridap_h = _assemble_jacobian(f,uh) f_(θ) = _change_input(f,θ,uh) jac_forwdiff_h = ForwardDiff.jacobian(f_,θ) test_array(jac_gridap_h,jac_forwdiff_h,≈) -f(uh) = j(uh,dv) -jac_gridap_j = assemble_matrix(jacobian(f,uh),U,V) -collect(get_array(jacobian(f,uh))) # just to check the working -f_(θ) = _change_input(f,θ,uh) -jac_forwdiff_j = ForwardDiff.jacobian(f_,θ) -test_array(jac_gridap_j,jac_forwdiff_j,≈) +# f(uh) = j(uh,dv) +# jac_gridap_j = _assemble_jacobian(f,uh) +# f_(θ) = _change_input(f,θ,uh) +# jac_forwdiff_j = ForwardDiff.jacobian(f_,θ) +# test_array(jac_gridap_j,jac_forwdiff_j,≈) end # module