Heat equation with homogeneous Neumann boundary condition

[urbainvaes]

Hi,

Together with a colleague, we have been trying for a while to solve the heat equation with homogeneous Neumann boundary conditions and a constant (in space and time) heat source.

With Dirichlet boundary conditions, the following code works well, but with homogeneous Neumann boundary condition, it does not work at all. Here is our code:

using Gridap
domain = (-1, +1, -1, +1)
partition = (20, 20)
model = CartesianDiscreteModel(domain, partition)
order = 1
reffe = ReferenceFE(lagrangian, Float64, order)
# V0 = TestFESpace(model, reffe, dirichlet_tags="boundary")   # Dirichlet boundary, this works perfectly
V0 = TestFESpace(model, reffe, dirichlet_tags=String[])   # Neumann boundary, does not work
u₀ = x -> 0.
g = x -> 0.
Ug = TrialFESpace(V0, g)

degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)

f = x -> 1.
m(t, dtu, v) = ∫(v * dtu)dΩ
a(t, u, v) = ∫(∇(v) ⋅ ∇(u))dΩ
l(t, v) = ∫(v * f)dΩ
op = TransientLinearFEOperator((m, a), l, Ug, V0, constant_forms=(true, true))

Δt = 0.05
θ = 0.5
solver = ThetaMethod(LUSolver(), Δt, θ)

t0, tF = 0.0, 10.0
uh0 = interpolate_everywhere(u₀, Ug)
uh = solve(solver, op, t0, tF, uh0)

To convince oneself that this does not work, one can print the evolution of the average temperature (which should increase linearly):

quad = CellQuadrature(Ω, 2)
createpvd("results") do pvd
  pvd[0] = createvtk(Ω, "$outdir/parabolic/results_0" * ".vtu", cellfields=["u" => uh0])
  for (i, dn) in enumerate(uh)
    tn, uhn = dn
    println(sum(integrate(uhn, quad)) / π)
    pvd[tn] = createvtk(Ω, "out/results_$i" * ".vtu", cellfields=["u" => uhn])
  end
end

Are we doing something wrong? Any help would be greatly appreciated.

[Antoinemarteau]

The first issue I see is that you don't provide a TransientFESpace to the TransientLinearFEOperator (see tuto 17), like:

g(t) = x -> 0
Ug = TransientTrialFESpace(V0, g)

@AlexandreMagueresse what about returning an error (or warning) if a non transient FE space is given to TransientLinearFEOperator ?

[urbainvaes]

@Antoinemarteau many thanks!

In fact, tutorial 17 was our starting point. In our initial, code, we used the approach you suggested, but it does not make any difference; with your modification, the example code above still does not work.

Also, in this case I think it would be reasonable to allow use of a TriaFESpace, since the space is not time-dependent. And indeed, it does work as intended with a Dirichlet boundary condition.

[Antoinemarteau]

@urbainvaes After some testing, I agree with all this.

The issue is... a magnificent typo in Tutorial 17: the stiffness matrix goes first and the mass second in the op constructor (a, m) :')

I opened a PR.

[urbainvaes]

With the change a <-> m, the code works, thank you!