Burgers' equation - Revision history

VagifO: Added link to repository

2024-08-26T16:11:19Z

Added link to repository

VagifO: /* Solution in 1D */

2024-04-23T14:55:23Z

‎Solution in 1D

VagifO: Boundary conditions for 2D

2024-04-23T14:38:37Z

Boundary conditions for 2D

VagifO: /* Convergence */

2024-03-18T15:26:46Z

‎Convergence

VagifO at 15:25, 18 March 2024

2024-03-18T15:25:38Z

VagifO at 15:09, 15 March 2024

2024-03-15T15:09:18Z

VagifO: /* Convergence */

2024-03-07T15:16:50Z

‎Convergence

VagifO: /* Solution in 1D */

2024-03-07T15:15:01Z

‎Solution in 1D

VagifO at 15:11, 7 March 2024

2024-03-07T15:11:33Z

Show changes

VagifO at 17:15, 6 March 2024

2024-03-06T17:15:40Z

@@ Line 19: / Line 19: @@
 \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). First we construct our 1D domain with length 2l, discretize it with a given step dx and find the closest $n$ nodes of each node which are called stencil nodes or support nodes or simply the stencil.
+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 on our [https://gitlab.com/e62Lab/medusa/-/tree/dev/examples?ref_type=heads repository]). First we construct our 1D domain with length 2l, discretize it with a given step dx and find the closest $n$ nodes of each node which are called stencil nodes or support nodes or simply the stencil.
 <syntaxhighlight lang="c++" line>

@@ Line 58: / Line 58: @@
      //Set initial condition
      for (int i : discretization.all()){
-         u(i) = 2*visc*a*k*std::sin(k*discretization.pos(i,0))/(b+a*std::cos(k*discretization.pos(i,0)));
+         u(i) = 2*nu*a*k*std::sin(k*discretization.pos(i,0))/(b+a*std::cos(k*discretization.pos(i,0)));
 </syntaxhighlight>
@@ Line 75: / Line 75: @@
      // Setting the equation interior
      for (int i : interior) {
-         op.value(i) + (u[i] * dt) * op.der1(i, 0) + (-dt*visc) * op.der2(i, 0, 0) = u[i];
+         op.value(i) + (u[i] * dt) * op.der1(i, 0) + (-dt*nu) * op.der2(i, 0, 0) = u[i];
+}
 </syntaxhighlight>

@@ Line 151: / Line 151: @@
 GeneralFill<Vec2d> fill; fill.seed(0);
 discretization.fill(fill, step);
-// Find support nodes
+// Find support nodes, excluding boundary nodes from supports of boundary nodes for better stability of neumann BC
-discretization.findSupport(FindClosest(45));
+    discretization.findSupport(FindClosest(n).forNodes(discretization.interior()));
 int domain_size = discretization.size();
 </syntaxhighlight>
@@ Line 161: / Line 162: @@
 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.
+u(0,y,t)=u(1,y,t)=v(x,0,t)=v(x,1,t)=0,\\
 \end{gather*}

@@ Line 120: / Line 120: @@
 We fix the time step to $dt = 10^{-6}$, vary the monomial augmentation order $m \in \{2, 4, 6, 8\}$ and analyze the convergence with increasing number of discretization nodes $N$.
-After some initial oscillation the maximum error decreases approximately as $\left(\frac{1}{N}\right)^k$, which is proportional to $dx^k$, where $k$ is the fit parameter for the function: $\log_{10}e_\infty = k\log_{10}N + c$ for every $m$. But eventually the error plateaus because for finer spatial discretizations we are bounded by the time step.
+After some initial oscillation the maximum error decreases approximately as $\left(\frac{1}{N}\right)^k$, which is proportional to $dx^k$, where $k$ is the fit parameter for the function: $\log_{10}e_\infty = -k\log_{10}N + c$ for every $m$. But eventually the error plateaus because for finer spatial discretizations we are bounded by the time step.
 [[File:Error_step6.png|400px|center]]

@@ Line 19: / Line 19: @@
 \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). First we construct our 1D domain with length 2l, discretize it with a given step dx and find the closest n nodes of each node which are called stencil nodes or support nodes or simply the stencil.
+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). First we construct our 1D domain with length 2l, discretize it with a given step dx and find the closest $n$ nodes of each node which are called stencil nodes or support nodes or simply the stencil.
 <syntaxhighlight lang="c++" line>
@@ Line 42: / Line 42: @@
 </syntaxhighlight>
-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>
+Next we set the initial condition to
 \begin{equation}
-u(x, 0) = \frac{2\nu ake^{-\nu k^2t}\sin(kx)}{b + ae^{-\nu k^2t}\cos(kx)},
+u(x) = \frac{2\nu ak\sin(kx)}{b + a\cos(kx)},
 \end{equation}
@@ Line 65: / Line 70: @@
 </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.
+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 u[i] represents $\b{u}_i^{n-1}$ - result of the previous step.
 <syntaxhighlight lang="c++" line>
@@ Line 74: / Line 79: @@
 </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.
+Boundary conditions are handled in the same way. The analytical solution is defined on an infinite domain, but we see from equation (5) 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
@@ Line 106: / Line 111: @@
 =Results in 1D=
-For the simulation bellow we used the following parameters: $\nu = 0.05, l = 1, a = 4, b = 4.1, k = 2\pi$, n = 13, $m = 4$, $dx = 2\cdot 10^{-3}$, $dt = 10^{-6}$.
+For the simulation bellow we used the following parameters: $\nu = 0.05, l = 1, a = 4, b = 4.1, k = 2\pi$, $n = 13$, $m = 4$, $dx = 2\cdot 10^{-3}$, $dt = 10^{-6}$.
 [[File:Burgers_animation.mp4|400px|center]]
@@ Line 115: / Line 120: @@
 We fix the time step to $dt = 10^{-6}$, vary the monomial augmentation order $m \in \{2, 4, 6, 8\}$ and analyze the convergence with increasing number of discretization nodes $N$.
-After some initial oscilation the maximum error decreases approximately as $(\frac{1}{N})^k$, which is proportional to $dx^k$, where $k$ is the fit parameter for the function: $\log_{10}e_\infty = k\log_{10}N + c$ for every $m$. But eventualy the error plateaus because for finer spatial discretizations we are bounded by the time step.
+After some initial oscillation the maximum error decreases approximately as $\left(\frac{1}{N}\right)^k$, which is proportional to $dx^k$, where $k$ is the fit parameter for the function: $\log_{10}e_\infty = k\log_{10}N + c$ for every $m$. But eventually the error plateaus because for finer spatial discretizations we are bounded by the time step.
 [[File:Error_step6.png|400px|center]]
@@ Line 203: / Line 208: @@
 =Results in 2D=
-We use the RBF-FD method with PHS radial basis function $\phi(r) = r^3$ augmented with monomials. For the simulation bellow we used parameters $\nu = 1$, n = 45, $m = 4$, $dr = 0.02$, $dt = 10^{-5}$. The animations below show the velocities in both directions and the magnitude of the velocity.
+We use the RBF-FD method with PHS radial basis function $\phi(r) = r^3$ augmented with monomials. For the simulation bellow we used parameters $\nu = 1$, $n$ = 45, $m = 4$, $dr = 0.02$, $dt = 10^{-5}$. The animations below show the velocities in both directions and the magnitude of the velocity.
 [[File:Burgers2D_u.mp4|360px]] [[File:Burgers2D_v.mp4|360px]] [[File:Burgers2D_norm.mp4|360px]]
@@ Line 210: / Line 215: @@
 <b>Spatial discretization</b>
-We fix the time step to $dt = 10^{-6}$ and change the augmentation monomial order. The rate at which the maximum error decreases is close to $(\frac{1}{\sqrt{N}})^{k}$ as shown on the graph bellow. The distance between neighbouring nodes in 2D is approximately proportional to $\frac{1}{\sqrt{N}}$ if the node distribution is uniform. For finer discretizations and bigger $m$ the error becomes approxmately constant as we are again bounded by the time step.
+We fix the time step to $dt = 10^{-6}$ and change the augmentation monomial order. The rate at which the maximum error decreases is close to $\left(\frac{1}{\sqrt{N}}\right)^{k}$ as shown on the graph bellow. The distance between neighbouring nodes in 2D is approximately proportional to $\frac{1}{\sqrt{N}}$ if the node distribution is uniform. For finer discretizations and bigger $m$ the error becomes approximately constant as we are again bounded by the time step.
 [[File:Burgers2D_error_step.png|400px||center]]

@@ Line 113: / Line 113: @@
 <b>Spatial discretization</b>
-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.
+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.

@@ Line 41: / Line 41: @@
 </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>
+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}
@@ Line 65: / Line 65: @@
 </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.
+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>
@@ Line 102: / Line 102: @@
          ++t_save;
          hdf_out.writeDoubleArray("u", u);
-−
+    }
 </syntaxhighlight>

@@ Line 114: / Line 114: @@
 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]]
+[[File:Burgers_animation.mp4|400px|center]]
-=Convergence=
+==Convergence==
 <b>Spatial discretization</b>
@@ Line 123: / Line 123: @@
 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]]
+[[File:Error_step6.png|400px|center]]
 <b>Temporal discretization</b>
@@ Line 129: / Line 129: @@
 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]]
+[[File:Burgers1D_time_err.png|400px|center]]
 =Solution in 2D=
@@ Line 205: / Line 205: @@
 [[File:Burgers2D_u.mp4|360px]] [[File:Burgers2D_v.mp4|360px]] [[File:Burgers2D_norm.mp4|360px]]
 =References=