Skip to content

Instantly share code, notes, and snippets.

@wence-
Created January 15, 2014 11:21
Show Gist options
  • Select an option

  • Save wence-/8434549 to your computer and use it in GitHub Desktop.

Select an option

Save wence-/8434549 to your computer and use it in GitHub Desktop.

Solving PDEs

Problem setup

Consider a weak variational problem

a(u, v) = L(v) \; \forall v \in V \mathrm{on} \Omega
u = u_0 \; \mathrm{on} \partial\Omega

The bilinear and linear parts of this form are represented as UFL objects, call them a and L respectively. Strongly imposed boundary conditions are represented as DirichletBC objects:

bc = DirichletBC(V, u_0, subdomain)

Where subdomain may either be an integer id describing the region over which the boundary condition should be applied, or, for extruded meshes, the strings "top" or "bottom" indicating that the condition applies at the top and bottom of the extruded column respectively. u_0 may either be a numeric value or an Expression.

Now that we have all the pieces of our variational problem, we can move forward to solving it.

Solving the variational problem

The function used to solve PDEs defined as above is solve. This is a unified interface for solving both linear and non-linear variational problems along with linear systems (where the arguments are already assembled matrices and vectors, rather than UFL forms). We will treat the variational interface first.

Linear variational problems

If the problem is linear, that is a is linear in both the test and trial functions and L is linear in the test function, we can use the linear variational problem interface to solve. To start, we need a Function to hold the value of the solution:

s = Function(V)

We can then solve the problem, placing the solution in s with:

solve(a == L, s)

To apply boundary conditions, one passes a list of DirichletBC objects using the bcs keyword argument. For example, if there are two boundary conditions, in bc1 and bc2, we write:

solve(a == L, s, bcs=[bc1, bc2])

Nonlinear variational problems

For nonlinear problems, the interface is similar. In this case, we solve a problem:

F(u; v) = 0 \; \forall v \in V \mathrm{on} \Omega
u = u_0 \; \mathrm{on} \partial\Omega

where the residual F is linear in the test function v but possibly non-linear in the unknown Function u. To solve this problem we write, if F is the residual form:

solve(F == 0, u)

to apply strong boundary conditions, as before, we provide a list of DirichletBC objects using the bcs keyword:

solve(F == 0, u, bcs=[bc1, bc2])

nonlinear problems in Firedrake are solved using Newton-like methods. That is, we compute successive approximations to the solution using

u_{k+1} = u_{k} - J(u_k)^{-1} F(u_k) \; k = 0, 1, \dots

where u_0 is an initial guess for the solution and J(u_k) = \frac{\partial F(u_k)}{\partial u_k} is the Jacobian of the residual, which should be non-singular at each iteration.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment