INF5620 - Numerical methods for partial differential equations

$$ \newcommand{\uex}{u_{\mbox{\footnotesize e}}} \newcommand{\uexd}[1]{u_{\mbox{\footnotesize e}, #1}} \newcommand{\vex}{v_{\mbox{\footnotesize e}}} \newcommand{\vexd}[1]{v_{\mbox{\footnotesize e}, #1}} \newcommand{\Aex}{A_{\mbox{\footnotesize e}}} \newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\Ddt}[1]{\frac{D #1}{dt}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\Std}[1]{\hbox{Std}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \newcommand{\X}{\boldsymbol{X}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\w}{\boldsymbol{w}} \newcommand{\V}{\boldsymbol{V}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\strain}{\boldsymbol{\varepsilon}} \newcommand{\stressc}{{\sigma}} \newcommand{\strainc}{{\varepsilon}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\ii}{\boldsymbol{i}} \newcommand{\jj}{\boldsymbol{j}} \newcommand{\kk}{\boldsymbol{k}} \newcommand{\ir}{\boldsymbol{i}_r} \newcommand{\ith}{\boldsymbol{i}_{\theta}} \newcommand{\iz}{\boldsymbol{i}_z} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\If}{\mathcal{I}_s} % for FEM \newcommand{\Ifd}{{I_d}} % for FEM \newcommand{\Ifb}{{I_b}} % for FEM \newcommand{\setb}[1]{#1^0} % set begin \newcommand{\sete}[1]{#1^{-1}} % set end \newcommand{\setl}[1]{#1^-} \newcommand{\setr}[1]{#1^+} \newcommand{\seti}[1]{#1^i} \newcommand{\sequencei}[1]{\left\{ {#1}_i \right\}_{i\in\If}} \newcommand{\basphi}{\varphi} \newcommand{\baspsi}{\psi} \newcommand{\refphi}{\tilde\basphi} \newcommand{\psib}{\boldsymbol{\psi}} \newcommand{\sinL}[1]{\sin\left((#1+1)\pi\frac{x}{L}\right)} \newcommand{\xno}[1]{x_{#1}} \newcommand{\Xno}[1]{X_{(#1)}} \newcommand{\yno}[1]{y_{#1}} \newcommand{\Yno}[1]{Y_{(#1)}} \newcommand{\xdno}[1]{\boldsymbol{x}_{#1}} \newcommand{\dX}{\, \mathrm{d}X} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integerp}{\mathbb{N}} \newcommand{\Integer}{\mathbb{Z}} $$

Default choice of the final project in INF5620

Project 1: Nonlinear diffusion equation

Project size: approximately 40 h. This particular project will have significant overlap with topics for the final exam.

The goal of this project is to discuss various numerical aspects of a nonlinear diffusion model:

$$ \varrho u_t = \nabla\cdot (\alpha (u)\nabla u) + f(\pmb{x}, t),$$ with initial condition $ u(\pmb{x},0)=I(x) $ and boundary condition $ \partial u/\partial n=0 $. The coefficiet $ \varrho $ is constant and $ \alpha(u) $ is a known function of $ u $.

a) Introduce some finite difference approximation in time that leads to an implicit scheme (say a Backward Euler, Crank-Nicolson, or 2-step backward scheme). Derive a variational formulation of the initial condition and the spatial problem to be solved at each time step (and at the first time step in case of a 2-step backward scheme).

b) Formulate a Picard iteration method at the PDE level, using the most recently computed $ u $ function in the $ \alpha(u) $ coefficient. Derive general formulas for the entries in the linear system to be solved in each Picard iteration. Use the solution at the previous time step as initial guess for the Picard iteration.

c) Restrict the Picard iteration to a single iteration. That is, simply use a $ u $ value from the previous time step in the $ \alpha(u) $ coefficient. Implement this method with the aid of the FEniCS software (in a dimension-independent way).

d) The first verification of the FEniCS implementation may assmue $ \alpha (u)=1 $, $ f=0 $, $ \Omega = [0,1]\times [0,1] $, P1 elements, and $ I(x,y)= \cos (\pi x) $. The exact solution is then $ u(x,y,t)=e^{-\pi^2 t}\cos (\pi x) $. The error in space is then $ {\cal O}(\Delta x^2) + {\cal O}(\Delta y^2) $, while the error in time is $ {\cal O}(\Delta t^p) $, with $ p=1 $ for the Backward Euler scheme and $ p=2 $ for the Crank-Nicolson or the 2-step backward schemes. Set $ h=\Delta t^p = \Delta x^2 $ and perform a convergence rate study as $ h $ is decreased. A measure of the error $ E $ can be computed by

e = u_e.vector().array() - u.vector().array()
E = numpy.sqrt(numpy.sum(e**2)/u.vector().array().size)

for some fixed point of time, say $ t=0.05 $. Here, u_e is a projection of the exact solution onto the function space used for u.

Show that $ E/h $ remains approximately constant as the mesh in space and time is simultaneously refined.

e) To get an indication whether the implementation of the nonlinear diffusion PDE is correct or not, we can use the method of manufactured solutions. Say we restrict the problem to one space dimension, $ \Omega=[0,1] $, and choose $$ u(x,t) = t\int_0^x q(1-q)dq = tx^2\left(\frac{1}{2} - \frac{x}{3}\right)$$ and $ \alpha(u) = 1+u^2 $. The following sympy session computes an $ f(x,t) $ such that the above $ u $ is a solution of the PDE problem:

>>> from sympy import *
>>> x, t, rho, dt = symbols('x t rho dt')
>>>
>>> def a(u):
...     return 1 + u**2
...
>>> def u_simple(x, t):
...     return x**2*(Rational(1,2) - x/3)*t
...
>>> # Show that u_simple satisfies the BCs
>>> for x_point in 0, 1:
...     print 'u_x(%s,t):' % x_point,
...     print diff(u_simple(x, t), x).subs(x, x_point).simplify()
...
u_x(0,t): 0
u_x(1,t): 0
>>> print 'Initial condition:', u_simple(x, 0)
Initial condition: 0
>>>
>>> # MMS: full nonlinear problem
>>> u = u_simple(x, t)
>>> f = rho*diff(u, t) - diff(a(u)*diff(u, x), x)
>>> print f.simplify()
-rho*x**3/3 + rho*x**2/2 + 8*t**3*x**7/9 - 28*t**3*x**6/9 +
7*t**3*x**5/2 - 5*t**3*x**4/4 + 2*t*x - t

Compare the FEniCS solution and the $ u $ given above as a function of $ x $ for several values of $ t $.

f) List the difference sources of numerical errors in the FEniCS program.

g) (Can be postponed.) To verify the nonlinear PDE implementation in FEniCS, we can eliminate the error due to a single Picard iteration by computing a manufactured solution corresponding to $ \alpha (u_1) $, where $ u_1 $ is the finite element function representation of the solution at the previous time step.

>>> u_1 = u_simple(x, t-dt)
>>> f = rho*diff(u, t) - diff(a(u_1)*diff(u, x), x)
>>> print f.simplify()
rho*x**2*(-2*x + 3)/6 -
(-12*t*x + 3*t*(-2*x + 3))*(x**4*(-dt + t)**2*(-2*x + 3)**2 + 36)/324
- (-6*t*x**2 + 6*t*x*(-2*x + 3))*(36*x**4*(-dt + t)**2*(2*x - 3)
+ 36*x**3*(-dt + t)**2*(-2*x + 3)**2)/5832

Perform a convergence rate test by decreasing $ h=\Delta t^p = \Delta x^2 $.

h) Simulate the nonlinear diffusion of Gaussian function. Due to symmetry, it suffices to simulate one quarter of the domain:

$$ I(x,y) = \exp{\left(-\frac{1}{2\sigma^2}\left( x^2 + y^2 \right)\right)}, \quad (x,y)\in\Omega = [0,1]\times [0,1], $$ where $ \sigma $ measures the width of the Gaussian function. Choose $ \alpha(u) = 1+\beta u^2 $, where $ \beta $ is a constant one can play around with.

The rest of the exercise is theoretical and not directly related to the FEniCS implementation. In 2013 several of these questions coincide with the exam and it is therefore not necessary to address them in this compulsory project.

i) Formulate a group finite element method for the $ \alpha (u) $ coefficient and derive formulas for eqution number $ i $ in the nonlinear algebraic system to be solved.

j) Derive expressions for the entries in the Jacobian matrix and the right-hand side to be used in a Newton method.

k) Reduce the problem to one space dimension. Derive formulas for the entries in the Jacobian matrix and the right-hand side in Newton's method using P1 elements and the Trapezoidal rule for integration.

l) Compare the expressions in the previous subexercise with the expressions arising from a finite difference method applied to PDE problem (using the same time discretization).

Kontakt

Institutt for informatikk

Faglærere

Hans Petter Langtangen