Commit 59c595fd authored by Lloyd Connellan's avatar Lloyd Connellan
Browse files

add cells to nvtest and fix errors

parent e87ce887
......@@ -9,6 +9,7 @@ USER root
# Install packages
RUN apt-get update && apt-get dist-upgrade --yes --no-install-recommends
RUN apt-get install --yes --no-install-recommends \
ca-certificates wget tar \
libc6-dev gcc g++ binutils make cmake pkg-config patch \
......@@ -54,38 +55,32 @@ RUN /root/install-dune.sh
# Set up Jupyter config
ADD jupyter_config.py /etc/
RUN chmod 644 /etc/jupyter_config.py
ADD jupyter_config.py /etc/
RUN chmod 644 /etc/jupyter_config.py
# Set up dune as user
RUN adduser --disabled-password --home /dune --gecos "" --uid 50000 dune
ENV LD_LIBRARY_PATH "/usr/lib:/usr/local/lib"
ENV DUNE_CMAKE_FLAGS="CMAKE_BUILD_TYPE=Release BUILD_SHARED_LIBS=TRUE -DPETSC_DIR=/usr/local/petsc-32"
# Copy dune-fempy notebooks and change ownership
COPY dune-fempy.ipynb /dune/
COPY crystal.ipynb /dune/
COPY mcf.ipynb /dune/
COPY laplace-adaptive.ipynb /dune/
RUN chown dune:dune /dune/dune-fempy.ipynb
RUN chown dune:dune /dune/crystal.ipynb
RUN chown dune:dune /dune/mcf.ipynb
RUN chown dune:dune /dune/laplace-adaptive.ipynb
COPY dune-fempy.ipynb crystal.ipynb mcf.ipynb laplace-adaptive.ipynb /dune/
# Copy dune-femnv files and change ownership
COPY nvtest.ipynb /dune/
COPY problems.py /dune/
COPY schemes.py /dune/
COPY results.py /dune/
COPY norms.py /dune/
RUN mkdir /dune/results
RUN mkdir /dune/pickle
RUN mkdir /dune/plots
RUN chown dune:dune /dune/nvtest.ipynb
#WORKDIR /dune
COPY nvtest.ipynb problems.py schemes.py results.py norms.py matrix.py /dune/
RUN mkdir /dune/results /dune/pickle /dune/plots /dune/pickle1 /dune/pickle3
RUN chown dune:dune /dune/results /dune/pickle /dune/plots /dune/pickle1 /dune/pickle3
WORKDIR /dune
# due to the 'MPI Finalaize' error we have to let this run without checking
# for errors - make sure to go through the output of the docker build command to
# check that the notebook ran until the end
RUN chown dune:dune /dune/dune-fempy.ipynb
USER dune
RUN jupyter nbconvert --config=/etc/jupyter_config.py --to notebook\
--ExecutePreprocessor.timeout=-1 --execute dune-fempy.ipynb --output dune-fempy.ipynb \
2>&1 > /dune/build.log || echo "Execution of dune-fempy.ipynb!"
ENV LD_LIBRARY_PATH "/usr/lib:/usr/local/lib"
ENV DUNE_CMAKE_FLAGS="CMAKE_BUILD_TYPE=Release BUILD_SHARED_LIBS=TRUE"
EXPOSE 8888
VOLUME ["/dune"]
......
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Python code from the Dune-FemPy paper\n",
"======================================\n",
"This Jupyter notebook contains the code examples shown in the\n",
"Dune-Python paper.\n",
"The sectioning reflects that of the paper and, while some introductory\n",
"comments are given, we refer to it for the details.\n",
"\n",
"__Note__: Some of the code cells depend on previous ones and cannot be\n",
"executed in arbitrary order."
]
},
{
"cell_type": "code",
"execution_count": null,
......@@ -32,60 +47,169 @@
},
"outputs": [],
"source": [
"from dune.fem.plotting import plotPointData as plot\n",
"from dune.grid import structuredGrid\n",
"grid = structuredGrid([0, 0], [1, 1], [4, 4])\n",
"grid.plot()\n",
"grid.hierarchicalGrid.globalRefine(1)\n",
"grid.plot()\n",
"grid.hierarchicalGrid.globalRefine(-1) # return grid to original state\n",
"\n",
"grid.hierarchicalGrid.globalRefine(-1) # revert grid refinement"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Setting up a discrete function space and some grid function"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import dune.create as create\n",
"space = create.space('lagrange', grid, dimrange=1, order=2)\n",
"\n",
"from ufl import SpatialCoordinate\n",
"x = SpatialCoordinate(space)\n",
"\n",
"initial = 1/2*(x[0]**2 + x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1\n",
"u_h = space.interpolate(initial, name='u_h')\n",
"\n",
"u_h_n = u_h.copy(name=\"previous\")\n",
"\n",
"initial = 1/2*(x[0]**2 + x[1]**2) - 1/3*(x[0]**3 - x[1]**3) + 1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can easily easily integrate grid function"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from dune.fem.function import integrate\n",
"from ufl import sqrt\n",
"norm = sqrt( integrate(grid, initial, 5)[0] )\n",
"print(norm)\n",
"\n",
"mass = integrate(grid, initial, 5)[0]\n",
"print(mass)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and plot them using matplotlib or write a vtk file for postprocessing"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from dune.fem.plotting import plotPointData as plot\n",
"plot(initial, grid=grid)\n",
"grid.writeVTK('initial', pointdata={'initial': initial})\n",
"grid.writeVTK('initial', pointdata={'initial': initial})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So far we used grid functions defined globally. An important subclass of\n",
"grid functions are discrete functions over a given discrete function space.\n",
"The easiest way to construct such functions is to use the interpolate\n",
"method on the discrete function space:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"u_h = space.interpolate(initial, name='u_h')\n",
"\n",
"from ufl import TestFunction, TrialFunction, triangle\n",
"u_h_n = u_h.copy(name=\"previous\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can set up our PDE model\n",
"\\begin{equation}\n",
" \\int_{\\Omega} u^{n+1} v + \\Delta t K(\\nabla u) \\nabla u^{n+1} \\cdot \\nabla v\\mathop{dx}\n",
" = \\int_{\\Omega} u^n v + \\Delta t f v\\ \\mathop{dx}\n",
"\\end{equation}\n",
"where the diffusion tensor is given by\n",
"\\begin{equation}\n",
" K(\\nabla u) = \\frac{2}{1+\\sqrt{1+4\\nabla u}}\n",
"\\end{equation}\n",
"and we force the system by defining $f$ so that\n",
"\\begin{equation}\n",
" u(x,t) = e^{-2t}\\left(\\frac{1}{2}(x_0^2 + x_1^2) - \\frac{1}{3}(x_0^3 - x_1^3)\\right) + 1\n",
"\\end{equation}\n",
"becomes the exact solution."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from ufl import TestFunction, TrialFunction\n",
"from dune.ufl import NamedConstant\n",
"u = TrialFunction(space)\n",
"v = TestFunction(space)\n",
"dt = NamedConstant(triangle, \"dt\") # time step\n",
"t = NamedConstant(triangle, \"t\") # current time\n",
"dt = NamedConstant(space, \"dt\") # time step\n",
"t = NamedConstant(space, \"t\") # current time\n",
"\n",
"from ufl import dx, grad, inner\n",
"from ufl import dx, grad, inner, sqrt\n",
"abs_du = sqrt(inner(grad(u), grad(u)))\n",
"K = 2/(1 + sqrt(1 + 4*abs_du))\n",
"a = (inner((u - u_h_n)/dt, v) + inner(K*grad(u), grad(v)))*dx\n",
"\n",
"from ufl import exp\n",
"exact = exp(-2*t)*(initial - 1) + 1\n",
"\n",
"from ufl import as_vector, replace\n",
"exact = as_vector([exp(-2*t)*(initial - 1) + 1])\n",
"b = replace(a, {u: exact})\n",
"from ufl import as_vector, exp\n",
"exact = lambda t: as_vector([exp(-2*t)*(initial - 1) + 1])\n",
"\n",
"from ufl import replace\n",
"b = replace(a, {u: exact(t)})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the model described as a ufl form, we can construct a scheme class\n",
"that provides the solve method which we can use to evolve the solution from\n",
"one time step to the next:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"scheme = create.scheme(\"galerkin\", a == b, solver='cg')\n",
"\n",
"scheme.model.dt = 0.05\n",
"\n",
"exact_end = as_vector( [exp(-2)*(initial - 1) + 1] )\n",
"l2error_fn = inner(u_h - exact_end, u_h - exact_end)\n",
"h1error_fn = inner(grad(u_h - exact_end), grad(u_h - exact_end))\n",
"\n",
"def evolve(scheme, u_h, u_h_n):\n",
" time = 0\n",
" endTime = 1.0\n",
......@@ -93,7 +217,31 @@
" scheme.model.t = time\n",
" u_h_n.assign(u_h)\n",
" scheme.solve(target=u_h)\n",
" time += scheme.model.dt\n",
" time += scheme.model.dt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since we have forced the system towards a given solution, we can compute\n",
"the discretization error. First we define ufl expressions for the $L^2$\n",
"and $H^1$ norms and will use those to compute the experimental order of\n",
"convergence of the scheme by computing the time evolution on different grid\n",
"levels."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"exact_end = exact(1)\n",
"l2error_fn = inner(u_h - exact_end, u_h - exact_end)\n",
"h1error_fn = inner(grad(u_h - exact_end), grad(u_h - exact_end))\n",
"\n",
"from math import log\n",
"error = 0\n",
......@@ -123,6 +271,13 @@
"packages and python functionality."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Implementing a simple Newton Krylov solver using the dune-fem linear solvers"
]
},
{
"cell_type": "code",
"execution_count": null,
......@@ -157,11 +312,27 @@
"\n",
"scheme_cls = Scheme(scheme)\n",
"\n",
"grid.hierarchicalGrid.globalRefine(-2) # return grid to original state\n",
"u_h.interpolate(initial) # reset u_h to original state\n",
"grid.hierarchicalGrid.globalRefine(-2) # revert grid refinement\n",
"u_h.interpolate(initial) # reset u_h to initial\n",
"evolve(scheme_cls, u_h, u_h_n)\n",
"plot(u_h)\n",
"\n",
"plot(u_h)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using a non linear solver from the Scipy package"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from scipy.optimize import newton_krylov\n",
"from scipy.sparse.linalg import LinearOperator\n",
"\n",
......@@ -197,7 +368,6 @@
" def __init__(self, scheme):\n",
" self.scheme = scheme\n",
" self.model = scheme.model\n",
"\n",
" def solve(self, target=None):\n",
" sol_coeff = target.as_numpy\n",
" # call the newton krylov solver from scipy\n",
......@@ -208,20 +378,53 @@
"scheme2_cls = Scheme2(scheme)\n",
"u_h.interpolate(initial)\n",
"evolve(scheme2_cls, u_h, u_h_n)\n",
"plot(u_h)\n",
"\n",
"import petsc4py, sys\n",
"petsc4py.init(sys.argv)\n",
"plot(u_h)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Switching to a storage based on the PETSc solver package and solving the\n",
"sysyem using the dune-fem bindings"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"space = create.space(\"lagrange\", grid, dimrange=1, order=2, storage='petsc')\n",
"scheme = create.scheme(\"galerkin\", a == b, space=space,\n",
" parameters={\"petsc.preconditioning.method\":\"sor\"})\n",
" parameters={\"petsc.preconditioning.method\":\"sor\"})\n",
"# first we will use the petsc solver available in the `dune-fem` package (using the sor preconditioner)\n",
"u_h = space.interpolate(initial, name='u_h')\n",
"u_h_n = u_h.copy(name=\"previous\")\n",
"scheme.model.dt = 0.05\n",
"evolve(scheme, u_h, u_h_n)\n",
"plot(u_h)\n",
"\n",
"plot(u_h)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Implementing a Newton Krylov solver using the binding provided by petsc4py"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import petsc4py, sys\n",
"petsc4py.init(sys.argv)\n",
"from petsc4py import PETSc\n",
"ksp = PETSc.KSP()\n",
"ksp.create(PETSc.COMM_WORLD)\n",
......@@ -231,32 +434,46 @@
"ksp.getPC().setType(\"icc\")\n",
"\n",
"class Scheme3:\n",
" def __init__(self, scheme):\n",
" self.model = scheme.model\n",
"\n",
" def solve(self, target=None):\n",
" res = target.copy(name=\"residual\")\n",
" sol_coeff = target.as_petsc\n",
" res_coeff = res.as_petsc\n",
"\n",
" n = 0\n",
" while True:\n",
" scheme(target, res)\n",
" absF = math.sqrt( res_coeff.dot(res_coeff) )\n",
" if absF < 1e-10:\n",
" break\n",
" matrix = scheme.assemble(target).as_petsc\n",
" ksp.setOperators(matrix)\n",
" ksp.setFromOptions()\n",
" ksp.solve(res_coeff, res_coeff)\n",
" sol_coeff -= res_coeff\n",
" n += 1\n",
" def __init__(self, scheme):\n",
" self.model = scheme.model\n",
" def solve(self, target=None):\n",
" res = target.copy(name=\"residual\")\n",
" sol_coeff = target.as_petsc\n",
" res_coeff = res.as_petsc\n",
" n = 0\n",
" while True:\n",
" scheme(target, res)\n",
" absF = math.sqrt( res_coeff.dot(res_coeff) )\n",
" if absF < 1e-10:\n",
" break\n",
" matrix = scheme.assemble(target).as_petsc\n",
" ksp.setOperators(matrix)\n",
" ksp.setFromOptions()\n",
" ksp.solve(res_coeff, res_coeff)\n",
" sol_coeff -= res_coeff\n",
" n += 1\n",
"\n",
"u_h.interpolate(initial)\n",
"scheme3_cls = Scheme3(scheme)\n",
"evolve(scheme3_cls, u_h, u_h_n)\n",
"plot(u_h)\n",
"\n",
"plot(u_h)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the petsc4py bindings for the non linear KSP solvers from PETSc"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def f(snes, X, F):\n",
" inDF = space.petscFunction(X)\n",
" outDF = space.petscFunction(F)\n",
......@@ -281,7 +498,6 @@
" snes = PETSc.SNES().create()\n",
" snes.setMonitor(lambda snes, i, r:print())\n",
" snes.setFunction(f, res_coeff)\n",
" # snes.setUseMF(True)\n",
" matrix = self.scheme.assemble(target).as_petsc\n",
" snes.setJacobian(Df, matrix, matrix)\n",
" snes.getKSP().setType(\"cg\")\n",
......
from scipy.sparse import coo_matrix
import numpy as np
def is_symmetric(m):
"""Check if a sparse matrix is symmetric
Parameters
----------
m : array or sparse matrix
A square matrix.
Returns
-------
check : bool
The check result.
"""
if m.shape[0] != m.shape[1]:
raise ValueError('m must be a square matrix')
if not isinstance(m, coo_matrix):
m = coo_matrix(m)
r, c, v = m.row, m.col, m.data
tril_no_diag = r > c
triu_no_diag = c > r
if triu_no_diag.sum() != tril_no_diag.sum():
return False
rl = r[tril_no_diag]
cl = c[tril_no_diag]
vl = v[tril_no_diag]
ru = r[triu_no_diag]
cu = c[triu_no_diag]
vu = v[triu_no_diag]
sortl = np.lexsort((cl, rl))
sortu = np.lexsort((ru, cu))
vl = vl[sortl]
vu = vu[sortu]
check = np.allclose(vl, vu)
return check
def eigen(A):
from scipy.sparse.linalg import eigs
matrix = A.transpose(copy=False)
evals_large, evecs_large = eigs(matrix, 1, which='LM')
evals_small, evecs_small = eigs(matrix, 1, sigma=0, which='LM')
print('largest evalue:', evals_large[0], ', smallest evalue:', evals_small[0], ', condition number:', evals_large[0]/evals_small[0])
#try:
# evals_small, evecs_small = eigs(matrix, 1, which='SM')
# print(" ", evals_large, evals_small, evals_large[0]/evals_small[0])
#except:
# print("non transposed version failed!")
def compare(scheme1, scheme2, uh):
# print(m1.todense())
# print(m2.todense())
m1 = scheme1.assemble(uh)
m2 = scheme2.assemble(uh)
if not isinstance(m1, coo_matrix):
m1 = coo_matrix(m1)
if not isinstance(m2, coo_matrix):
m2 = coo_matrix(m2)
#for row, col, value in zip(m1.row, m1.col, m1.data):
# print( "({0}, {1}) {2}".format(row, col, value) )
#for row, col, value in zip(m2.row, m2.col, m2.data):
# print( "({0}, {1}) {2}".format(row, col, value) )
m12 = coo_matrix(m1 - m2)
#for row, col, value in zip(m12.row, m12.col, m12.data):
# print( "({0}, {1}) {2}".format(row, col, value) )
print(m12)
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Nonvariational Test Program"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
"collapsed": true
},
"outputs": [],
"source": [
"from argparse import ArgumentParser\n",
"from dune.grid import cartesianDomain\n",
"from math import sqrt\n",
"from norms import compute_norms\n",
......@@ -16,52 +22,133 @@
"import dune.create as create\n",
"import os.path,sys\n",
"import pickle\n",
"import time\n",
"\n",
"problemName = 'laplace' #fullHessian, varyingCoeff, advection, nonD, nonD2, nonSmooth\n",
"schemeName = 'nvdg' #h1H, h1D2, l2H, l2D2, mu, feng, var\n",
"solutionName = 0# 0=smooth, 1=rough, 2=exp\n",
"use_cg = False #True\n",
"polOrder = 2\n",
"\n",
"import time"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start off by selecting the problem and method used, as well as the solution type and polynomial order. These options can be exchanged with the commented out options if desired."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"problemName = 'laplace' # advection2, nonD\n",
"schemeName = 'nvdg' # h1H, h1D2, l2H, l2D2, mu, feng, var\n",
"solutionName = 0 # 0=smooth, 1=rough, 2=exp\n",
"polOrder = 2 # or 1, 3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We set up a simple Cartesian grid with default storage and solution `uh`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"grid = create.grid(\"ALUSimplex\", cartesianDomain([0, 0], [1, 1], [4, 4]))\n",
"space = create.space(\"Lagrange\", grid, dimrange=1, order=polOrder, storage='fem')\n",
"uh = space.interpolate([0], name=\"solution\")\n",
"\n",
"uh = space.interpolate([0], name=\"solution\")"
]
},
{
"cell_type": "markdown",
"metadata": {},