gridap / Gridap.jl

Grid-based approximation of partial differential equations in Julia
MIT License
710 stars 97 forks source link

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

Open kishore-nori opened 1 year ago

kishore-nori commented 1 year ago

The Laplacian of product of CellFields (formed using the same FESpace and Triangulation) fails. I came across the following, https://github.com/gridap/Gridap.jl/blob/672e3d9027246bfa677c07b4e9fabb9e1fde9eb4/src/Fields/FieldsInterfaces.jl#L94-L102 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)

carlodev commented 1 year ago

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 commented 1 year ago

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.