Boundary Conditions with Heat Equation

[carlos-adir]

Hello! Recently I've been working on the heat equation (which is linear and easy) to learn how to use Firedrake. I'm not an expert on finite elements, so I appreciate your patience. If you like the text, you may use it as an example to put in the documentation, as you like.

I was inspired by the code Mixed formulation for the Poisson equation to make both codes below (I'm using colab.research.google.com):

• HeatEquation1Field
• HeatEquation2Fields

A more detailed description is shown within each file. I will put only essential information. If you have any question, please ask me.

The simulation is mainly to find the temperature field T(scalar). To do so, I have two boundaries called GamD(Dirichlet Boundary) and GamN(Neumann Boundary). In GamD I know T0 while in GamN I know the flux q0.

---

1. In the first example: HeatEquation1Field

• I rewrite the equation to get the laplacian of T
• I use only one test function v
• The variational term L has the integral over GamN
• I apply DirichletBC: T0 on GamD
• and it works fine!

I did not like this alternative because a boundary condition is in the variational form and I would like to apply it using DirichletBC function (It seems faster to assemble the matrix). So, I've tried the second example.

---

2. In the second example: HeatEquation2Fields

• I kept T and q in separated equations as independent variables
• I use two test functions v and w
• The variational term L doesn't have any integral over the boundary
• I apply the DirichletBC: T0 on GamD and q0 on GamN
• It doesn't work!

I thought that I could have a system of equations such that every boundary condition could be a Dirichlet condition of its respective variable.

In the bottom part of the code, there is a section Old equations where I found an expression with the integral over both boundaries, but it doesn't work either.

---

The main point is about the boundary conditions. Why did it not work in the second case?

PS1: In the Mixed formulation for the Poisson equation, there is a u0 where it was supposed to be non-zero (as example), but it's not even mentioned in the code.
PS2: I could be content with the first code, but further I will work with more complex equations in which I cannot replace a term in another like the static model in continuum mechanics:

div(tensor(sigma)) + vec(f) = vec(0)
tensor(sigma) = lambda * trace( tensor(epsilon) ) * tensor(Identity) + 2 * mu * tensor(epsilon)
tensor(epsilon) = (1/2) * grad_sym( vec(u))

I appreciate your help!

[ksagiyam]

If I copy/paste your code and run locally, it runs fine:

from firedrake import *
  

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, family = "CG", degree = 3)  # Scalar
W = VectorFunctionSpace(mesh, family = "CG", degree = 2, dim=2)  # Vetorial
MixedField = V * W
T, q = TrialFunctions(MixedField)
v, w = TestFunctions(MixedField)
x, y = SpatialCoordinate(mesh)
f = Function(V)
T0 = Function(V)  # Identicamente nulo
q0 = Function(V)  # Identicamente nulo
a = (v * div(q) + inner(w, q) + inner(w, grad(T)))*dx
L = f * v * dx
bc0 = DirichletBC(MixedField.sub(0), g=0, sub_domain=1)  # Left Line
bc1 = DirichletBC(MixedField.sub(0), g=1, sub_domain=2)  # Right Line
bc2 = DirichletBC(MixedField.sub(1), g=as_vector([0, 0]), sub_domain=3)  # Bottom line
bc3 = DirichletBC(MixedField.sub(1), g=as_vector([0, 0]), sub_domain=4)  # Top Line
bcs = [bc0, bc1, bc2, bc3]
u = Function(MixedField)
solve(a == L, u, bcs=bcs)

Maybe you forgot to run some cells after some changes?

[carlos-adir]

Does it give good results? For me, it solves without any error, but the temperature field T is like that:

While it should be like this one:

[ScottMacLachlan]

It isn't clear to me that the problem that you're solving is well-posed as written. You're using H^1 vector spaces for both T and q, but your mixed formulation is almost certainly not well-posed in the product space of scalar H^1 and vector H^1. The standard solution for this is the one in the mixed Poisson example that you point to, where we solve for T in an L^2 space (DG) and q in an H(div) space (BDM, in that example, or RT). The BC on T enters the weak formulation when you integrate by parts on inner(w, grad(T))*dx which gives you -inner(div(w),T)*dx plus a boundary term. That boundary term should vanish over GamN but not over GamD, where it contributes to the linear form L.

In the example, u0 = 0 is specified, so this term vanishes altogether. If you include it in your form in place of the DirichletBC terms that you now use, you should get the correct solution.

[carlos-adir]

It worked! In fact ScottMacLachlan, you were right about the space functions: W bases' functions must be vectorial functions, not the weights (as I was doing). I corrected it and I put the same functions (V as DG and W as BDM) as in the Mixed Poisson example.

But I still have a question about the weak formulation: I can rewrite the weak formulation in some ways, as shown in the file HeatEquation2Fields. So, I can get 4 assemblies of (a, L) like (a1, L1), (a2, L2), (a3, L3), (a4, L4):

• L1 has no integral over the boundary
• L2 has a integral over GamN
• L3 has a integral over GamD
• L4 have a integral over GamD and another one over GamN
  I tested with all of them but only the third converged (when the integral over GamN vanish), but I still don't know why, cause I always put the boundary conditions on GamD and GamN using the function DirichletBC.
  In other problems, should I make all the possibilities of weak formulations and test them one by one? Or there's a way to know which weak formulation is the "best"?

About the Mixed Poisson Example, if u0 = 1, the modifications would be:

u0 = Function(DG) + 1  # new
n = FacetNormal(mesh)  # new
L = - f*v*dx 
L += u0 * inner(tau, n) * ds  # new
bc0 = DirichletBC(W.sub(0), as_vector([0.0, -sin(5*x)]), 3)
bc1 = DirichletBC(W.sub(0), as_vector([0.0, sin(5*x)]), 4)
bc2 = DirichletBC(W.sub(1), u0, 1)  # new
bc3 = DirichletBC(W.sub(1), u0, 2)  # new
# solve(a == L, w, bcs=[bc0, bc1])
solve(a == L, w, bcs=[bc0, bc1, bc2, bc3])  # new 

[ScottMacLachlan]

The key here is that your formulations 1, 2, and 4 aren't allowed, because they involve gradients of T and/or v, and these come from the discontinuous space, so their gradients aren't in L2 over the mesh. Firedrake assembles the forms element-wise, but doesn't do anything to handle the discontinuities over the element boundaries, so you aren't actually consistently discretizing your PDE. Form 3, in contrast, is fine, since the only derivatives it contains are div(q) and div(w), and q and w are discretized in H(div), so their divergences are in L2.

Note also that Firedrake does not strongly impose boundary conditions on DG fields, even if you tell it to do so. So, your bc0 and bc1 in the colab notebook (bc2 and bc3 in the snippet above, with the subspaces in the opposite order) actually don't do anything. You should see the same results with equation = 3 and those two lines omitted. Canonically, DG fields have no DoFs on the boundary of the domain (as all DG DoFs are interior to an element), so there are no DoFs on which to impose strong (Dirichlet) BCs.

[ScottMacLachlan]

To fully answer this question:
In other problems, should I make all the possibilities of weak formulations and test them one by one? Or there's a way to know which weak formulation is the "best"?
is difficult, but the essence is that you need to choose your discretization space(s) in a way that is compatible with one of the possible forms of the equations, ensuring that all terms that appear in a bilinear form are properly integrable over these spaces. How, precisely, to do this depends greatly on the problem at hand. For many problems, there are known "correct" approaches, documented in the research literature and/or textbooks or monographs, but there are also problems for which these are not yet known.