PDE Solver DY

This PDE solver is part of a project for the Advanced Programming course. It is concerned with the solution of the homogenous 2D heat transfer equation (LaPlace’s equation) for the steady-state case. It is based on C++ and is currently supported on Linux OS.

Installation

Prerequisites

The program assumes that the following tools and libraries are setup properly on the Linux OS:

• Git
• Make
• CMake
• Python (v3)

sudo apt-get install python3-dev

• Python pip

sudo apt install python3-pip

• Python numpy

sudo apt-get install python-numpy

• Matplotlib

pip install matplotlib

• Python GUI

sudo apt-get install python3-tk

• Install Eigen libray (v3.0)

sudo apt install libeigen3-dev

Git cloning

Clone the git repository using the following code

git clone https://gitlab.lrz.de/00000000014AE223/pde-solver-dy.git

Building and Running

Inside the /build folder, run cmake ../src command. Then inside the /build folder run make command. Finally, open the executable file using ./pde_solver command inside the /build folder.

Sprint 1 (v1.0)

Introduction

The solver implementation in this version is limited to a simple rectangular 2D domain in the cartesian coordinate system (i.e. x and y). The problem definition is summarized in the figure below.

In this version, only Dirichlet boundary conditions are allowed for continuous homogeneous materials (only a single material for the whole domain is expected).

The user inputs and their corresponding restrictions are:

1. Domain specification (double, double)
   • Major dimensions: the width and height of the rectangle (W and H, respectively). Must be greater than 0 and with compatible units.
2. Mesh specification (int, int)
   • Number of nodes in each direction (n_x and n_y), step size is calculated accordingly. Must be integers and greater than 1.
3. Boundary conditions (double, double, double, double)
   • Currently, only Dirichlet BCs are supported. Therefore, 4 temperature values are required according to the diagram above (T1, T2, T3, and T4). Temperature values are in degree Celsius.

The 5-point finite difference scheme, demonstrated in the figure below, is used to discretize the 2D heat transfer equation. The resulting linear system is solved using the LU-factorization method.

Outputs

After execution of the simulation with all input parameters, a CSV file (results.csv) can be found in the /results folder of the project containing all of the information required. The format of this comma-separated file is as follows:

• Node number, 1st coordinate, 2nd coordinate, BC (1/0), T_sim, T_analytic, error

Also, if dependencies are installed correctly, a 3D plot for the simulated temperature values will be generated. Finally, the code outputs to the console the maximum, minimum, and average values of the error percentage. (if the test case option is selected).

Test case

The simple unit test implemented for the testing of the numerical simulation accuracy is defined as special case of the general problem presented before. Only two temperature values (T1 and T2) are expected as BCs, and the inhomogeneous case is not supported. This special case is presented below:

For this problem, the analytical solution for the temperature value (T) for any point inside the domain is given by:

Results from the analytical solution is compared with the grid values obtained from the numerical solution. The error percentage for each node is calculated and the maximum, minimum, and average values are reported to the console.

Example usage

In the /build folder, use the command make to build the executable file. Then, use the command ./pde_solver.

Test case 1:

Set Width to: 1

Set Height to: 1

Set Nodes in X direction to: 11

Set Nodes in Y direction to: 11

Set Run unit test to: 1

Set Temperature T1 to: 0

Set Temperature T2 to: 100

Check ../results/results.csv file with the details of the results.

Basic case

Set Width to: 2

Set Height to: 5

Set Nodes in X direction to: 25

Set Nodes in Y direction to: 18

Set Run unit test to: 0

Set Temperature T1 to: -20

Set Temperature T2 to: 100

Set Temperature T3 to: 50

Set Temperature T4 to: 0

Check ../results/results.csv file with the details of the results.

Sprint 2 (v2.0)

What's new?

The solver from Sprint 1 (v1.0) is now extented with more functionality using the OOP structure. The following is a summary of the newly added features:

1. New coordinate system, Polar coordinates. Now it's possible to choose from the Cartesian (eg. square, rectangle) or the Polar coordinates (eg. circle or oval).
2. Neumann Boundary conditions (d_u/d_n = constant). The user can now choose either one of BC types (Dirichlet or Neumann) for each boundary in the problem. Note: at least one Dirichlet boundary condition should be present to arrive at a unique solution.
3. Inhomogenous problem (i.e. Poisson's equation).
4. New numerical solvers.
   • Special LU scheme to utilize the sparse structure of the matrix (much faster than classical)
   • Gauss-seidel (iterative solver requiring max number of iterations, relaxation factor, and max value of the residual norm)

Main user inputs

An updated list for the main user inputs and the corresponding restrictions:

1. Coordinate system (int - Cartesian = 1, Polar = 2)
2. Domain type (int - rectangle, square [only one input required], circle, or oval). Each coordinate system restricts the possible domain types to choose from.
3. Domain specification (double, double)
   • Major dimensions: the width and/or height of the rectangle/square, or the radius (for a circle), or minor and major radii (for oval). Must be greater than 0 and with compatible units.
4. (for inhomogenous case) Value and location of heat source/sink (double, double, double)
   • The location of the heat source/sink must be inside the specified domain. The provided location is then applied to the nearest node in the mesh.
5. Mesh specification (int, int)
   • Number of nodes in each direction (n_x and n_y), step size is calculated accordingly. Must be integers and greater than 1.
6. Boundary condition types (int, int, int, int)
   • Only two possible options (Dirichlet or Neumann). At least one Dirichlet BC should be specified for a unique solution.
7. Boundary condition values (double, double, double, double)
   • Four temperature or heat source/sink values are required for Dirichlet and Neumman BCs, respectively, according to the diagram below (BC1, BC2, BC3, and BC4 for Cartesian domain OR _BC1 and BC2 for Polar case).
8. Solver (int)
   • Current available solvers: classical LU factorization, LU-Sparse, Gauss-Seidel

Outputs

After execution of the simulation with all input parameters, a CSV file (results.csv) can be found in the /results folder of the project containing all of the information required. The format of this comma-separated file is as follows:

• Node number, 1st coordinate, 2nd coordinate, BC (1/0), T_sim, T_analytic (if test case chosen), error (if test case chosen)

Also, if dependencies are installed correctly, a 3D plot for the simulated temperature values will be generated. Finally, the code outputs to the console the maximum, minimum, and average values of the error percentage. (if the test case option is selected).

Extending the functunality

Due to the structure of the current implementation, a new feature can be added to extend the functionality of the solver. The following is a list of possible extensions that can be carried out to add more features to the current version:

1. Beyond Elliptic PDE Heat equations
   • Although the main purpose of this solver is to solve the 2D Poisson's heat equation, it has been implemented with the possibility to solve other types of PDE models in mind.
   • This can be done by creating a child class of the Initiation base class and providing the required paramters to the Initiation constructor.
   • The code is already equipped with three PDE Types (Elliptic, Hyperbolic, and Parabolic) with the Elliptic case being the default value. The pure virtual getter functions must be declared in the new child class.
2. More coordinate systems
   • The current solver considers Cartesian and Polar systems.
   • This can simply be extended by adding other systems such as cylinderical, spherical, barycentric, or curvilinear to the CoordinateSystem enum class defined in Initiation header.
   • The current coordinate systems are being used as part of the workflow logic in different classes and methods (w.o.l.g. setting inhomogenity values or, generally, boundary conditions in EllipticPDE, grid generation in Mesh, setting up the system matrix and the RHS b vector in Solver, possible shapes and domain sizing in Domain). Since each coordinate system will need special considerations for the aforementioned examples, the current logic need to be extended for the newly added coordinate system.
3. More shapes for each coordinate system
   • To add a new domain shape (eg. triangle, hollow circle/oval, or hollow square), each shpae need to be restricted to only a single coordinate system.
   • If a shape needs to be defined with more than 2 values, the _major_dimensions array size needs to be adjusted accordingly in Domain. In addition, the number of possible boundary conditions need to be defined using the currently implemented shaapes as an example in the Domain constructor. In the Mesh, the logic needs to incorporate how to handle the grid generation of this new shape.
4. Different meshing strategy
   • Structured rectangular mesh (consistent step sizes in all dimensions) is the only option being employed in this version of the solver.
   • To add/define other grid generation schemes such as unstructured mesh or triangle mesh, the MeshType enum class can be extended accordingly. The default option of __mesh_type is the rectangular structured grid. This can be explicity defined by the user and the program can handle this case independently in Mesh constructor.
   • Any new scheme must always abide by the format of the __mesh matrix. Any further information that needs to be stored should be stored in a new matrix/vector (such as the nodes and sides numbers of a certain triangular element), this still needs to be a protected class member.
5. More numerical solvers
   • Besides the implemented classical LU factorization, LU-sparse, and Gauss-Seidel solvers, other direct or iterative solvers can be defined such as Jacobi, Gauss elimination, or banded solver to name a few. The current structure of the software suggest that a new solver has the following features:
     • Private/Protected parameters where the data required for the solver is stored such as a vector for the right-hand side of the system, a matrix for the left-hand side of the system, the maximum number of iterations of a matrix-free solver and the error tolerance of an iterative solver.
     • Private/Protected methods that set the parameters defined previously as well as intermediate methods necessary for the solution of the system.
     • A public solve method that solves the system and sets it in the solution parameter of the Solver class. The signature of this method must be: void solve(const std::unique_ptr& pde, const std::unique_ptr& mesh)
     • A default constructor that takes the parameter test_case as input. The suggested signature of this method is: new_solver(const bool test_case)
     • Public methods that define the getters desired for the new solver.

Example Case: test

Set Test case to: 1

Set shape to: 1 (Square)

Set length to: 1

Set Temperature T1, T2, T3 to: 25

Set Temperature T4 to: 100

Set no. of nodes for X direction to: 115

Set no. of nodes for Y direction to: 93

Set solution method to: 1 (LU sparse)

Check ../results/results.csv file with the details of the results.

Example Case: Neumann boundary condition

Set Test case to: 0

Set Coordinate System to: 1 (Cartesian)

Set Homogeneous to: 1

Set shape to: 1 (Square)

Set length to: 1

Set Bottom boundary type to: 2 (Neumann)

Set Left boundary type to: 1 (Dirichlet)

Set Right boundary type to: 1 (Dirichlet)

Set Top boundary type to: 1 (Dirichlet)

Set Bottom heat flux to: 0 (isolated)

Set Left Temperature to: 75

Set Right Temperature to: 50

Set Top Temperature to: 100

Set no. of nodes for X direction to: 5

Set no. of nodes for Y direction to: 5

Set solution method to: 2 (Gauss Seidel)

Set maximum number of iterations to: 100

Set residual error tolerance to: 0.00001

Set relaxation parameter to: 1.5

Check ../results/results.csv file with the details of the results.

Sprint 3 (v3.0)

What's new?

After studying the performance of the code, optimization was carried out for the bottleneck functions/methods. A complete performance study and the new results after optimization can be found in performance.md. As a teaser, the program is now at least 4.8x faster than the original code (depends on the test case as well as other parameters, see a complete discussion in performance.md)!