Solving a diffusion PDE on a GAMer mesh using FEniCS

This tutorial will walk you through the process of setting up and solving a diffusion partial differential equation (PDE) using a mesh imported from GAMer into FEniCS

import pygamer
import math

Generating the mesh using GAMer

# generate the surface mesh for a cube
smesh = pygamer.surfacemesh.cube(5)

# Mark sides of the face to use later for boundary conditions
for fid in smesh.faceIDs:
    indices = fid.indices()
    z_coords = [smesh.getVertex([idx]).data().position[2] for idx in list(indices)]
    if all([z==1 for z in z_coords]):
        fid.data().marker = 2
    elif all([z==-1 for z in z_coords]):
        fid.data().marker = 3
    else:
        fid.data().marker = 4

# construct the volumetric mesh
vmesh = pygamer.makeTetMesh([smesh], "q1.1/15O10/2ACVY")

# write volumetric mesh into dolfin readable format
pygamer.writeDolfin('outputTetMesh.xml', vmesh)
pygamer.writeVTK('outputTetMesh.vtk', vmesh)

Number of vertices: 1538
Number of Faces: 3072
Number of Regions: 1
Number of Holes: 0
Region point: Tensor({0}:-0.875; {1}:0.333333; {2}:-0.333333)

Setup the model problem in FEniCS

Model problem: strong formulation

We will model diffusion in a \(2 \times 2 \times 2\) cube with an exponentially decaying source term centered at the origin, an influx boundary condition at the \(z=1\) surface, an outflux boundary condition at the \(z=-1\) surface, and a \(u=0\) Dirichlet boundary condition on the four other sides of the cube.

Find \(u(\mathbf{x},t)\) such that

\begin{align} \frac{\partial u}{\partial t} &= D\nabla^2 u + f ~~\text{in}~~ \Omega \\ D(n \cdot \nabla u) &= j_\text{in} ~~\text{on}~~ \partial\Omega_\text{in} \\ D(n \cdot \nabla u) &= j_\text{out} ~~\text{on}~~ \partial\Omega_\text{out} \\ u(\mathbf{x},t=0) &= u_0 \\ u(\mathbf{x},t) &= 0 ~~\text{on}~~ \partial\Omega_\text{dirichlet} \\ \partial\Omega &= \partial\Omega_\text{in} \cup \partial\Omega_\text{out} \cup \partial\Omega_\text{dirichlet} \end{align}

where \(f\) is the source term, \(f(\mathbf{x}) = 32 - e^{||\mathbf{x}||}\), \(D\) is the diffusion coefficient, \(j_\text{in}\) and \(j_\text{out}\) are the influx/outflux rates for the Neumann boundary conditions, and \(u_0\) is the initial condition.

import dolfin as d

# Model parameters
D    = 1 # diffusion coefficient
u0   = d.Constant(0.001) # initial condition
dt   = 0.1 # time step size
T    = 1.0 # final time
jin  = 10.0 # influx rate
jout = -50.0 # outflux rate
dist = d.Expression("sqrt(pow(x[0],2) + pow(x[1],2) + pow(x[2],2))", degree=1)
f    = d.Expression("32 - exp(dist)", dist=dist, degree=1) # source term
t    = 0.0 # current time

numsteps = round(T/dt)

Model problem: weak/variational formulation

The finite element method uses a variational formulation of the PDE which allows us to search for solutions from a larger function space. Equations input into FEniCS are done so through their variational formulation; we will now show how to convert our problem into a valid FEniCS input.

First we multiply the governing PDE by a test function, \(v\), coming from a function space \(V\), integrate over the domain and apply the divergence theorem to obtain the following variational formulation:

\begin{equation} \int_\Omega \frac{\partial u}{\partial t}v \,dx = -\int_\Omega D\nabla u \cdot \nabla v \,dx + \int_{\partial\Omega} D(n\cdot\nabla u)v \,ds + \int_\Omega fv \,dx ~~\text{in}~~ \Omega \end{equation}

We discretize the time-derivative with a backward Euler scheme due to its simplicity and unconditional stability properties

\begin{equation} \frac{\partial u}{\partial t} \approx \frac{u - u_n}{\Delta t} \end{equation}

where \(u_n\) represents the (known) solution computed at the previous timestep. We further simplify the variational formulation by inserting the Neumann boundary conditions and using the shorthand notation where \(\langle \cdot , \cdot \rangle\) represents the inner-product over \(\Omega\) and \(\langle \cdot , \cdot \rangle_{\partial\Omega_\text{in}}\) represents the inner-product on the boundary \(\partial\Omega_\text{in}\). Terms are separated onto the left or right hand sides by their dependence on \(u\):

\begin{equation} \langle u,v \rangle + D\Delta t \langle \nabla u,\nabla v \rangle = \Delta t\langle f,v \rangle + \Delta t \langle j_\text{in},v \rangle_{\partial\Omega_\text{in}} + \Delta t \langle j_\text{out},v \rangle_{\partial\Omega_\text{out}} + \langle u_n,v \rangle \end{equation}

Notice that since \(u_n\) is a known value it is not dependent on the unknown, \(u\), and therefore is placed on the right hand side. Also, since \(D(\mathbf{n}\cdot \nabla u)=0\) on the no-flux boundary those terms drop out of the variational formulation.

In the abstract form this is written as,

\begin{align} a(u,v) &= \langle u,v \rangle + D\Delta t \langle \nabla u,\nabla v \rangle \\ L(v) &= \Delta t\langle f,v \rangle + \Delta t \langle j_\text{in},v \rangle_{\partial\Omega_\text{in}} + \Delta t \langle j_\text{out},v \rangle_{\partial\Omega_\text{out}} + \langle u_n,v \rangle \end{align}

where \(a(u,v)\) is a bilinear form and \(L(v)\) is a linear functional. Note that the Dirichlet boundary condition does not appear in the variational form and therefore must be enforced separately.

# Import the mesh and construct a linear Lagrange function space over the mesh
dolfin_mesh = d.Mesh('outputTetMesh.xml')
V = d.FunctionSpace(dolfin_mesh,'P',1)

# Define trial and test functions
u = d.TrialFunction(V)
v = d.TestFunction(V)
# Known solution corresponding to previous timestep. Initialize with initial condition
un = d.interpolate(u0,V)

# Earlier we marked faces on the influx boundary with 2 and faces on the outflux boundaries as 3

# define integration domains
face_dimension = 2 # we marked faces which have a topological dimension of 2
meshBoundary = d.MeshFunction("size_t",dolfin_mesh,face_dimension,dolfin_mesh.domains())
dx = d.Measure('dx',domain=dolfin_mesh)
ds = d.Measure('ds',domain=dolfin_mesh,subdomain_data=meshBoundary)

# define bilinear and linear forms
a = u*v*dx + D*dt*d.inner(d.grad(u),d.grad(v))*dx
L = dt*f*v*dx + dt*jin*v*ds(2) + dt*jout*v*ds(3) + un*v*dx

Solve and store solution

Solutions are saved in the \(\texttt{.vtk}\) format and can be visualized using software such as Paraview

filename = 'solution/dolfinOut.pvd'
vtkfile = d.File(filename)

# store the initial condition
u = d.Function(V)
un.rename("u", "solution")
u.rename("u", "solution")
vtkfile << (un,t)

bc = d.DirichletBC(V,d.Constant(0),meshBoundary,4)

# Main loop
for idx in range(numsteps):
    # step forward in time
    t += dt

    d.solve(a==L,u,bcs=bc)

    # save solution to file
    vtkfile << (u,t)

    # assign just computed solution to un
    un.assign(u)

    print('Time step %d complete ...' % int(idx+1))

Time step 1 complete ...
Time step 2 complete ...
Time step 3 complete ...
Time step 4 complete ...
Time step 5 complete ...
Time step 6 complete ...
Time step 7 complete ...
Time step 8 complete ...
Time step 9 complete ...
Time step 10 complete ...