File:Burgers2D error step.png

2024-03-15T15:08:42Z

VagifO: VagifO uploaded a new version of File:Burgers2D error step.png

File:Burgers1D time err.png

2024-03-15T15:04:59Z

VagifO: VagifO uploaded a new version of File:Burgers1D time err.png

File:Error step6.png

2024-03-15T15:04:26Z

VagifO: VagifO uploaded a new version of File:Error step6.png

File:Error step6.png

2024-03-15T14:49:33Z

VagifO: VagifO uploaded a new version of File:Error step6.png

File:Burgers animation.mp4

2024-03-15T14:28:44Z

Burgers' equation

2024-03-07T15:15:01Z

VagifO: /* Solution in 1D */

Let us consider the Burgers' equation, which describes the dynamics of viscous fluid without the effects of pressure. Despite being unrealistic, it is the simplest description of advective flow with diffusive effects of viscosity. It has the following form:
\begin{equation}
\frac{\partial \b{u}}{\partial t} + (\b{u}\cdot\nabla)\b{u} = \nu \nabla^2 \b{u} ,
\end{equation}

where $\b{u}$ is a velocity field and $\nu$ is the dynamic viscosity. To solve the equation we will linearize the advection term. There are several tactics of linearization, but we choose the simplest one, which is using the previous velocity as an approximation of the current velociy. Using the implicit Euler method we arrive at the equation

\begin{equation}
\b{u}^{n}_i + \Bigl[\b{u}^{n-1}_i\cdot\nabla_i\b{u}^{n} - \nu\nabla^2_i \b{u}^{n}\Bigr]\Delta t = \b{u}^{n-1}_i,
\end{equation}

where the subscript $i$ is the index of the domain node, superscripts $n$ and $n-1$ denote the current and the previous time step respectively, $\Delta t$ is the length of the time step. The notation $\nabla_i$ means the gradient operator at node $i$, not the component of the gradient. Same holds for $\nabla^2_i$. Using the Medusa library, this equation can be solved for $\b{u}^{n}_i$.

=Solution in 1D=
Let's look at an example in 1D with Dirichlet boundary conditions. We are solving the equivalent of equation (2) but with only one spatial variable:
\begin{equation}
u^{n}_i + \Bigl[u^{n-1}_i\frac{\partial u^{n}}{\partial x} - \nu\frac{\partial^2 u^{n}}{\partial x^2}\Bigr]\Delta t = u^{n-1}_i ,\\
u_i = 0 \text{ on } \partial \Omega .
\end{equation}

This will be rewritten into a linear system of equations with Medusa and then solved with Eigen. Let's look at the code (find the full version in our repository).

<syntaxhighlight lang="c++" line>
// Build 1D domain
BoxShape<Vec1d> domain((-l), (l)); // the interval [-l, l]
auto discretization = domain.discretizeWithStep(dx);
int domain_size = discretization.size();
// Find support nodes
FindClosest find_support(5);
discretization.findSupport(find_support);
</syntaxhighlight>

First we construct our 1D domain with length 2l, discretize it with a given step dx and find the closest 5 nodes of each node which are called stencil nodes.

<syntaxhighlight lang="c++" line>
Polyharmonic<double, 3> p;
Monomials<Vec1d> mon(4);
RBFFD<Polyharmonic<double, 3>, Vec1d, ScaleToClosest> approx(p, mon);
auto storage = discretization.computeShapes(approx);
VectorXd u = VectorXd::Zero(domain_size); // vector for containing the current state
</syntaxhighlight>

We will use polyharmonic splines of order 3 as Radial Basis Functions (RBFs) augmented with monomials up to 4th order for our approximation engine. We use it in conjunction with stencil nodes to compute the stencil weights (also called shape functions) which are used in approximations of differential operators on our domain. Next we set the initial condition to $$u(x) = \frac{2\nu ak\sin(kx)}{b + a\cos(kx)},$$ where $a, b, k$ are constants and $\nu$ is viscosity. We chose this function because it has an analytical solution<ref name="Burgers' equation">A. Salih, Burgers' equation. Indian Institute of Space Science and Technology, Thiruvananthapuram, 2016.</ref>

\begin{equation}
u(x, 0) = \frac{2\nu ake^{-\nu k^2t}\sin(kx)}{b + ae^{-\nu k^2t}\cos(kx)},
\end{equation}

which we will compare our solution to. It is derived using the Cole-Hopf transformation of a solution of the diffusion equation with the initial profile $f(x) = b + a\cos(kx)$ for $b>a$.

<syntaxhighlight lang="c++" line>
//Set initial condition
for (int i : discretization.all()){
u(i) = 2*visc*a*k*std::exp(-visc*pow(k,2)*t_0)*std::sin(k*discretization.pos(i,0))/(b+a*std::exp(-
visc*pow(k,2)*t_0)*std::cos(k*discretization.pos(i,0)));
</syntaxhighlight>

Inside our time-stepping loop, we instantiate the matrix $M$ and the vector $rhs$ representing the linear system. We then construct the differential operators for implicit solving.
<syntaxhighlight lang="c++" line>
for (tt = 1; tt <= t_steps; ++tt) {
SparseMatrix<double> M(domain_size, domain_size);
VectorXd rhs = VectorXd::Zero(domain_size);
M.reserve(Range<int>(domain_size, n));
auto op = storage.implicitOperators(M, rhs);
</syntaxhighlight>

Now comes the crucial part - setting the equation. In essence, we rewrite equation (3) in Medusa's terms. Terms on the left hand side are written into matrix $M$ and the terms on the right into $rhs$. The terms with the prefix op. are matrices that represent the apropriate differential operators (op.value(i) is a matrix with $A_{ii}$ = 1 and zero everywhere else), vector u[i] represents $\b{u}_i^{n-1}$ - result of the previous step.

<syntaxhighlight lang="c++" line>
// Setting the equation
for (int i : interior) {
op.value(i) + (u[i] * dt) * op.der1(i, 0) + (-dt*visc) * op.der2(i, 0, 0) = u[i];
}
</syntaxhighlight>

Boundary conditions are handled in the same way. The analytical solution is defined on an infinite domain, but we see from equation (4) that $u = 0$ for $kx = z\pi, z \in \mathbb{Z}$. So we will be able to enforce Dirichlet boundary conditions on a finite domain, but this will limit our choice of $k$ to $\frac{z\pi}{l}, z \in \mathbb{Z}$, where $l$ is the length of the domain.
<syntaxhighlight lang="c++" line>
// Setting Dirichlet boundary conditions
for (int i : boundary) {
op.value(i) = 0;
}
</syntaxhighlight>

We then solve the equation and update u for use in the next time step. The solution is also saved into an output file with HDF5.

<syntaxhighlight lang="c++" line>
M.makeCompressed();
BiCGSTAB<SparseMatrix<double, RowMajor>> solver;
solver.compute(M);

// Solve the linear system and update the state
VectorXd u2 = solver.solve(rhs);
u = u2;

if (tt%(t_factor) == 0) {
// Saving state in output file
hdf_out.openGroup("/step" + std::to_string(t_save));
hdf_out.writeSparseMatrix("M", M);
hdf_out.writeDoubleAttribute("TimeStep", t_save);
std::cout<< "Time step " << t_save << " of " << t_steps/t_factor << "." << std::endl;
hdf_out.writeDoubleAttribute("time",(t_0 + tt * dt));
++t_save;
hdf_out.writeDoubleArray("u", u);
}
</syntaxhighlight>

=Results in 1D=
For our simulation we used the following parameters: $\nu = 1, l = 1, a = 4, b = 4, k = \frac{4\pi}{2l}$. As mentioned before the support of each node are its 5 closest nodes including the node itself. We used the RBF-FD method with polyharmonic splines of order 3 augmented with monomials of up to order $m = 4$. For stencils of minimal size - $m+1$, the method effectively becomes regular Finite Difference method (FD), which is true in our case.

[[File:Burgers_animation.mp4|400px|center]]

==Convergence==
Spatial discretization

We fix the time step to $dt = 10^{-6}$, vary the monomial augmentation order $m \in \{2, 4, 6, 8\}$, stencil size $n = m+1$ reducing our method to FD and analyze the convergence with increasing number of discretization nodes N.

After some initial oscilation the error decreases as $dr^k$, where $k$ is the fit parameter for every $m$ and $dr$ is the distance between consecutive nodes. But eventualy it becomes approximately constant. The reason is that for finer spatial discretizations we are bounded by the time step.

[[File:Error_step6.png|400px|center]]

Temporal discretization

Now we fix $N$ and $m$ to 1024 and 4 respectively and vary the size of the time step from $10^{-3}$ to $10^{-6}$. The error decreases linearly with the time step.

[[File:Burgers1D_time_err.png|400px|center]]

=Solution in 2D=
In 2D the Burgers' equation (1) gives a system of two coupled nonlinear equations:
\begin{gather}
\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \nu\nabla^2u,\\
\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = \nu\nabla^2v,
\end{gather}

where $u$ and $v$ are velocities in the $x$ and $y$ direction respectively. We discretize equations (4) and (5) in the same way as in the 1D case:
\begin{gather}
u^{n}_i + \Bigl[u^{n-1}_i\frac{\partial u^{n}}{\partial x} + v^{n-1}_i\frac{\partial u^{n}}{\partial y} - \nu\nabla^2u^{n}\Bigr]\Delta t = u^{n-1}_i ,\\
v^{n}_i + \Bigl[u^{n-1}_i\frac{\partial v^{n}}{\partial x} + v^{n-1}_i\frac{\partial v^{n}}{\partial y} - \nu\nabla^2v^{n}\Bigr]\Delta t = v^{n-1}_i .
\end{gather}

We will solve each equation separately every time step, this will enable us to partialy parallelize the algorithm. First we contruct a 2D square domain with side length of 1 and find 45 closest nodes of each node.

<syntaxhighlight lang="c++" line>
BoxShape<Vec2d> b((0, 0), (1, 1));
double step = 0.05;
auto discretization = b.discretizeBoundaryWithStep(step);
GeneralFill<Vec2d> fill; fill.seed(0);
discretization.fill(fill, step);
// Find support nodes
discretization.findSupport(FindClosest(45));
int domain_size = discretization.size();
</syntaxhighlight>

The approximation engine is created and shape functions are computed in the same way as in the 1D example, except the template parameter is <code>Vec2d</code> instead of <code>Vec1d</code>. We will be simulating the following initial and boundary conditions:

\begin{gather*}
u(x,y,0) = sin(\pi x)cos(\pi y),\\
v(x,y,0) = cos(\pi x)sin(\pi y),\\
u(0,y,t)=u(1,y,t)=v(x,0,t)=v(x,1,t)=0.
\end{gather*}

The analytical solution is given here <ref name="Analytical 2D">Q. Gao, M.Y. Zou, An analytical solution for two and three dimensional nonlinear Burgers' equation, Applied Mathematical Modelling, Volume 45, 2017, Pages 255-270, ISSN 0307-904X, https://doi.org/10.1016/j.apm.2016.12.018.</ref>. In the time stepping loop we use <code>#pragma omp parallel sections</code> to enable concurrent solving for $u$ and $v$. In each section we define separate matrices Mu or Mv, rhs vectors rhu and rhv, and implicit operators opu and opv which will write into their respective matrices. The code bellow includes only the section for $u$. The types -1 ,-2 , -3 ,-4 correspond to the left, right, bottom and top boundaries. At the end of each section the velocities are updated.
<syntaxhighlight lang="c++" line>
for (tt = 1; tt <= t_steps; ++tt) {

double t = tt * dt;

#pragma omp parallel sections shared(domain_size, discretization, storage, dt, u, v)
{
////Solving for u////
#pragma omp section
{
SparseMatrix<double, RowMajor> Mu(domain_size, domain_size);
VectorXd rhu = VectorXd::Zero(domain_size);
auto opu = storage.implicitOperators(Mu, rhu);
Mu.reserve(storage.supportSizes());

// Setting the equation
#pragma omp parallel for
for (int i : discretization.interior()) {
opu.value(i) + (u[i]*dt) * opu.der1(i, 0) + (v[i]*dt) * opu.der1(i, 1) + (-1/rndls*dt) * opu.lap(i) = u[i];
}
// Boundary conditions
for (int i : discretization.types() == -3){
opu.value(i) = analytical_u(discretization.pos(i,0), discretization.pos(i,1), t, 1.0/rndls);
}
for (int i : discretization.types() == -4){
opu.value(i) = analytical_u(discretization.pos(i,0), discretization.pos(i,1), t, 1.0/rndls);
}
for (int i : discretization.types() == -1){
opu.value(i) = 0.0;
}
for (int i : discretization.types() == -2){
opu.value(i) = 0.0;
}
Mu.makeCompressed();
BiCGSTAB<SparseMatrix<double, RowMajor>> solveru; // Creating solver object
solveru.compute(Mu); // Initialize the iterative solver
VectorXd u2 = solveru.solve(rhu);
u = u2;
}
</syntaxhighlight>

=Results in 2D=
We use the RBF-FD method with PHS radial basis function $\phi(r) = r^3$ augmented with monomials for our approximation engine. For the simulation we used $\nu = 1$. The animations below show the velocities in both directions as well as the magnitude of the velocity.

[[File:Burgers2D_u.mp4|360px]] [[File:Burgers2D_v.mp4|360px]] [[File:Burgers2D_norm.mp4|360px]]

==Convergence==
Spatial discretization

We fix the time step to $dt = 10^{-6}$ and change the augmentation polinomial order and analyze the convergence with increasing node density. The distance $dr$ between neighbouring nodes in 2D is approximately proportional to $\frac{1}{\sqrt{N}}$. The rate at which the error decreases is close to $dr^{-m}$ as show on the graph bellow. For finer discretizations and bigger $m$ the error becomes constant as we are again bounded by the time step.

[[File:Burgers2D_error_step.png|400px||center]]

Temporal discretization

We fix $N = 753$ and $m = 4$ and change the time step. The error decreases linearily as in the 1D case.

[[File:Burgers2D_time_error.png|400px||center]]
=References=

Burgers' equation

2024-03-07T15:11:33Z

VagifO:

Let us consider the Burgers' equation, which describes the dynamics of viscous fluid without the effects of pressure. Despite being unrealistic, it is the simplest description of advective flow with diffusive effects of viscosity. It has the following form:
\begin{equation}
\frac{\partial \b{u}}{\partial t} + (\b{u}\cdot\nabla)\b{u} = \nu \nabla^2 \b{u} ,
\end{equation}

where $\b{u}$ is a velocity field and $\nu$ is the dynamic viscosity. To solve the equation we will linearize the advection term. There are several tactics of linearization, but we choose the simplest one, which is using the previous velocity as an approximation of the current velociy. Using the implicit Euler method we arrive at the equation

\begin{equation}
\b{u}^{n}_i + \Bigl[\b{u}^{n-1}_i\cdot\nabla_i\b{u}^{n} - \nu\nabla^2_i \b{u}^{n}\Bigr]\Delta t = \b{u}^{n-1}_i,
\end{equation}

where the subscript $i$ is the index of the domain node, superscripts $n$ and $n-1$ denote the current and the previous time step respectively, $\Delta t$ is the length of the time step. The notation $\nabla_i$ means the gradient operator at node $i$, not the component of the gradient. Same holds for $\nabla^2_i$. Using the Medusa library, this equation can be solved for $\b{u}^{n}_i$.

=Solution in 1D=
Let's look at an example in 1D with Dirichlet boundary conditions. We are solving the equivalent of equation (2) but with only one spatial variable:
\begin{equation}
u^{n}_i + \Bigl[u^{n-1}_i\frac{\partial u^{n}}{\partial x} - \nu\frac{\partial^2 u^{n}}{\partial x^2}\Bigr]\Delta t = u^{n-1}_i ,\\
u_i = 0 \text{ on } \partial \Omega .
\end{equation}

This will be rewritten into a linear system of equations with Medusa and then solved with Eigen. Let's look at the code (find the full version in our repository).

<syntaxhighlight lang="c++" line>
// Build 1D domain
BoxShape<Vec1d> domain((-l), (l)); // the interval [-l, l]
auto discretization = domain.discretizeWithStep(dx);
int domain_size = discretization.size();
// Find support nodes
FindClosest find_support(5);
discretization.findSupport(find_support);
</syntaxhighlight>

First we construct our 1D domain with length 2l, discretize it with a given step dx and find the closest 5 nodes of each node which are called stencil nodes.

<syntaxhighlight lang="c++" line>
Polyharmonic<double, 3> p;
Monomials<Vec1d> mon(4);
RBFFD<Polyharmonic<double, 3>, Vec1d, ScaleToClosest> approx(p, mon);
auto storage = discretization.computeShapes(approx);
VectorXd u = VectorXd::Zero(domain_size); // vector for containing the current state
</syntaxhighlight>

We will use polyharmonic splines of order 3 as Radial Basis Functions(RBFs) augmented with monomials up to 4th order for our approximation engine. We use it in conjunction with stencil nodes to compute the stencil weights (also called shape functions) which are used in approximations of differential operators on our domain. Next we set the initial condition to $$u(x) = \frac{2\nu ak\sin(kx)}{b + a\cos(kx)},$$ where $a, b, k$ are constants and $\nu$ is viscosity. We chose this function because it has an analytical solution<ref name="Burgers' equation">A. Salih, Burgers' equation. Indian Institute of Space Science and Technology, Thiruvananthapuram, 2016.</ref>

\begin{equation}
u(x, 0) = \frac{2\nu ake^{-\nu k^2t}\sin(kx)}{b + ae^{-\nu k^2t}\cos(kx)},
\end{equation}

which we will compare our solution to. It is derived using the Cole-Hopf transformation of a solution of the diffusion equation with the initial profile $f(x) = b + a\cos(kx)$ for $b>a$.

<syntaxhighlight lang="c++" line>
//Set initial condition
for (int i : discretization.all()){
u(i) = 2*visc*a*k*std::exp(-visc*pow(k,2)*t_0)*std::sin(k*discretization.pos(i,0))/(b+a*std::exp(-
visc*pow(k,2)*t_0)*std::cos(k*discretization.pos(i,0)));
</syntaxhighlight>

Inside our time-stepping loop, we instantiate the matrix $M$ and the vector $rhs$ representing the linear system. We then construct the differential operators for implicit solving.
<syntaxhighlight lang="c++" line>
for (tt = 1; tt <= t_steps; ++tt) {
SparseMatrix<double> M(domain_size, domain_size);
VectorXd rhs = VectorXd::Zero(domain_size);
M.reserve(Range<int>(domain_size, n));
auto op = storage.implicitOperators(M, rhs);
</syntaxhighlight>

Now comes the crucial part - setting the equation. In essence, we rewrite equation (3) in Medusa's terms. Terms on the left hand side are written into matrix $M$ and the terms on the right into $rhs$. The terms with the prefix op. are matrices that represent the apropriate differential operators(op.value(i) is a matrix with $A_{ii}$ = 1 and zero everywhere else), vector E1[i] represents $\b{u}_i^{n-1}$ - result of the previous step.

<syntaxhighlight lang="c++" line>
// Setting the equation
for (int i : interior) {
op.value(i) + (u[i] * dt) * op.der1(i, 0) + (-dt*visc) * op.der2(i, 0, 0) = u[i];
}
</syntaxhighlight>

Boundary conditions are handled in the same way. The analytical solution is defined on an infinite domain, but we see from equation (4) that $u = 0$ for $kx = z\pi, z \in \mathbb{Z}$. So we will be able to enforce Dirichlet boundary conditions on a finite domain, but this will limit our choice of $k$ to $\frac{z\pi}{l}, z \in \mathbb{Z}$, where $l$ is the length of the domain.
<syntaxhighlight lang="c++" line>
// Setting Dirichlet boundary conditions
for (int i : boundary) {
op.value(i) = 0;
}
</syntaxhighlight>

We then solve the equation and update u for use in the next time step. The solution is also saved into an output file with HDF5.

<syntaxhighlight lang="c++" line>
M.makeCompressed();
BiCGSTAB<SparseMatrix<double, RowMajor>> solver;
solver.compute(M);

// Solve the linear system and update the state
VectorXd u2 = solver.solve(rhs);
u = u2;

if (tt%(t_factor) == 0) {
// Saving state in output file
hdf_out.openGroup("/step" + std::to_string(t_save));
hdf_out.writeSparseMatrix("M", M);
hdf_out.writeDoubleAttribute("TimeStep", t_save);
std::cout<< "Time step " << t_save << " of " << t_steps/t_factor << "." << std::endl;
hdf_out.writeDoubleAttribute("time",(t_0 + tt * dt));
++t_save;
hdf_out.writeDoubleArray("u", u);
}
</syntaxhighlight>

=Results in 1D=
For our simulation we used the following parameters: $\nu = 1, l = 1, a = 4, b = 4, k = \frac{4\pi}{2l}$. As mentioned before the support of each node are its 5 closest nodes including the node itself. We used the RBF-FD method with polyharmonic splines of order 3 augmented with monomials of up to order $m = 4$. For stencils of minimal size - $m+1$, the method effectively becomes regular Finite Difference method (FD), which is true in our case.

[[File:Burgers_animation.mp4|400px|center]]

==Convergence==
Spatial discretization

We fix the time step to $dt = 10^{-6}$, vary the monomial augmentation order $m \in \{2, 4, 6, 8\}$, stencil size $n = m+1$ reducing our method to FD and analyze the convergence with increasing number of discretization nodes N.

After some initial oscilation the error decreases as $dr^k$, where $k$ is the fit parameter for every $m$ and $dr$ is the distance between consecutive nodes. But eventualy it becomes approximately constant. The reason is that for finer spatial discretizations we are bounded by the time step.

[[File:Error_step6.png|400px|center]]

Temporal discretization

Now we fix $N$ and $m$ to 1024 and 4 respectively and vary the size of the time step from $10^{-3}$ to $10^{-6}$. The error decreases linearly with the time step.

[[File:Burgers1D_time_err.png|400px|center]]

=Solution in 2D=
In 2D the Burgers' equation (1) gives a system of two coupled nonlinear equations:
\begin{gather}
\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \nu\nabla^2u,\\
\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = \nu\nabla^2v,
\end{gather}

where $u$ and $v$ are velocities in the $x$ and $y$ direction respectively. We discretize equations (4) and (5) in the same way as in the 1D case:
\begin{gather}
u^{n}_i + \Bigl[u^{n-1}_i\frac{\partial u^{n}}{\partial x} + v^{n-1}_i\frac{\partial u^{n}}{\partial y} - \nu\nabla^2u^{n}\Bigr]\Delta t = u^{n-1}_i ,\\
v^{n}_i + \Bigl[u^{n-1}_i\frac{\partial v^{n}}{\partial x} + v^{n-1}_i\frac{\partial v^{n}}{\partial y} - \nu\nabla^2v^{n}\Bigr]\Delta t = v^{n-1}_i .
\end{gather}

We will solve each equation separately every time step, this will enable us to partialy parallelize the algorithm. First we contruct a 2D square domain with side length of 1 and find 45 closest nodes of each node.

<syntaxhighlight lang="c++" line>
BoxShape<Vec2d> b((0, 0), (1, 1));
double step = 0.05;
auto discretization = b.discretizeBoundaryWithStep(step);
GeneralFill<Vec2d> fill; fill.seed(0);
discretization.fill(fill, step);
// Find support nodes
discretization.findSupport(FindClosest(45));
int domain_size = discretization.size();
</syntaxhighlight>

The approximation engine is created and shape functions are computed in the same way as in the 1D example, except the template parameter is <code>Vec2d</code> instead of <code>Vec1d</code>. We will be simulating the following initial and boundary conditions:

\begin{gather*}
u(x,y,0) = sin(\pi x)cos(\pi y),\\
v(x,y,0) = cos(\pi x)sin(\pi y),\\
u(0,y,t)=u(1,y,t)=v(x,0,t)=v(x,1,t)=0.
\end{gather*}

The analytical solution is given here <ref name="Analytical 2D">Q. Gao, M.Y. Zou, An analytical solution for two and three dimensional nonlinear Burgers' equation, Applied Mathematical Modelling, Volume 45, 2017, Pages 255-270, ISSN 0307-904X, https://doi.org/10.1016/j.apm.2016.12.018.</ref>. In the time stepping loop we use <code>#pragma omp parallel sections</code> to enable concurrent solving for $u$ and $v$. In each section we define separate matrices Mu or Mv, rhs vectors rhu and rhv, and implicit operators opu and opv which will write into their respective matrices. The code bellow includes only the section for $u$. The types -1 ,-2 , -3 ,-4 correspond to the left, right, bottom and top boundaries. At the end of each section the velocities are updated.
<syntaxhighlight lang="c++" line>
for (tt = 1; tt <= t_steps; ++tt) {

double t = tt * dt;

#pragma omp parallel sections shared(domain_size, discretization, storage, dt, u, v)
{
////Solving for u////
#pragma omp section
{
SparseMatrix<double, RowMajor> Mu(domain_size, domain_size);
VectorXd rhu = VectorXd::Zero(domain_size);
auto opu = storage.implicitOperators(Mu, rhu);
Mu.reserve(storage.supportSizes());

// Setting the equation
#pragma omp parallel for
for (int i : discretization.interior()) {
opu.value(i) + (u[i]*dt) * opu.der1(i, 0) + (v[i]*dt) * opu.der1(i, 1) + (-1/rndls*dt) * opu.lap(i) = u[i];
}
// Boundary conditions
for (int i : discretization.types() == -3){
opu.value(i) = analytical_u(discretization.pos(i,0), discretization.pos(i,1), t, 1.0/rndls);
}
for (int i : discretization.types() == -4){
opu.value(i) = analytical_u(discretization.pos(i,0), discretization.pos(i,1), t, 1.0/rndls);
}
for (int i : discretization.types() == -1){
opu.value(i) = 0.0;
}
for (int i : discretization.types() == -2){
opu.value(i) = 0.0;
}
Mu.makeCompressed();
BiCGSTAB<SparseMatrix<double, RowMajor>> solveru; // Creating solver object
solveru.compute(Mu); // Initialize the iterative solver
VectorXd u2 = solveru.solve(rhu);
u = u2;
}
</syntaxhighlight>

=Results in 2D=
We use the RBF-FD method with PHS radial basis function $\phi(r) = r^3$ augmented with monomials for our approximation engine. For the simulation we used $\nu = 1$. The animations below show the velocities in both directions as well as the magnitude of the velocity.

[[File:Burgers2D_u.mp4|360px]] [[File:Burgers2D_v.mp4|360px]] [[File:Burgers2D_norm.mp4|360px]]

==Convergence==
Spatial discretization

We fix the time step to $dt = 10^{-6}$ and change the augmentation polinomial order and analyze the convergence with increasing node density. The distance $dr$ between neighbouring nodes in 2D is approximately proportional to $\frac{1}{\sqrt{N}}$. The rate at which the error decreases is close to $dr^{-m}$ as show on the graph bellow. For finer discretizations and bigger $m$ the error becomes constant as we are again bounded by the time step.

[[File:Burgers2D_error_step.png|400px||center]]

Temporal discretization

We fix $N = 753$ and $m = 4$ and change the time step. The error decreases linearily as in the 1D case.

[[File:Burgers2D_time_error.png|400px||center]]
=References=

Burgers' equation

2024-03-06T17:15:40Z

2024-03-05T17:06:20Z

VagifO:

Let us consider the Burgers' equation, which describes the dynamics of viscous fluid without the effects of pressure. Despite being unrealistic, it is the simplest description of advective flow with diffusive effects of viscosity. It has the following form:
\begin{equation}
\frac{\partial \b{u}}{\partial t} + (\b{u}\cdot\nabla)\b{u} = \nu \nabla^2 \b{u} ,
\end{equation}

where $\b{u}$ is a velocity field. For time stepping the advection term must be linearized. There are several tactics of linearization, but we choose the simplest one, which is using the previous velocity as an approximantion of the current velociy. Using the implicit Euler method we arrive at the equation

\begin{equation}
\b{u}^{n}_i + \Bigl[\b{u}^{n-1}_i\cdot\nabla_i\b{u}^{n} - \nu\nabla^2_i \b{u}^{n}\Bigr]\Delta t = \b{u}^{n-1}_i,
\end{equation}

where the subscript $i$ is the index of the domain node, superscripts $n$ and $n-1$ denote the current and the previous time step respectively, $\Delta t$ is the length of the time step. The notation $\nabla_i$ means the gradient operator at node $i$, not the component of the gradient. Same holds for $\nabla^2_i$. Using the Medusa library, this equation can be solved for $\b{u}^{n}_i$.

=Solution in 1D=
Let's look at an example in 1D with Dirichlet boundary conditions. We are solving the equivalent of equation (2) but with only one spatial variable:
\begin{equation}
u^{n}_i + \Bigl[u^{n-1}_i\frac{\partial u^{n}}{\partial x} - \nu\frac{\partial^2 u^{n}}{\partial x^2}\Bigr]\Delta t = u^{n-1}_i ,\\
u_i = 0 \text{ on } \partial \Omega .
\end{equation}

This will be rewritten into a linear system of equations with Medusa and then solved with Eigen. Let's look at the code(find the full version in our repository).

<syntaxhighlight lang="c++" line>
// Build 1D domain
BoxShape<Vec1d> domain((-l), (l)); // the interval [-l, l]
auto discretization = domain.discretizeWithStep(dx);
int domain_size = discretization.size();
// Find support nodes
FindClosest find_support(5);
discretization.findSupport(find_support);
</syntaxhighlight>

First we construct our 1D domain with lenght 2l, discretize it with a given step dx and find the closest 5 nodes of each node which are called stencil nodes.

<syntaxhighlight lang="c++" line>
Polyharmonic<double, 3> p;
Monomials<Vec1d> mon(4);
RBFFD<Polyharmonic<double, 3>, Vec1d, ScaleToClosest>
approx(p, mon);
auto storage = discretization.computeShapes(approx);
VectorXd E1 = VectorXd::Zero(domain_size); // vector for containing the current state
</syntaxhighlight>

We will use polyharmonic splines of order 3 as Radial Basis Functions(RBFs) augmented with monomials up to 4th order for our approximation engine. We use it in conjunction with stencil nodes to compute the stencil weights(also called shape functions) which are used in approximations of differential operators on our domain. Next we set the initial condition to $$u(x, 0) = \frac{2\nu ake^{-\nu k^2t}\sin(kx)}{b + ae^{-\nu k^2t}\cos(kx)},$$ where $a, b, k$ are constants. We chose this function because it has an analytical solution which we will compare our solution to. It is derived using the Cole-Hopf transformation of a solution of the diffusion equation with the initial profile $f(x) = b + a\cos(kx)$ for $b>a$.

<syntaxhighlight lang="c++" line>
//Set initial condition
for (int i : discretization.all()){
E1(i) = 2*visc*a*k*std::exp(-visc*pow(k,2)*t_0)*std::sin(k*discretization.pos(i,0))/(b+a*std::exp(-
visc*pow(k,2)*t_0)*std::cos(k*discretization.pos(i,0)));
</syntaxhighlight>

Inside our time-stepping loop, we instantiate the matrix $M$ and the vector $rhs$ representing the linear system. We then construct the differential operators for implicit solving.
<syntaxhighlight lang="c++" line>
for (tt = 1; tt <= t_steps; ++tt) {
SparseMatrix<double> M(domain_size, domain_size); // construct matrix of apropreate size
VectorXd rhs = VectorXd::Zero(domain_size);
M.reserve(Range<int>(domain_size, n));
auto op = storage.implicitOperators(M, rhs); // compute the implicit operators
</syntaxhighlight>

Now comes the crucial part - setting the equation. In essence, we rewrite equation (3) in Medusa's terms. Terms on the left hand side are written into matrix $M$ and the terms on the right into $rhs$. The terms with the prefix op. are matrices that represent the apropriate differential operators(op.value(i) is a matrix with $A_{ii}$ = 1 and zero everywhere else), vector E1[i] represents $\b{u}_i^{n-1}$ - result of the previous step.

<syntaxhighlight lang="c++" line>
// Setting the equation
for (int i : interior) {
op.value(i) + (E1[i] * dt) * op.der1(i, 0) + (-dt*visc) * op.der2(i, 0, 0) = E1[i];
}
</syntaxhighlight>

Boundary conditions are handled in the same way. Here we set Dirichlet boundary conditions.
<syntaxhighlight lang="c++" line>
// Setting Dirichlet boundary conditions
for (int i : boundary) {
op.value(i) = 0;
}
</syntaxhighlight>

We then solve the equation and update E1 for use in the next time step. The solution is also saved into an output file with HDF5.

<syntaxhighlight lang="c++" line>
M.makeCompressed(); // Optimization
BiCGSTAB<SparseMatrix<double, RowMajor>> solver; // Creating solver object
solver.compute(M); // Initialize the iterative solver

// Solve the linear system and update the state
VectorXd E2 = solver.solve(rhs);
E1 = E2;

if (tt%(t_factor) == 0) {
// Saving state in output file
hdf_out.openGroup("/step" + std::to_string(t_save));
hdf_out.writeSparseMatrix("M", M);
hdf_out.writeDoubleAttribute("TimeStep", t_save);
std::cout<< "Time step " << t_save << " of " << t_steps/t_factor << "." << std::endl;
hdf_out.writeDoubleAttribute("time",(t_0 + tt * dt));
++t_save;
hdf_out.writeDoubleArray("E", E1);
}
</syntaxhighlight>

=Results in 1D=
For our simulation we used the following parameters:
\begin{align*}
\nu &= 1,\\
l &= 1,\\
a &= 4,\\
b &= 4,\\
k &= \frac{4\pi}{2l},\\
\phi(r) &= r^3,\\
n &= 5,\\
m &= 4.
\end{align*}
As mentioned before the support of each node are its 5 closest nodes including the node itself. We used the RBF-FD method with polyharmonic splines of order 3 augmented with polynomials of up to order 4. Animation of the result is shown below.

[[File:Burgers_animation.mp4|400px]]

=Convergence=
Spatial discretization

We use the RBF-FD method with PHS radial basis function $\phi(r) = r^3$ and fix the time step to $dt = 10^{-6}$. We vary the polynomial augmentation order $m \in \{2, 4, 6, 8\}$ and analize the convergence with increasing number of discretization nodes N.

After some initial oscilation the error decreases as $O(h^k)$, where $k$ is the fit parameter for every $m$ and $h$ is the distance between consecutive nodes. But eventualy it becomes constant value. The reason is that for finer spatial discretizations we are bounded by the time step because of linearization.

[[File:Error_step6.png|400px]]

Temporal discretization

Now we fix $N$ and $m$ to 1024 and 4 and vary the size of the time step from $10^{-3}$ to $10^{-6}$. The error decreases linearly with the time step.

[[File:Burgers_time_err.png|400px]]

=Solution in 2D=
In 2D the Burgers' equation (1) gives a system of two coupled nonlinear equations:
\begin{gather}
\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = \nu\nabla^2u,\\
\frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = \nu\nabla^2v,
\end{gather}
where $u$ and $v$ are ve;ocities in the $x$ and $y$ direction respectively. We discretize equations (4) and (5) in the same way as in the 1D case:
\begin{gather}
u^{n}_i + \Bigl[u^{n-1}_i\frac{\partial u^{n}}{\partial x} + v^{n-1}_i\frac{\partial u^{n}}{\partial y} - \nu\nabla^2u^{n}\Bigr]\Delta t = u^{n-1}_i ,\\
v^{n}_i + \Bigl[u^{n-1}_i\frac{\partial v^{n}}{\partial x} + v^{n-1}_i\frac{\partial v^{n}}{\partial y} - \nu\nabla^2v^{n}\Bigr]\Delta t = v^{n-1}_i .
\end{gather}
We will solve each equation separately every time step, this will enable us to partialy parallelize the algorithm. First we contruct a 2D square domain with side length of 1 and find n closest nodes of each node.
<syntaxhighlight lang="c++" line>
BoxShape<Vec2d> b((0, 0), (a, a));
double step = std::stod(argv[1]);
auto discretization = b.discretizeBoundaryWithStep(step);
GeneralFill<Vec2d> fill; fill.seed(0);
discretization.fill(fill, step);
// Find support nodes
discretization.findSupport(FindClosest(n));
int domain_size = discretization.size();
</syntaxhighlight>
We will be simulating the following initial and boundary conditions:
\begin{gather*}
u(x,y,0) = sin(\pi x)cos(\pi y),\\
v(x,y,0) = cos(\pi x)sin(\pi y),\\
u(0,y,t)=u(1,y,t)=v(x,0,t)=v(x,1,t)=0.
\end{gather*}
The analytical solution is given here <ref name="Analytical 2D">Q. Gao, M.Y. Zou, An analytical solution for two and three dimensional nonlinear Burgers' equation, Applied Mathematical Modelling, Volume 45, 2017, Pages 255-270, ISSN 0307-904X, https://doi.org/10.1016/j.apm.2016.12.018.</ref>. In the time stepping loop we use <code>#pragma omp parallel sections</code> to enable concurrent solving for $u$ and $v$. In each section we define separate matricies Mu or Mv, rhs vectors rhu and rhv, and implicit operators opu and opv which will write into their respective matricies. The code bellow includes only the section for $u$. The types -1 ,-2 , -3 ,-4 correspond to the left, right, bottom and top boundaries. At the end of each section the velocities are updated.
<syntaxhighlight lang="c++" line>
for (tt = 1; tt <= t_steps; ++tt) {

double t = t_0 + tt * dt;

#pragma omp parallel sections shared(domain_size, discretization, storage, dt, u, v)
{
////Solving for u////
#pragma omp section
{
SparseMatrix<double, RowMajor> Mu(domain_size, domain_size); // construct matrix of apropreate size
VectorXd rhu = VectorXd::Zero(domain_size); // set empty vector for rhs
auto opu = storage.implicitOperators(Mu, rhu);
Mu.reserve(storage.supportSizes());

// Setting the equation
#pragma omp parallel for
for (int i : discretization.interior()) {
opu.value(i) + (u[i]*dt) * opu.der1(i, 0) + (v[i]*dt) * opu.der1(i, 1) + (-1/rndls*dt) * opu.lap(i) = u[i];
}
// Boundary conditions
for (int i : discretization.types() == -3){
opu.value(i) = analytical_u(discretization.pos(i,0), discretization.pos(i,1), t, 1.0/rndls);
}
for (int i : discretization.types() == -4){
opu.value(i) = analytical_u(discretization.pos(i,0), discretization.pos(i,1), t, 1.0/rndls);
}
for (int i : discretization.types() == -1){
opu.value(i) = 0.0;
}
for (int i : discretization.types() == -2){
opu.value(i) = 0.0;
}
Mu.makeCompressed();
BiCGSTAB<SparseMatrix<double, RowMajor>> solveru; // Creating solver object
solveru.compute(Mu); // Initialize the iterative solver
VectorXd u2 = solveru.solve(rhu);
u = u2;
}
</syntaxhighlight>

=Results in 2D=
For all of the simulation we used $\nu = 1$. The animations below show the velocities in both directions as well as the magnitude of the velocity.
[[File:Burgers2D_u.mp4|360px]] [[File:Burgers2D_v.mp4|360px]] [[File:Burgers2D_norm.mp4|360px]]

=References=

File:Burgers2D v.mp4

2024-03-05T17:03:37Z

VagifO:

File:Burgers2D u.mp4

2024-03-05T17:03:22Z

VagifO:

2024-01-23T11:12:28Z

VagifO:

Let us consider the Burger's equation, which describes the dynamics of viscous fluid without the effects of pressure. Despite being unrealistic, it is the simplest description of advective flow with diffusive effects of viscosity. It has the following form:
\begin{equation}
\frac{\partial \b{u}}{\partial t} + (\b{u}\cdot\nabla)\b{u} = \nu \nabla^2 \b{u} ,
\end{equation}

where $\b{u}$ is a velocity field. For time stepping the advection term must be linearized. There are several tactics of linearization, but we choose the simplest one, which is using the previous velocity as an approximant of the current velociy. Using the implicit Euler method we arrive at the equation

\begin{equation}
\b{u}^{n}_i + \Bigl[\b{u}^{n-1}_i\cdot\nabla_i\b{u}^{n} - \nu\nabla^2_i \b{u}^{n}\Bigr]\Delta t = \b{u}^{n-1}_i,
\end{equation}

where the subscript $i$ is the index of the domain node, superscripts $n$ and $n-1$ denote the current and the previous time step respectively, $\Delta t$ is the length of the time step. The notation $\nabla_i$ means the gradient operator at node $i$, not the component of the gradient. Same holds for $\nabla^2_i$. Using the Medusa library, this equation can be solved for $\b{u}^{n}_i$.

=Solution in 1D=
Let's look at an example in 1D with Dirichlet boundary conditions. We are solving the equivalent of equation (2) but with only one spatial variable:
\begin{equation}
u^{n}_i + \Bigl[u^{n-1}_i\frac{\partial u^{n}}{\partial x} - \nu\frac{\partial^2 u^{n}}{\partial x^2}\Bigr]\Delta t = u^{n-1}_i ,\\
u_i = 0 \text{ on } \partial \Omega .
\end{equation}

This will be rewritten into a linear system of equations with Medusa and then solved with Eigen. Let's look at the code(find the full version in our repository)

<syntaxhighlight lang="c++" line>
// Build 1D domain
BoxShape<Vec1d> domain((-l), (l)); // the interval [-l, l]
auto discretization = domain.discretizeWithStep(dx);
int domain_size = discretization.size();
// Find support nodes
FindClosest find_support(5);
discretization.findSupport(find_support);
</syntaxhighlight>

First we construct our 1D domain with lenght 2l, discretize it with a given step dx and find the closest 5 nodes of each node which are called stencil nodes.

<syntaxhighlight lang="c++" line>
Polyharmonic<double, 3> p;
Monomials<Vec1d> mon(4);
RBFFD<Polyharmonic<double, 3>, Vec1d, ScaleToClosest>
approx(p, mon);
auto storage = discretization.computeShapes(approx);
VectorXd E1 = VectorXd::Zero(domain_size); // vector for containing the current state
</syntaxhighlight>

We will use polyharmonic splines of order 3 as Radial Basis Functions(RBFs) augmented with monomials up to 4th order for our approximation engine. We use it in conjunction with stencil nodes to compute the stencil weights(also called shape functions) which are used in approximations of differential operators on our domain. Next we set the initial condition to $$u(x, 0) = \frac{2\nu ake^{-\nu k^2t}\sin(kx)}{b + ae^{-\nu k^2t}\cos(kx)},$$ where $a, b, k$ are constants. We chose this function because it has an analytical solution which we will compare our solution to. It is derived using the Cole-Hopf transformation of a solution of the diffusion equation with the initial profile $f(x) = b + a\cos(kx)$ for $b>a$.

<syntaxhighlight lang="c++" line>
//Set initial condition
for (int i : discretization.all()){
E1(i) = 2*visc*a*k*std::exp(-visc*pow(k,2)*t_0)*std::sin(k*discretization.pos(i,0))/(b+a*std::exp(-
visc*pow(k,2)*t_0)*std::cos(k*discretization.pos(i,0)));
</syntaxhighlight>

Inside our time-stepping loop, we instantiate the matrix $M$ and the vector $rhs$ representing the linear system. We then construct the differential operators for implicit solving.
<syntaxhighlight lang="c++" line>
for (tt = 1; tt <= t_steps; ++tt) {
SparseMatrix<double> M(domain_size, domain_size); // construct matrix of apropreate size
VectorXd rhs = VectorXd::Zero(domain_size);
M.reserve(Range<int>(domain_size, n));
auto op = storage.implicitOperators(M, rhs); // compute the implicit operators
</syntaxhighlight>

Now comes the crucial part - setting the equation. In essence, we rewrite equation (3) in Medusa's terms. Terms on the left hand side are written into matrix $M$ and the terms on the right into $rhs$. The terms with the prefix op. are matrices that represent the appropriate differential operators(op.value(i) is a matrix with $A_{ii}$ = 1 and zero everywhere else), vector E1[i] represents $\b{u}_i^{n-1}$ - result of the previous step.

<syntaxhighlight lang="c++" line>
// Setting the equation
for (int i : interior) {
op.value(i) + (E1[i] * dt) * op.der1(i, 0) + (-dt*visc) * op.der2(i, 0, 0) = E1[i];
}
</syntaxhighlight>

Boundary conditions are handled in the same way. Here we set Dirichlet boundarz conditions.
<syntaxhighlight lang="c++" line>
// Setting Dirichlet boundary conditions
for (int i : boundary) {
op.value(i) = 0;
}
</syntaxhighlight>

We then solve the equation and update E1 for use in the next time step. The solution is also saved into an output file with HDF5.

<syntaxhighlight lang="c++" line>
M.makeCompressed(); // Optimization
BiCGSTAB<SparseMatrix<double, RowMajor>> solver; // Creating solver object
solver.compute(M); // Initialize the iterative solver

// Solve the linear system and update the state
VectorXd E2 = solver.solve(rhs);
E1 = E2;

if (tt%(t_factor) == 0) {
// Saving state in output file
hdf_out.openGroup("/step" + std::to_string(t_save));
hdf_out.writeSparseMatrix("M", M);
hdf_out.writeDoubleAttribute("TimeStep", t_save);
std::cout<< "Time step " << t_save << " of " << t_steps/t_factor << "." << std::endl;
hdf_out.writeDoubleAttribute("time",(t_0 + tt * dt));
++t_save;
hdf_out.writeDoubleArray("E", E1);
}
</syntaxhighlight>

=Results=
For our simulation we used the following parameters:
\begin{align*}
\nu &= 1\\
l &= 1\\
a &= 4\\
b &= 4\\
k &= \frac{4\pi}{2l}\\
\phi(r) &= r^3\\
n &= 5\\
m &= 4
\end{align*}
As mentioned before the support of each node are its 5 closest nodes including the node itself. We used the RBF-FD method with polyharmonic splines of order 3 augmented with polynomials of up to order 4. Animation of the result is shown below.

[[File:Burgers_animation.mp4|400px]]

=Convergence=
Spatial discretization

We use the RBF-FD method with PHS radial basis function $\phi(r) = r^3$ and fix the time step to $dt = 10^{-6}$. We vary the polynomial augmentation order $m \in {2, 4, 6, 8}$ and analize the convergence with increasing number of discretization nodes N.

After some initial oscilation the error decreases as $O(h^k)$, where $k$ is the fit parameter for every $m$. But eventualy it reaches a constant value. We can see that for bigger $N$ we are bounded by the time step because of linearization.

[[File:Error_step6.png|400px]]

Temporal discretization

Now we fix $N = 1024$ and $m = 4$ and vary the size of the time step from $10^{-3}$ to $10^{-6}$. The error decreases linearly with the time step.

[[File:Burgers_time_err.png|400px]]