Solution to a 1D Laplacian equation using FDM

This repository proves a 1D solution to the Laplacian equation using PETSc. It also supplies a build in Nix development environment using the flake.nix file.

The main purpose of this repository is to provide a simple learning example for interfacing with the C bindings in the PETSc library.

How to Compile the simulation

Included in this directory is a makefile. This make file is about as simple as possible for a PETSc document and is reproduced in its entirety below:

include ${PETSC_DIR}/lib/petsc/conf/variables
include ${PETSC_DIR}/lib/petsc/conf/rules

This uses the users installation of PETSC_DIR to include some standard variables and build targets for a PETSC application. Alternatively we could specify each included library in our makefile.

To compile all we need to do is:

make main

Which will compile the main.c file and link the appropriate PETSc libraries

How to run the simulation

This is a uniprocessor example. In general we could run a PETSc progam as follows:

mpiexec -n 8 ./petsc_program_name petsc_options

Therefore, the simplest way of running the main.c file in this repository is:

mpiexec -n 1 ./main

In this example, on line 35 we have checked if the user has provided the value of n where n is the dimension of the problem. That means we can supply this value at runtime as follows:

mpiexec -n 1 ./main -n 100

This will override the default value of n = 10 and run the simulation with n = 100

Visualizing the solution

Generally, in scientific simulations like this, we can visualize the solution using a visualization software like Paraview. In this case, sicne we have a simple 1D example, it doesn’t make much sense to use Paraview.

Instead we can print out the vector using a function from PETSc called VecView().

// View the solution vector x
 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_SELF));

We can view the matrix, A in a similar way using the following line of code:

// Print the matrix A
PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF));

This gives the following result (when written in LaTeX):

$$ A = \begin{bmatrix} 2 & -1 & 0 & 0 & \cdots & 0 \\ -1 & 2 & -1 & 0 & \cdots & 0 \\ 0 & -1 & 2 & -1 & \cdots & 0 \\ 0 & 0 & -1 & 2 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & 2 \end{bmatrix} $$

What this program shows

This program essentially solves a linear system of equation, Ax = b where A is the tridiagonal matrix shown above. We can think of this matrix as the discrete Laplacian operator in one dimension. This is used to approximate the second derivative in differential equations.

When using a finite difference discretization, the second derivative at a point $x_i$ can be approximated by:

$$ \frac{d^2 u}{dx^2} ≈ \frac{u_ {i-1} - 2u_i + u_ {i+1}}{h^2} $$

If we let $h=1$ then each row of $Au$ corresponds to the finite difference equation for a point in the domain.

This point of this program is to carry out this approximation (solving the one-dimensionsal Laplacian equation) using a PETSc linear solver, and to compare the result with the known solution. The norm of the error in the solution is then printed to the screen, along with the number of iterations taken by the solver.