28.2 Solution to the Poisson Equation

28.2.1 Introduction

In this section we will discuss numerical methods for solving the Poisson equation in Cartesian coordinates. If you recall section 15, the Poisson equation is the underlying PDE that we need to solve in order to derive the velocity profiles, e.g., in rectangular channels. We have already derived the velocity profiles in rectangular channels for Poiseuille flow in section 16.4. One of the most commonly found equations for calculating these profiles is Eq. 16.53. Obviously, now that we have studied various numerical methods, we want to apply this knowledge to finding a numerical solution to the Poisson equation (see Eq. 15.14) in rectangular channels which is given by

$\frac{{\partial }^2{v}_x}{\partial z^2}+\frac{{\partial }^2{v}_x}{\partial y^2}=\frac{1}{\eta }\frac{dp}{dx}$

when neglecting the volume forces k_x. Obviously, the right-hand side of Eq. 28.1 is a constant. We will now apply a finite-distance scheme to the left-hand side. For this we use the Taylor series expansion for the second (partial) derivative, which is given by Eq. 4.10 to approximate the two partial derivatives as

$\begin{array}{l}\frac{{\partial }^2{v}_x}{\partial z^2}=\frac{d^2{v}_x}{\partial z^2}=\frac{v_x\left(y,z+\Delta z\right)+v_x\left(y,z-\Delta z\right)-2v_x\left(y,z\right)}{{\left(\Delta z\right)}^2}\hfill \\ \frac{{\partial }^2{v}_x}{\partial y^2}=\frac{d^2{v}_x}{d{y}^2}=\frac{v_x\left(y+\Delta y,z\right)+v_x\left(y-\Delta y,z\right)-2v_x\left(y,z\right)}{{\left(\Delta y\right)}^2}\hfill \end{array}$

Here we have made use of the fact that the finite change along the y-axis is independent of the change along the z-axis. If we now assume the finite differences to be identical in both directions, i.e., ∆y = ∆z = h, we can rewrite Eq. 28.1 as

$\begin{array}{ll}{v}_x\left(y,z+h\right)+v_x\left(y,z-h\right)-2v_x\left(y,z\right)+v_x\left(y+h,z\right)+v_x\left(y-h,z\right)-2v_x\left(y,z\right)\hfill & =\frac{h^2}{\eta }\frac{dp}{dx}\hfill \\ v_x\left(y,z+h\right)+v_x\left(y,z-h\right)+v_x\left(y+h\right)+v_x\left(y-h\right)-4v_x\left(y,z\right)\hfill & =\frac{h^2}{\eta }\frac{dp}{dx}\hfill \end{array}$

or in standard notation

$\begin{array}{l}{F}^{\left(y,z+1\right)}+F^{\left(y,z-1\right)}+F^{\left(y+1,z\right)}+F^{\left(y-1,z\right)}-4F^{\left(y,z\right)}=\frac{h^2}{\eta }\frac{dp}{dx}\hfill \\ F^{\left(y,z\right)}=\frac{1}{4}\left(F^{\left(y,z+1\right)}+F^{\left(y,z-1\right)}+F^{\left(y+1,z\right)}+F^{\left(y-1,z\right)}-\frac{h^2}{\eta }\frac{dp}{dx}\right)\hfill \end{array}$

This is a second-order method from which we can expect decent numerical stability. Eq. 28.2 tells us that we are able to calculate the value F^(y,z) from the neighboring values in the positive and negative y-direction and z-direction, respectively, and the constant on the right-hand side.

28.2.2 A Simple 4 × 4 Mesh

As a first example, we consider the simple example of a microfluidic channel with square cross-section which we divide into a simple 4 × 4 mesh as shown in Fig. 28.1. We can apply Eq. 28.2 for each element of the mesh and thus we have a total of 16 equations for 16 unknowns. Obviously, the microfluidic channel has boundaries on all sides for which we know the values v_x. In the case of no-slip boundary conditions, the velocity is zero at these points. Therefore we have

$\begin{array}{l}{v}_x\left(1,1\right)=v_x\left(1,2\right)=v_x\left(1,3\right)=v_x\left(1,4\right)=v_x\left(2,1\right)=v_x\left(2,4\right)=\hfill \\ v_x\left(3,1\right)=v_x\left(3,4\right)=v_x\left(4,1\right)=v_x\left(4,2\right)=v_x\left(4,3\right)=v_x\left(4,4\right)=0\hfill \end{array}$

which leaves us with only the unknowns v_x (2, 2), v_x (2, 3), v_x (3, 2), and v_x (3, 3). We therefore have to solve a system of four linear equations for four unknowns. This can easily be done using any method discussed in section 25.3. In the following, we will use the sequential dense LU-decomposition (see section 25.3.6) and the two functions DoLUDecompositionSequentialDensePivot and DoLUSubstitutionDensePivot that we implemented in listing 25.10 and listing 25.11.

28.2.3 Sequential Dense LU-Decomposition

Rewriting to the Correct Format. Before turning to a numerical solution, we first rewrite Eq. 28.2 such that we can extract the matrix A and the right-hand side $\overrightarrow{b}$ . For this we find

$F^{\left(y,z+1\right)}+F^{\left(y,z-1\right)}+F^{\left(y+1,z\right)}+F^{\left(y-1,z\right)}-4F^{\left(y,z\right)}=\frac{h^2}{\eta }\frac{dp}{dx}$

As we can see, the right-hand side is a constant, which is why we can simply normalize it to 1. We then find the following four equations

$\begin{array}{l}{v}_x\left(2,3\right)+v_x\left(2,1\right)+v_x\left(3,2\right)+v_x\left(1,2\right)-4v_x\left(2,2\right)=1\hfill \\ v_x\left(2,4\right)+v_x\left(2,2\right)+v_x\left(3,3\right)+v_x\left(1,3\right)-4v_x\left(2,3\right)=1\hfill \\ v_x\left(3,3\right)+v_x\left(3,1\right)+v_x\left(4,2\right)+v_x\left(2,2\right)-4v_x\left(3,2\right)=1\hfill \\ v_x\left(3,4\right)+v_x\left(3,2\right)+v_x\left(4,3\right)+v_x\left(2,3\right)-4v_x\left(3,3\right)=1\hfill \end{array}$

Obviously, several of these values are zero because they are boundary values. If we remove these, we obtain

$\begin{array}{l}{v}_x\left(2,3\right)+v_x\left(3,2\right)-4v_x\left(2,2\right)=1\hfill \\ v_x\left(2,2\right)+v_x\left(3,3\right)-4v_x\left(2,3\right)=1\hfill \\ v_x\left(3,3\right)+v_x\left(2,2\right)-4v_x\left(3,2\right)=1\hfill \\ v_x\left(3,2\right)+v_x\left(2,3\right)-4v_x\left(3,3\right)=1\hfill \end{array}$

which are four unknowns and four equations. We construct the following linear system from this

$\begin{array}{lll}\left(\begin{array}{llll}-4\hfill & 1\hfill & 1\hfill & 0\hfill \\ 1\hfill & -4\hfill & 0\hfill & 1\hfill \\ 1\hfill & 0\hfill & -4\hfill & 1\hfill \\ 0\hfill & 1\hfill & 1\hfill & -4\hfill \end{array}\right)\cdot \hfill & \left(\begin{array}{l}{v}_x\left(2,3\right)\hfill \\ v_x\left(2,3\right)\hfill \\ v_x\left(3,2\right)\hfill \\ v_x\left(3,3\right)\hfill \end{array}\right)\hfill & =\left(\begin{array}{l}1\hfill \\ 1\hfill \\ 1\hfill \\ 1\hfill \end{array}\right)\hfill \\ \text{A}\hfill & \overrightarrow{x}\hfill & =\overrightarrow{b}\hfill \end{array}$

General Rule for H × W Meshes. If you look at A in Eq. 28.4, you can see that there is a certain pattern. In general, the matrix for a H × W mesh with H being the number of elements in the y-direction and W being the number of elements in the z-direction, will have a total number of H · W rows and columns. It consists of two types of submatrices of type

${\text{A}}_1=\left(\begin{array}{lllllllll}-4\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ 1\hfill & -4\hfill & 1\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & 1\hfill & -4\hfill & 1\hfill & 0\hfill & \cdots \hfill & 0\hfill & 0\hfill & 0\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill & ⋮\hfill & ⋮\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 1\hfill & -4\hfill & 1\hfill \\ 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill & 1\hfill & -4\hfill \end{array}\right)$

${\text{I}}_{\text{W}}=\left(\begin{array}{lllll}1\hfill & 0\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ 0\hfill & 1\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ 0\hfill & 0\hfill & 1\hfill & \cdots \hfill & 0\hfill \\ ⋮\hfill & ⋮\hfill & ⋮\hfill & \ddots \hfill & ⋮\hfill \\ 0\hfill & 0\hfill & 0\hfill & \cdots \hfill & 1\hfill \end{array}\right)$

where the second submatrix is the identity matrix I_W. The first matrix A₁ is the matrix that implements Eq. 28.3 for column y in row z. Obviously, we still require the entries for the row above y − 1 and below y + 1, which we obtain by placing the two identity matrices before and after the submatrix A₁. The matrix A is built using the following rule

• Initialize an empty matrix with H · W × H · W entries. Set all values to zero.

• Using the row index z and the column index y, place the submatrix A₁ at (y = 1, z = 1) in the matrix A.

• Place an identity matrix I_W as (y = W + 1, z = 1) just next to the submatrix A₁. This matrix accounts for the entry just below the mesh value (y, z). In the first line, there is no value to account for above (y, z), so this completes the first line of the mesh.

• Place the submatrix A₁ at (y = 1, z = W + 1) in the matrix A.

• Place an identity matrix “before” the submatrix A₁ (y = 1, z = W + 1) in the matrix A. This accounts for the value above the mesh entry (y, z).

• Place an identity matrix “after” the submatrix A₁ (y = 2 W + 1, z = W + 1) in the matrix A. This accounts for the value below the mesh entry (y, z).

• Repeat this process until the last line.

• In the last line, omit the identity matrix after submatrix A₁, as there is no value to account for below (y, z = H).

Vector$\overrightarrow{b}$ . As we normalized the right-hand side, the vector $\overrightarrow{b}$ has a value of 1 for each of its W · H entries.

28.2.4 Maple Worksheet

Preparing A and$\overrightarrow{\mathbf{b}}$ . As we can already see from Eq. 28.5, A will be a good matrix for numerical schemes because it has non-zero entries on the downwards-diagonal, which also happen to be the largest entries in each column. Therefore, it is not even required to perform pivoting. Furthermore, the matrix is sparse and has a lot of zero entries off the downwards-diagonal, which makes it good for numerical schemes.

We will quickly write a Maple worksheet which implements the LU-decomposition. For this we will make use of the functions DoLUDecompositionSequentialDensePivot and DoLUSubstitutionDensePivot from CoreFunctions.txt. The listing is shown in listing 28.1. The listing implements the creation of the matrix A (see line 25) using the two helper functions MakeA1Matrix, which creates the submatrix A₁ (see line 3) and the helper function MakeIdentityMatrix, which creates the identity matrix I_W (see line 16). Both functions paste their data directly into the array provided as parameter A. They require the parameter W, which is the number of columns of the mesh. The function MakeMatrix additionally requires the number of rows. It returns the array A, the vector $\overrightarrow{b}$ and the values of H and W .

 1  restart:read "Core.txt": read "CoreFunctions.txt":
 2  
 3  MakeA1Matrix:=proc(A,offseth,offsetw,W) local h,w:
 4   for h from 1 to W do
 5    for w from 1 to W do
 6     if h=w then
 7      A[offseth+h,offsetw+w]: = -4:
 8     elif w=h-1 then
 9      A[offseth+h,offsetw+w]:=1:
10     elif w=h+1 then
11      A[offseth+h,offsetw+w]:=1:
12     end if:
13    end do:
14   end do:
15  end proc:
16  MakeIdentityMatrix:=proc(A,offseth,offsetw,W) local h,w:
17   for h from 1 to W do
18    for w from 1 to W do
19     if h=w then
20      A[offseth+h,offsetw+w]: = 1:
21     end if:
22    end do:
23   end do:
24  end proc:
25  MakeMatrix:=proc(H,W) local A,h,w,b:
26   A:=Matrix(H*W,H*W):
27   b:=Vector(1..H*W):
28   for h from 1 to H do
29    MakeA1Matrix(A,(h-1)*W,(h-1)*W,W):
30    for w from 1 to W do
31     if h>1 then
32      MakeIdentityMatrix(A,(h-1)*W,(h-2)*W,W):
33     end if:
34     if h<H then
35      MakeIdentityMatrix(A,(h-1)*W,(h)*W,W):
36     end if:
37     b[(h-1)*W+w]:=1:
38    end do:
39   end do:
40   return [A,b,H,W]:
41 end proc:

Performing the LU-Decomposition. The second listing we require is shown in listing 28.2. This listing implements the function SolvePoissonLU, which accepts the matrix created by MakeMatrix in addition to scalingfactor, which is the scaling factor applied to the right-hand side. The function extracts A, b, H, and W from the return value of MakeMatrix (see line 2) and prepares the output matrix mout that we use for displaying the results (see line 6). The actual LU-decomposition is performed in line 7 using two functions that we previously exported to CoreFunctions.txt. The loops around line 10 copy the result of the LU-decomposition and LU-substitution into the output array, which is plotted in line 13.

 1  SolvePoissonLU:=proc(mesh,scalingfactor) local h,w,H,W,b,A,solution,mout:
 2   A:=mesh[1]:
 3   b:=mesh[2]:
 4   H:=mesh[3]:
 5   W:=mesh[4]:
 6   mout:=Matrix(H,W):
 7   solution:=DoLUSubstitutionDensePivot(DoLUDecompositionSequentialDensePivot(A),b):
 8   for h to H do
 9     for w to W do
10       mout[h,w]:=scalingfactor*solution[(h-1)*W+w]:
11     end do:
12   end do:
13   return matrixplot(mout);
14  end proc:
15 
16  SolvePoissonLU(MakeMatrix(20,10),-0.25);

Line 16 applies the function SolvePoissonLU to a 20 × 10 mesh using a given prefactor. This prefactor was calculated from Eq. 28.3 using the values of water for the viscosity, a pressure drop of $\frac{dp}{dx}=-0.1$ and a stepwidth of h = 5 µ m. From Eq. 28.3 we learn that the prefactor in this case is

$\begin{array}{ll}prefactor\hfill & =\frac{h^2}{\eta }\frac{dp}{dx}\hfill \\ \hfill & =-\frac{{\left(5\text{μm}\right)}^2}{1\text{mPa}\text{s}}0.1\text{mbar}{\text{mm}}^{-1}\hfill \\ \hfill & =-0.25\text{mm}{\text{s}}^{-1}\hfill \end{array}$

The output created from listing 28.2, line 16 is shown in Fig. 28.2. Assuming a stepwidth of 5 µ m, the overall dimensions of the channel amount to a width of 100 µ m and a height of 50 µ m. This is a problem we have studied before analytically and for which we showed the flow profiles in Fig. 16.11. As you can see, the velocity profile looks a lot like the profiles we have derived in section 16.4. However, if we look in detail, we can see that our numerical solution is slightly off as we obtain center velocities of around 3.35 mm s⁻¹ which is bigger than the values we found for the analytical solution. The correct value for the center velocity should be around 2.8 mm s⁻¹. This incorrectness is due to the fact that the mesh we choose is still reasonably coarse, i.e., the increments are too big. This makes our results incorrect. However, you will see that increasing the values for H and W in MakeMatrix will quickly lead to very long computational time.

Making the Solution More Exact. As stated, in order to make our solution more exact, we would need to speed up the computation somehow. We can do this by using a less computationally expensive method (e.g., SOR) or by using a faster implementation of our methods.

28.2.5 Summary

In this section we have solved the Poisson equation that describes the velocity profile in a rectangular cross-section microfluidic channel using LU-decomposition, making use of the methods we previously implemented. We found that the methods provide us with solutions, albeit on very rough meshes. The inherent disadvantage of LU-decomposition is the fact that it is a somewhat slow method, especially when implementing pivoting.

28.3 Solution to the Poisson Equation Using SOR

28.3.1 Introduction

We will now provide a slightly modified version of SolvePoissonLU using SOR (see section 25.3.4) instead of LU-decomposition. As already stated, one the strengths of LU-decomposition is the fact that we can reuse solutions, i.e., the decomposed matrices L and U for different right-hand side $\overrightarrow{b}$ . However, in our case, this is not really necessary, as the right-hand side of the Poisson equation is a constant (see Eq. 28.1).

For this we will use the function DoSOR, which we also exported to CoreFunctions.txt. The modified implementation is shown in listing 28.3.

 1  SolvePoissonSOR:=proc(mesh,scalingfactor) local h,w,H,W,b,A,solution,mout:
 2   A:=mesh[1]:
 3   b:=mesh[2]:
 4   H:=mesh[3]:
 5   W:=mesh[4]:
 6   mout:=Matrix(H,W):
 7   solution:=DoSOR(A,b,Vector(1..W*H),1.6,100):
 8   for h to H do
 9    for w to W do
10     mout[h,w]:=scalingfactor*solution[(h-1)*W+w]:
11    end do:
12   end do:
13   return matrixplot(mout);
14  end proc:
15 
16  SolvePoissonSOR(MakeMatrix(40,20),-0.0625);

The main difference between the implementation for SOR and LU-decomposition is the fact that we require some additional parameters for the function DoSOR. We use an empty vector as an initial guess and assume a relaxation parameter ω = 1.6. This value is chosen arbitrarily and may be further optimized. In addition, we set a total number of iterations of 100. As the function is faster than LU-decomposition, we use a finer mesh using a stepwidth h = 2.5 µ m, which also changes our scaling factor according to Eq. 28.3. You can see that we now use a mesh of 40 × 20 elements that is already significantly more refined. The output produced by the function SolvePoissonSOR is shown in Fig. 28.3. As you can see, the output is close to the exact solution shown in Fig. 16.11. The maximum velocities obtained from line 16 are in the range of 2.9 mm s⁻¹, which is definitely close to the correct value.

Disadvantages of the SOR Implementation. As stated, SOR allows us to work with more refined meshes, but it still is not a fast implementation. Also, the fact that we have to indicate the total number of iterations is disadvantageous. Ideally, we would like to give the required precision as parameter and let the algorithm proceed until either the accuracy is obtained, i.e., the values from one iteration to the next change for values smaller than the precision, or a given maximum number of iterations has been reached. Obviously, these are implementation details that are not difficult to work out.

Additionally, it seems reasonable to implement SOR in a programming language that allows for faster running code. Maple, being a symbolic high-level programming language, is inherently slower because it is not optimized for running numerical schemes efficiently. This is why we will now reimplement the SOR algorithm in the programming language C.

28.3.2 A SOR Implementation in C

Interaction of Low- and High-Level Languages. We are about to move some of the implementations we made in Maple into a so-called low-level language. This is a very common method. Low-level languages usually run significantly faster, as they are closer to the machine architecture and require less computationally expensive overhead, such as type casting and overflow checking. High-level languages obviously offer significant advantages, as they support the user with many features without which complex state-of-the-art programs could not run, including graphical user interfaces (GUIs), parallel computing, and object-oriented programming.

However, the numerical schemes used for solving ODEs or PDEs do not require this overhead in computation. Most of the time, the program will simply loop through arrays of data, improving the current estimations until a given accuracy is reached. These routines are best implemented in a way that they can be run fast and (ideally) in a parallel manner on a computer. Usually, inputting data into these routines, as well as interpretation of the results of these routines, is best programmed in a high-level programming language, because they provide better interfaces, e.g., for displaying the results.

This separation of implementation is actually very common. Whenever you start the numerical solver of a modern 3D computer aided design (3D-CAD) system, you are likely to see smaller, usually less fancy-looking, interfaces pop up which give you rough information about the calculation in the background. At this point, you know that low-level code is running a numerical scheme, which will be passed back to the main program once the calculation is finished.

28.3.3 Implementation

The implementation of our SOR routine in C is called SolvePoissonSOR and is shown in listing 28.4. It accepts the number of columns W and the number of rows H in the mesh, as well as the required accuracy epsilon and the right-hand side of Eq. 28.1inhomogeneous. The relaxation parameter ω is given as the parameter omega. The initial data of the mesh (which are treated as the initial conditions) is given as array in the parameter data. The array BC is an array that is expected to be of the same dimension as data. If an entry at the location (i, j) in this array is set to the value 1, the algorithm will treat the respective value in data as a boundary condition, i.e., it will not change this value during the calculation. The last parameter log is a string array into which the algorithm will write additional information for the user.

 1  int DLLEXPORTIMPORT SolvePoissonSOR(int W, int H, double epsilon, double inhomogeneous, double omega, double* data, char* BC, char* log){
 2   char isconverged;
 3   double store;
 4   int count=0;
 5   int currenty;
 6   double prefactor1=1-omega;
 7   double prefactor2=0.25*omega;
 8   do{
 9    currenty=0;
10    isconverged=1;
11    for(int y=0;y<H;y++){
12     for(int z=0;z<W;z++){
13      if(BC[currenty+z]!=1){
14       if(y==0) data[currenty+z]=data[currenty+W+z];
15       else if(y==H-1) data[currenty+z]=data[currenty-W+z];
16       else if(z==0) data[currenty+z]=data[currenty+z+1];
17       else if(z==W-1) data[currenty+z]=data[currenty+z-1];
18       else{
19          store=prefactor1*data[currenty+z]+prefactor2*(data[currenty+z-1]+data[currenty+z+1]+data[currenty-W+z]+data[currenty+W+z]-inhomogeneous);
20         if(isconverged && fabs(store-data[currenty+z])>epsilon){isconverged=0;}
21         data[currenty+z]=store;
22       }
23      }
24     }
25     currenty+=W;
26    }
27    count++;
28    if (count==MAXITERATIONS) {
29     sprintf(log, "Reached maximum number of iterations");
30     return count;
31    }
32   }while(!isconverged);
33   sprintf(log, "Converged! Number of iterations: %i",count);
34   return count;
35  }

The function SolvePoissonSOR will first initialize a number of variables. isconverged is a boolean value which will be set to 1 once the algorithm has not changed any value in the mesh by a value greater than epsilon between two iterations. This variable is therefore used as the termination condition of the main loop in line 8. The convergence is checked in line 32. The value store is used as an intermediate storage. The variable currenty is used to speed up accessing the data of row y of the mesh in data. count is the total number of iterations. The two variables prefactor1 and prefactor2 are used to speed up the calculation of the next parameter, according to Eq. 25.16. They are used in line 19 where the next value at location (y, z) is calculated from the last value at the location (y, z) which is stored in data at the location (y · W + z) = (currenty + z), and the current value at location (y, z) which is calculated according to Eq. 28.2.

As stated, the core of the function is the main loop starting from line 8. This loop first resets currenty as we start with the first row. We then assume that the problem is already converged by setting isconverged to 1 (see line 10). Obviously, this is not the case and (as we will see), isconverged is reset within the main loop. We then iterate over all rows and columns starting from line 11. At each position, we check if the value given is a boundary condition (see line 13). If this is the case, we skip this entry in the mesh as we know this value must not change. The next lines check if we are at the edge of the mesh, i.e., if y = 0 or y = H − 1 in y-direction and/or if z = 0 or z = W − 1 in z-direction (see line 14). If this is the case, we are left with an open surface for which we have a Neumann boundary condition. In this case, we simply copy down the value next to the boundary toward the center of the mesh. We will look at one example of such a mesh in just a moment. In all of these cases, we do not require a new calculation for the current value.

The code starting from line 18 is only called for values which are neither boundary conditions, nor at the edge of the mesh. Here, we perform SOR for calculating a better approximation of the value at the location (y, z). We store this better approximation in store because we need to perform a comparison with the last value to find out if this value has changed to a degree larger than the tolerance given by epsilon (see line 20). If this is not the case, the flag isconverged is reset. We only perform the comparison operation if the flag was not already reset previously. After the comparison, we store the newly calculated value in the array data (see line 21). At the end of the loop, we increase the current array index currenty (see line 25). This index pointer allows us to access the current row in the array data more quickly.

Before going into the next loop, the function checks if the maximum number of iterations has been reached (see line 28). If this is the case, the function was not able to find a converged solution to the given problem within the given maximum number of iterations. In the file NeptunLib.c this value is set to 10 000. This value is reasonably high, but there may be applications where higher numbers are required, especially for good precision. If the function has not reset the flag isconverged the function exits by logging the required number of iterations.

28.3.4 Calling SolvePoissonSOR From Maple

If you are somewhat experienced with programming, you may have realized that SolvePoissonViaSORinC has the export directive DLLEXPORTIMPORT, which means that this function is exported to a DLL. We can call this function from Maple using a suitable import directive. The listing is shown in listing 28.5. The Maple function define_external is used to import functions from a DLL (see line 1). It requires the name of the function to import, the name of the DLL, as well as the arguments expected and their respective data type. We pass W and H as 4-byte integer values. epsilon, inhomogeneous, and omega are passed as 8-byte floating point value. The two arrays data and BC are passed as ARRAY, for which we must inform Maple about the size of the array. The fact that we set the names of the two arrays in single quotation marks means that we expect Maple to pass a pointer of the arrays to the DLL. This way, all modifications on the array done by the DLL are performed on the actual data, i.e., Maple will later see these modifications as well. The last argument passed is the string log which we pass (again) as a pointer to a string array of fixed length with 256 characters. The return value we assume is, again, a 4-byte integer value.

 1  C_SolvePoissonSOR:=define_external('SolvePoissonSOR',LIB="NeptunLib.dll",W::(integer[4]),H::(integer[4]),epsilon::(float[8]),inhomogeneous::(float[8]),omega::(float[8]),'data'::(ARRAY(1.. W*H,float[8])),'BC'::(ARRAY(1..W*H,boolean[1])),'log'::(string[256]),RETURN::(integer[4]));
 2
 3  MakeBCIV_rectangular:=proc(W,H,data) local BC,z,y:
 4   BC:=Matrix(W,H,datatype=boolean[1]):
 5   for y to H do
 6    for z to W do
 7     if y=1 or z=1 or y=H or z=W then
 8      BC[z,y]:=1:
 9     end if:
10    end do:
11   end do:
12   return BC:
13  end proc:
14
15  SolvePoissonSORinC:=proc(W,H,MakeBCIV,epsilon,inhomogeneous,omega) local data,BC,z,y,log:="":
16   data:=Matrix(W,H,datatype=float[8]):
17   BC:=MakeBCIV(W,H,data):
18   C_SolvePoissonSOR(W,H,epsilon,inhomogeneous,omega,data,BC,'log'):
19   print(log);
20   return matrixplot(data);
21  end proc:
22
23  SolvePoissonSORinC(40,20,10E-5,-0.0625,1.8);
24  SolvePoissonSORinC(200,100,10E-5,-0.0025,1.8);

A detailed explanation of all the features of the define_external and native calls from Maple to DLLs are beyond the scope of this book. The interested reader is referred to the excellent online documentation of Maple or to the Maple Advanced Programming Guide[1].

The rest of the listing is quite simple. We first define the function SolvePoissonSORinC which again accepts the number of columns W and rows H of our mesh, as well as the required accuracy epsilon, the right-hand side inhomogeneous and the relaxation parameter omega (see line 15). It then initializes the two arrays data and BC that hold the data and the boundary conditions, respectively. The boundary conditions are actually created using the function passed as the argument MakeBCIV, which in our examples will be the function MakeBCIV_rectangular. This function makes a mesh of correct size with fixed boundary conditions at the outer edges of the mesh by setting all values for (1, z), (H, z), (y, 1), and (y, W ) to one, which tells the C-function SolvePoissonSOR to treat these values as the boundary values. In line 18 we call the external function passing the variable log as a pointer.

After the function completes, we print out log and display the results. The function passed as MakeBCIV also received a pointer to data, as it is also responsible for setting the initial values.

Line 23 repeats the example we already used in listing 28.3. You will see that the function completes almost instantaneously. It will require around 60 iterations to achieve an accuracy of 10 × 10⁻⁵. You will see that the results look significantly better, and that the profile obtained is very close to the analytical result with a center velocity of around 2.58 mm s⁻¹. Given that the function completes so fast, we can now use an even finer mesh with 200 × 100 entries (see line 24). Even this large number of mesh points will not require more than a few seconds on an ordinary computer. The total number of iterations will be around 1000 and strongly dependent on the value of omega. Depending on the mesh and the right-hand side, this value allows significant fine-tuning of the performance of the function.

The results of line 24 are shown in Fig. 28.4. As you can see, the obtained maximum center velocity is 2.8 mm s⁻¹, which is the correct value.

Assessment. We have just seen a very striking example of the speed benefit that an optimized algorithm in a fast programming language can obtain. Please note that the functions we used are essentially the same as the ones we implemented in Maple before. However, C-code runs significantly faster being closer to machine language than any high-level programming language can ever be. Interfacing these functions from Maple is easy and allows using the visualization tools we have available in Maple, which are somewhat more challenging to implement in C.

The function we have just implemented is generic, i.e., we can use it for a wide range of applications. Analytically, we have written a solver for Poisson and Laplace equations, as the Laplace equation is just the homogeneous version of the Poisson equation. We have demonstrated its use for solving the Poisson flow in rectangular cross-sections, i.e., for calculating the velocity profiles. However, the solver can be used for any type of problem as long as the underlying PDE is a Laplace or a Poisson equation.

28.3.5 Applying the Solver

Couette Flow. Having written such a useful solver, we now want to apply it to various application scenarios, checking if the results obtained with it correspond to the results found analytically. The first case we will check is the Couette flow (see section 15.2). Here, we will give the solver different boundary conditions to work with

• The top of our mesh is expected to move at a constant velocity, i.e., 5 mm s⁻¹.

• The left and the right boundary are open surfaces, i.e., we expect Neumann boundary conditions (see section 3.1.7.2). If you remember our implementation in listing 28.4, we can do so by simply letting the values of BC be zero on the right and left boundary.

• The bottom of our mesh is expected to be static, i.e., the values in BC should be zero.

As you can see, the only thing we need to change for this application scenario is passing a different function for making the boundary conditions as parameter MakeBC to SolvePoissonSORinC. For the Couette flow, the function used is called MakeBCIV_Couette (see listing 28.6). It sets the boundary values for the top and the bottom edge of the mesh, but leaves values of 0 on the right and the left edge. Furthermore, it sets the initial value of the upper edge to 5. The output created from line 15 is shown in Fig. 28.5a. As expected, we see the linear Couette flow profile. Please note that the left and the right edge of the mesh fulfill the Neumann boundary conditions as expected.

 1  MakeBCIV_Couette:=proc(W,H,data) local BC,z,y:
 2   BC:=Matrix(W,H,datatype=boolean[1]):
 3   for y to H do
 4     for z to W do
 5      if y=1 then
 6       BC[z,y]:=1:
 7       data[z,y]:=5:
 8      elif y=H then
 9       BC[z,y]:=1:
10      end if:
11     end do:
12    end do:
13    return BC:
14  end proc:
15  SolvePoissonSORinC(40,40,MakeBCIV_Couette,10E-5,0,1.8);
16
17  MakeBCIV_Elliptical:=proc(W,H,data) local BC,z,y:
18   BC:=Matrix(W,H,datatype=boolean[1]):
19   for y to H do
20     for z to W do
21      if ((z-W/2)/(W/2))^2 + ((y-H/2)/(H/2))^2 > 0.9 then
22        BC[z,y]:=1:
23      end if:
24     end do:
25   end do:
26   return BC:
27  end proc:
28  SolvePoissonSORinC(40,40,MakeBCIV_Elliptical,10E-5,-0.0625,1.8);

Hagen-Poiseuille Flow. In the next step, we will apply the solver to a circular cross-section that we have derived analytically as the Hagen-Poiseuille flow (see section 16.3). For this, we use the function MakeBCIV_Elliptical (see listing 28.6), which simply sets all values outside of the ellipse enclosed by the mesh with width W and height H to 1, therefore making them boundary conditions. The output produced by line 28 is shown in Fig. 28.5b. We observe the familiar Hagen-Poiseuille parabola. The maximum velocity in the center is around 5.7 mm s⁻¹, which is the same value we find analytically using Eq. 16.20.

28.3.6 Summary

In this section, we have implemented our first real solver using SOR in the C programming language. As we have seen, the solver outperforms all functions we have currently written, which is not surprising given that it is implemented in a very fast low-level language. We have used this solver to solve numerous cases of microfluidic flows for which we have derived analytical solutions in section 15. As we have seen, such numerical solvers can be very handy to solve cases for which analytical solutions are difficult to obtain.