PDE Project Course

4. DOLFIN, Systems of PDEs

Johan Jansson

johanjan@math.chalmers.se

Department of Computational Mathematics

Lecture plan

• Goal of DOLFIN

• Solving a PDE with DOLFIN • Systems of PDEs

Goal of DOLFIN

PDE Project Course 04/05 – p. 3

Goal of DOLFIN

• Automatic solution of PDE (variational formulation)

• Why?

• Efficiency

• Execution time

• Development time

• Safety - reduce chance of human error (bugs)

• Solving PDEs:

• Modeling - problem-specific

• Method - standard methods (FEM, FD)

• General solver - still open problem

• How to construct an efficient general solver for PDEs? (instead of making a special solver for every problem)

• How to construct an efficient general solver for PDEs?

• DOLFIN, Puffin try to answer question with:

(A_{h})_{ij} = a(φ_{j}, φ̂_{i}).

• DOLFIN: algebra of basis functions to compute general form

• Likely answer: form compilation (precompute basis function expression and integrals)

Solving a PDE with DOLFIN

PDE Project Course 04/05 – p. 7

Solving Poisson’s equation

int main()
{
Mesh mesh("mesh.xml.gz");
Problem poisson("poisson", mesh);
poisson.set("source", f);
poisson.set("boundary condition", mybc);
poisson.solve();
return 0;
}

Implementing a solver

void PoissonSolver::solve()
{
Galerkin      fem;
Matrix        A;
Vector        x, b;
Function      u(mesh, x);
Function      f(mesh, "source");
Poisson       poisson(f);
KrylovSolver solver;
File          file("poisson.m");
fem.assemble(poisson, mesh, A, b);
solver.solve(A, x, b);
u.rename("u", "temperature");
file << u;
}

Automatic assembling

class Poisson : public PDE {
...
real lhs(const ShapeFunction& u, const ShapeFunction& v)
{
  return (grad(u),grad(v)) * dK;
}
real rhs(const ShapeFunction& v)
{
  return f*v * dK;
}
...
};

Automatic assembling

class ConvDiff : public PDE {
...
real lhs(const ShapeFunction& u, const ShapeFunction& v)
{
  return (u*v + k*((b,grad(u))*v + a*(grad(u),grad(v))))*dK;
}
real rhs(const ShapeFunction& v)
{
  return (up*v + k*f*v) * dK;
}
...
};

Abstract assembly algorithm

(A_{h})_{ij} = a(φ_{j}, φ̂_{i}) = ∫_{Ω} A(φ_{j})φ̂_{i} dx
= ∑_{K∈T} ∫_{K} A(φ_{j})φ̂_{i} dx = ∑_{K∈T} a(φ_{j}, φ̂_{i})_{K}.

Iterate over all elements K and for each element K compute the contributions to all (A_{h})_{ij}, for which φ_{j} and φ̂_{i} are supported within K.

Assembly in DOLFIN

void FEM::assembleInterior(PDE& pde, Mesh& mesh, Matrix& A,
FiniteElement::Vector& element, Map& map,
Quadrature& interior_quadrature,
Quadrature& boundary_quadrature)
{
// Iterate over all cells in the mesh
for (CellIterator cell(mesh); !cell.end(); ++cell)
{
// Update map
map.update(*cell);
// Update element
for (unsigned int i = 0; i < element.size(); ++i)
  element(i)->update(map);
// Update equation
pde.updateLHS(element, map, interior_quadrature, boundary_quadrature);
// Iterate over test and trial functions
for (FiniteElement::Vector::TestFunctionIterator v(element); !v.end(); ++v)
  for (FiniteElement::Vector::TrialFunctionIterator u(element); !u.end(); ++u)
     A(v.dof(*cell), u.dof(*cell)) += pde.lhs(*u, *v);
}
}

Demonstration

•...

Systems of PDEs

PDE Project Course 04/05 – p. 15

Systems of PDEs

A, u and f have n components:

A_{1}(u_{1}) = f_{1}
A_{2}(u_{2}) = f_{2}
...
A_{n}(u_{n}) = f_{n}

Variational formulation

• V^{n} space of functions with n components.

• V_{h}^{n} finite dimensional space of functions with n components.

• Assume n = 2 for compactness of notation.

• u = (u_{1}, u_{2}), v = (v_{1}, v_{2}), f = (f_{1}, f_{2})

A(u) = f ⇒
∫_{Ω} A(u)v = ∫_{Ω} fv ⇒

PDE Project Course 04/05 – p. 17

Variational formulation

\begin{aligned}∫_{Ω}(A_{1}(u_{1}), A_{2}(u_{2}))(v_{1}, v_{2}) & = & ∫_{Ω}(f_{1}, f_{2})(v_{1}, v_{2}) ⇒ \\ ∫_{Ω} A_{1}(u_{1})v_{1} + A_{2}(u_{2})v_{2} & = & ∫_{Ω} f_{1}v_{1} + f_{2}v_{2} ⇒ \\ a(u, v) & = & l(v) \end{aligned}

PDE Project Course 04/05 – p. 18

Finite element basis

• Every component of V_{h}^{n} is a scalar V_{h}.

• What is a basis for V_{h}^{n}?

• Standard basis for a vector space (n = 2):

ϕ̂_{1} = (ϕ_{1}, 0)
ϕ̂_{2} = (0, ϕ_{1})
ϕ̂_{3} = (ϕ_{2}, 0)
ϕ̂_{4} = (0, ϕ_{2})
...
ϕ̂_{2N−1} = (ϕ_{N} , 0)
ϕ̂_{2N} = (0, ϕ_{N})

PDE Project Course 04/05 – p. 19

Ansatz

• The ansatz is just as before:

U(x) = ∑_{j=1}^{N} ξ_{j}φ_{j}(x).

• Concretely:

• ξ_{j} is the dof (value of U) for component 1 of node j

• ξ_{j+1} is the dof for component 2 of node j + 1

Example: linear elasticity

The equations of linear elasticity:

u̇ − v = 0 in Ω^{0},
v̇ − ∇ · σ = f in Ω^{0},
σ = Eϵ(u) = E(∇u^{⊤} + ∇u)
v(0, ·) = v^{0}, u(0, ·) = u^{0} in Ω^{0}.

Variational formulation

• ϵ(u) = ∇u^{⊤} + ∇u

• σ = Eϵ(u), E linear transformation (material properties)

• Turns out that:

\begin{aligned}a(u, v) & = & ∫_{Ω} σ(u)ϵ(v) \\ l(v) & = & ∫_{Ω} fv \end{aligned}