Laplacian of product of `CellField`s doesn't work. #875

kishore-nori · 2023-02-12T23:04:11Z

The Laplacian of product of CellFields (formed using the same FESpace and Triangulation) fails. I came across the following,

Gridap.jl/src/Fields/FieldsInterfaces.jl

Lines 94 to 102 in 672e3d9

    
           function push_∇∇(∇∇a::Field,ϕ::Field) 
        
             @notimplemented """\n 
        
             Second order derivatives of quantities defined in the reference domain not implemented yet. 
        
             This is a feature that we want to have at some point in Gridap. 
        
             If you are ready to help with this implementation, please contact the 
        
             Gridap administrators. 
        
             """ 
        
           end

But the error doesn't come from here, the evaluation doesn't hit the above, at least from preliminary debugging observation. The following is the MWE:

using Gridap
using Test

domain = (0.,1.,0.,1.)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

g(x) = exp(-norm(x))
reffe = ReferenceFE(lagrangian,Float64,2)
V = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary")
U = TrialFESpace(V,g)

n = num_free_dofs(U)
uh = FEFunction(U,rand(n))
qh = FEFunction(U,rand(n))

p = Point(0.1,0.2)

cf1 = ∇(uh + qh)
@show @test cf1(p) ≈ ∇(uh)(p) + ∇(qh)(p) # works 

cf2 = ∇(uh*qh)
@show @test cf2(p) ≈ (qh*∇(uh))(p) + (uh*∇(qh))(p) # works 

cf3 = Δ(uh + qh) # equiv to tr(∇∇(⋅)), making it an OperationCellField
@show @test cf3(p) ≈ Δ(uh)(p) + Δ(qh)(p) # works 

cf4 = ∇∇(uh + qh)
@show @test cf4(p) ≈ ∇∇(uh)(p) + ∇∇(qh)(p) # works 

# Although the CellField is built, when evaluated it results in error
cf5 = ∇∇(uh*qh) # works 
cf5(p) # fails 

cf6 = Δ(uh*qh) # even the construction fails

I would be happy to fix the related issue and implement the second order derivatives based on your guidance and suggestions, as the @notimplemented message directs.

cc @amartinhuertas, @oriolcg (I think latest commits along these lines where your contributions)

The text was updated successfully, but these errors were encountered:

carlodev · 2023-03-09T11:32:25Z

I don't know if this is completely correct. I looked also as push_∇(∇a::Field,ϕ::Field) = pinvJt(∇(ϕ))⋅∇a:

function push_∇∇(∇∇a::Field,ϕ::Field)
  s = pinvJt(∇(ϕ))
  tr(s⋅s⋅∇∇a)
end

kishore-nori · 2023-03-10T03:04:17Z

Hi @carlodev, thank you for the possible solution. I'll work it out on paper and get back. At first sight it seems to be alright to me, although I am not sure on whether the above is only for a ϕ which is affine linear. An other thing is, like I mentioned above, the error doesn't seem to be from the unavailability of above, at least from my very preliminary debugging experience, rather, it seems to me that the higher order chain rules are unavailable as well. I am in middle of something else right now, will start working on it after 10 days and get back here, sorry. For now, until the above is fully resolved, we will have to rely on manual product rule expansion of such operators.

kishore-nori self-assigned this Mar 10, 2023

kishore-nori linked a pull request Nov 17, 2024 that will close this issue

Adding rules and methods for hessian of sum of scalar CellFields and product of two scalar CellFields #1053

Draft

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Laplacian of product of `CellField`s doesn't work. #875

Laplacian of product of `CellField`s doesn't work. #875

kishore-nori commented Feb 12, 2023 •

edited

Loading

carlodev commented Mar 9, 2023

kishore-nori commented Mar 10, 2023

Laplacian of product of CellFields doesn't work. #875

Laplacian of product of CellFields doesn't work. #875

Comments

kishore-nori commented Feb 12, 2023 • edited Loading

carlodev commented Mar 9, 2023

kishore-nori commented Mar 10, 2023

Laplacian of product of `CellField`s doesn't work. #875

Laplacian of product of `CellField`s doesn't work. #875

kishore-nori commented Feb 12, 2023 •

edited

Loading