Help with the biharmonic equation

[paolopiers]

Hello Everyone,
I am new to this forum and to Firedrake. I am trying to implement the non homogeneous biharmonic equation with Dirichlet and Neumann boundary conditions (i.e., the solution is sought in the Sobolev space H^2_0).

I manage to obtain a result but the boundary conditions are not correctly implemented. The boundary conditions are implemented via the penalty method.

Could you please take a look at the enclosed code and tell me what is wrong? Sorry, I was not able to format the code correctly using the quotes.

Thanks in advance for your help.

from ast import Expression
from pyclbr import Function
import math
import sympy as sym
import matplotlib as plt
from matplotlib import pyplot
import numpy as np
from firedrake import *
mesh2 = RectangleMesh(100, 100, math.pi, 2.0)
f = Constant(0.1)
W = FunctionSpace(mesh2,"Argyris",5)
#W = FunctionSpace(mesh2,"Morley",2)
nu = FacetNormal(mesh2)
u_tr = TrialFunction(W)
u_test = TestFunction(W)
tgv=10000
a=inner(div(grad(u_tr)),div(grad(u_test)))*dx+tgv*inner(grad(u_tr),nu)*u_test*ds+tgv*u_tr*u_test*ds
L=f*u_test*dx
u=Function(W)
#bc = boundary_L(W, Constant(0), "on_boundary")
#bc = DirichletBC(W, Constant(0.0),"on_boundary")
solve(a == L, u)
File("helmholtz.pvd").write(u)

[dham]

I cut this back to remove all the unnecessary imports:

from firedrake import *
mesh2 = RectangleMesh(100, 100, pi, 2.0)
f = Constant(0.1)
W = FunctionSpace(mesh2,"Argyris",5)
#W = FunctionSpace(mesh2,"Morley",2)
nu = FacetNormal(mesh2)
u_tr = TrialFunction(W)
u_test = TestFunction(W)
tgv=10000
a=inner(div(grad(u_tr)),div(grad(u_test)))*dx+tgv*inner(grad(u_tr),nu)*u_test*ds+tgv*u_tr*u_test*ds
L=f*u_test*dx
u=Function(W)
#bc = boundary_L(W, Constant(0), "on_boundary")
#bc = DirichletBC(W, Constant(0.0),"on_boundary")
solve(a == L, u)
File("helmholtz.pvd").write(u)

It runs for me and the answer doesn't look obviously wrong. There is some error in achieving the boundary conditions but given that you have a uniform forcing function which will be pushing against the boundary conditions, I'm not sure this is obviously incorrect. This is the danger with trying to assess the correctness of a solution without an anaytic solution.

I suggest you use the method of manufactured solutions (MMS) to create a forcing function for a solution which satisfies the equation and BCs. I might choose $(\cos(\pi x[0])+1) (\cos(\pi x[1])+1)$ as my solution and drop that into the equation to work out the correct forcing. For a little more on MMS see: https://finite-element.github.io/6_finite_element_problems.html#the-method-of-manufactured-solutions

[pefarrell]

@rckirby and @wence- have archived code for solving the biharmonic equation in a bunch of ways here: https://zenodo.org/record/3381901/files/code-data.tar.gz?download=1

(relevant links: https://arxiv.org/abs/1808.05513 , https://zenodo.org/record/3381901 )

[paolopiers]

Thank you very much, @dham and @pefarrell , for your help.

@dham : The theory of the biharmonic equation (cf., e.g., Section 6.8 of P. G. Ciarlet - Linear and Nonlinear Functional Analysis with Applications (2013)) states that under the sole assumption that the force is of class L^2, there exists a unique solution with Dirichlet and Neumann boundary conditions on the contour. Therefore, it should work.

@pefarrell Thank you for the links. I have the folder with the code discussed in Kirby's paper. My idea is to extend the biharmonic equation code to a more complex one describing deformations in elasticity. Unfortunately, the code does not function well. When I tried to implement the boundary conditions of Dirichlet type in the standard way, I am prompted a message stating that this feature is not implemented yet and that I have to resort to the penalization machinery, that we saw is not exactly as precise as it ought to be.

Any ideas to improve my code?

[wence-]

The code @pefarrell references in Kirby & Mitchell (2018) applies Dirichlet conditions using Nitsche's method which your simple example does not do, so you may want to start there. The problem with not being able to apply Dirichlet conditions directly to the function space (given the particular implementation of H^2 elements in Firedrake) is discussed in the paper: effectively the set of nodes one must set to values is not unique.

When you say the code from our paper "does not function well", what does that mean?

[paolopiers]

Thank you, I reviewed the code and now it works! I can see the Neumann and Dirichlet boundary conditions very clearly in my new code, that I enclose below

from firedrake import *
mesh2 = UnitSquareMesh(50,50)
#mesh2 = RectangleMesh(100, 100, math.pi, 2.0)
h = CellSize(mesh2)
f = Constant(12.0)
W = FunctionSpace(mesh2,"Argyris",5)
nu = FacetNormal(mesh2)
u_tr = TrialFunction(W)
u_test = TestFunction(W)
a=inner(div(grad(u_tr)),div(grad(u_test)))*dx+h*inner(grad(u_tr),nu)*inner(grad(u_test),nu)*ds+pow(h,3)*u_tr*u_test*ds+inner(grad(div(grad(u_tr))),nu)*u_test*ds-div(grad(u_tr))*inner(grad(u_test),nu)*ds+inner(grad(div(grad(u_test))),nu)*u_tr*ds-div(grad(u_test))*inner(grad(u_tr),nu)*ds
L=f*u_test*dx
u=Function(W)
solve(a == L, u)
File("helmholtz.pvd").write(u)

Thank you very much for the prompt help to everyone!

[paolopiers]

One more question (I don't know if I had better ask in a different thread): How do you implement boundary conditions on a portion of the boundary in Firedrake. For instance, only on two edges of the unit square?
Thanks in advance!

[dham]

The Python help on UnitSquareMesh will tell you what the boundary IDs of the different sides of the square are. You supply these as arguments to the ds. E.g. to apply an integral only over boundary segment 2, use ds(2).