View source for Poisson's equation

__FORCETOC__

== Dirtchlet boundary conditions ==
It is fitting that examples start with a simple case, and we will gradually make our way towards more complicated cases with different domain shapes, boundary conditions and approximation types.
Consider the solution of the 2D Poisson equation on a unit square with Dirichlet boundary conditions:
<math>
\begin{align*}
    	\Delta u &= f      &&\text{in } \Omega, \\
    	  u &= 0           &&\text{on } \partial \Omega,
	\end{align*}
</math>
where $u(x,y)$ is the solution to the problem, $\Omega = [0, 1] \times [0, 1]$ and in the case we will consider $f(x,y) = -2\pi^2\sin(\pi x)\sin(\pi y)$, as it makes for 
a simple solution $u(x,y) = \sin(\pi x)\sin(\pi y)$.


First, we construct our domain and discretize it, then find support (neighborhood) for each node.

<syntaxhighlight lang="c++" line>
BoxShape<Vec2d> box(0.0, 1.0);
double dx = 0.01;
DomainDiscretization<Vec2d> domain = box.discretizeWithStep(dx);

int N = domain.size();
domain.findSupport(FindClosest(9));
</syntaxhighlight>

We constructed a box shape, discretized it with a constant step dx (both inside and the boundary) and choose the closest 9 nodes of each node as its support. The next step is to construct our approximation engine. 
Simply put, approximation engines are classes responsible for [[Computation of shape functions|computing shape functions]] for given operator and points. While many different setups are supported in Medusa we will start off simple,
with an approximation that reproduces standard Finite Difference method.

<syntaxhighlight lang="c++" line>
int m = 2
WLS<Monomials<Vec2d>, NoWeight<Vec2d>, ScaleToFarthest> wls(Monomials<Vec2d>::tensorBasis(m));
</syntaxhighlight>

We constructed a [[Moving Least Squares (MLS)| WLS]] engine using monomial tensor basis of order 2. Finally, we write the equations for our problem by directly translating
the mathematical formulation above into code. 

<syntaxhighlight lang="c++" line>
for (int i : domain.interior()) {
    double x = domain.pos(i, 0);
    double y = domain.pos(i, 1);
    op.lap(i) = -2*PI*PI*std::sin(PI*x)*std::sin(PI*y);
}
for (int i : domain.boundary()) {
    op.value(i) = 0.0;
}
</syntaxhighlight>

Here is the plot of the solution $u(x, y)$.

[[File:dirichlet.png|800px]]

The whole example can be found as [https://gitlab.com/e62Lab/medusa poisson_dirichlet_2D.cpp] along with the plot script that can be run by Matlab or Octave [https://gitlab.com/e62Lab/medusa poisson_dirichlet_2D.m].
Same goes for similar cases solved in 1D and 3D and their respective plot scripts.

== Mixed boundary conditions  == 

We can obtain the solution to the above problem by solving it on just a quarter of the original domain and using mixed boundary conditions. By using Dirichlet boundary conditions on the bottom and the left side of the square domain, and Neumann boundary conditions on the right and upper side of the domain we are effectively solving the same problem but on a quarter of the domain. Our problem is now:

<math>
\begin{align*}
    	\Delta u &= f      &&\text{in } \Omega, \\
    	  u &= 0           &&\text{on } \partial \Omega_1 \cup \Omega_3,\\
          \frac{du}{d \b{\hat{n}}} &= 0           &&\text{on } \partial \Omega_2 \cup \Omega_4 
	\end{align*}
</math>

Where $\partial \Omega_1$, $\partial \Omega_2$, $\partial \Omega_3$ and $\partial \Omega_4$ are left, right, bottom and top boundary respectively. The domain is $\Omega = [0, 0.5] \times [0, 0.5]$, and $\frac{du}{d \b{\hat{n}}}$ is the normal derivative (on the boundary).

When constructing approximation engines, there are a lot of different options, you can choose between different basis functions, weights, scale functions and each setup comes with its advantages and disadvantages. Some might work for certain problems and some might not. When writing the examples we tried to present the reader with as many different options as possible. 

Let us solve the problem using Gaussian RBF (Radial Basis Functions) instead of Monomials:

<syntaxhighlight lang="c++" line>
    // shape = 30
    // m = 9
    WLS<Gaussians<Vec2d>, NoWeight<Vec2d>, ScaleToFarthest,
            Eigen::LLT<Eigen::MatrixXd>> wls({9, 30});
</syntaxhighlight>

While the case is the same for the nodes inside the domain, we need to change the boundary conditions:

<syntaxhighlight lang="c++" line>
    int BOTTOM = -3;
    int TOP = -4;
    int LEFT = -1;
    int RIGHT = -2;

    for (int i : (domain.types() == LEFT) + (domain.types() == BOTTOM)) {
        op.value(i) = 0.0;  // Dirichlet boundary conditions on the left and bottom edge of the box
    }
    for (int i : (domain.types() == TOP) + (domain.types() == RIGHT)) {
        // Neumann boundary conditions on upper and right edge of the box
        op.neumann(i, domain.normal(i)) = 0.0;
    }
</syntaxhighlight>

When solving a PDE implicitly a system of equations needs to be solved to obtain the solution. For this, we use [http://eigen.tuxfamily.org/index.php?title=Main_Page Eigen].

<syntaxhighlight lang="c++" line>
    Eigen::BiCGSTAB<decltype(M), Eigen::IncompleteLUT<double>> solver;
    solver.compute(M);
    ScalarFieldd u = solver.solve(rhs);
</syntaxhighlight>

Matrix M and rhs are the sparse matrix and vector in which we wrote (with the help of Medusa) the system of equations that when solved gives us an approximate solution of our problem. We solve the system using an Eigen linear equation solver, in this case, the BiCGStab algorithm with ILUT preconditioning. The solution is plotted below:

[[File:mixed.png|800px]]

The whole example can be found as [https://gitlab.com/e62Lab/medusa poisson_mixed_2D.cpp] along with the plot script that can be run by Matlab or Octave [https://gitlab.com/e62Lab/medusa poisson_mixed_2D.m]. Similar examples such as mixed boundary conditions on a circular domain, mixed boundary conditions on an irregular domain and a solution to a similar problem using monomial augmentation of radial basis functions can be found [https://gitlab.com/e62Lab/medusa here].

== Full Neumann boundary conditions ==
When solving the Navier-Stokes equation with Explicit Pressure correction, a Poisson's equation with full Neumann boundary conditions needs to be solved. In order for this type of problem to be solved, it needs to satisfy the compatibility condition. Inside the domain we still have:

$$ 
    \Delta u = f 
$$

The Poisson equation has a solution for the requisite function $u(x,y)$, when the integral of the source function $f(x,y)$ over our domain needs to be equal to the net flow rate, that is expressed by the boundary integral of the normal derivative of $u(x,y)$. In mathematical formulation, for our 2D case:

$$
    \int_{\Omega} f(x,y) dxdy = \int_{\partial \Omega} \frac{du}{d \b{\hat{n}}} dl
$$

Where again $\Omega$ is our domain, $\partial \Omega$ is the boundary of our domain, $dl$ is the differential arc length along the boundary. When the compatibility condition is met, the Poisson equation has a solution that can be determined up to an arbitrary constant. To obtain only one solution another constraint needs to be implemented. One option is to select a node and set it to a constant, e.g. p(0, 0) = 0, however much more stable approach is to enforce solution with an additional condition, also referred to as a regularization:

$$
    \int_{\Omega }^{{}}{u(x,y)d}\Omega = 0
$$

On our unit square domain, $\Omega = [0, 1] \times [0, 1]$ an example case that satisfies the compatibility condition would be:

<syntaxhighlight lang="c++" line>
    for (int i : domain.interior()) {
        op.lap(i) = 2;
    }

    int BOTTOM = -3;
    int TOP = -4;
    int LEFT = -1;
    int RIGHT = -2;

    for (int i : (domain.types() == LEFT)) {
        op.der1(i, 0) = 0.0;
    }
    for (int i : (domain.types() == RIGHT)) {
        op.der1(i, 0) = 1.0;
    }
    for (int i : (domain.types() == BOTTOM)) {
        op.der1(i, 1) = 0.0;
    }
    for (int i : (domain.types() == TOP)) {
        op.der1(i, 1) = 1.0;
    }
</syntaxhighlight>

Which is simply, $\Delta u = 2$ inside $\Omega$ while $\frac{du}{d \b{\hat{n}}} = 1$ on top and right boundary and $\frac{du}{d \b{\hat{n}}} = 0$ on the other two. With a little mental math we see that the compatibility condition is met, $2 * 1 = 0 + 1 + 0 + 1$. With that out of the way let us see how we implement regularization:

Well the idea add another line and column to matrix M holding the shape function s, and enforce the condition. We do this by adding a Lagrange multiplier $\alpha$:

$$
\left[ \begin{matrix}
   {{M}_{11}} & .. & {{M}_{1n}} & 1  \\
   .. & .. & .. & 1  \\
   {{M}_{n1}} & ... & {{M}_{nn}} & 1  \\
   1 & 1 & 1 & 0  \\
\end{matrix} \right]\left[ \begin{matrix}
   {{u}_{1}}  \\
   ...  \\
   {{u}_{n}}  \\
   \alpha   \\
\end{matrix} \right]=\left[ \begin{matrix}
   rhs_1  \\
   ...  \\
   rhs_n  \\
   0  \\
\end{matrix} \right]
$$

When the problem is well-posed $\alpha$ should be equal to zero, or at least very close. Let us fill our matrix M:

<syntaxhighlight lang="c++" line>
    // regularization
    // set the last row and column of the matrix
    for (int i = 0; i < N; ++i) {
        M.coeffRef(N, i) = 1;
        M.coeffRef(i, N) = 1;
    }
    // set the sum of all values
    rhs[N] = 0.0;
</syntaxhighlight>

If we take a look at the matrix we will see the extra row and column that enforce our constraint. They should be the last column and the last row of the matrix:

[[File:sparse_neumann_problem.png|800px]]

With this done we can solve the problem, just as we did in the examples above. 

[[File:poisson_neumann_2D.png|800px]]

For $dx = 0.05$ we get $\alpha = -1.03846 \cdot 10^{-13}$, which is expected from a valid solution. 

The whole example can be found as [https://gitlab.com/e62Lab/medusa poisson_neumann_2D.cpp] along with the plot script that can be run by Matlab or Octave [https://gitlab.com/e62Lab/medusa poisson_neumann_2D.m].

== 3D model case ==

A 3D domain <code>triceratops.h5</code> was created from the Mathematica's triceratops model. The point coordinates can be imported to medusa and used as a domain. In the example given in <code>triceratops.cpp</code> 
the equation $-\nabla^u = 1$ with $u=0$ on the boundary is solved.

The plot of the solution obtained using <code>triceratops.m</code> is shown below: