Burgers' equation - Revision history

VagifO: Added link to repository

2024-08-26T16:11:19Z

Added link to repository

← Older revision		Revision as of 16:11, 26 August 2024
Line 19:		Line 19:
	\end{equation}		\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>		<syntaxhighlight lang="c++" line>

VagifO: /* Solution in 1D */

2024-04-23T14:55:23Z

Solution in 1D

← Older revision		Revision as of 14:55, 23 April 2024
Line 58:		Line 58:
	//Set initial condition		//Set initial condition
	for (int i : discretization.all()){		for (int i : discretization.all()){
	u(i) = 2~~visc~~akstd::sin(kdiscretization.pos(i,0))/(b+astd::cos(k*discretization.pos(i,0)));		u(i) = 2nuakstd::sin(kdiscretization.pos(i,0))/(b+astd::cos(k*discretization.pos(i,0)));
	</syntaxhighlight>		</syntaxhighlight>

Line 75:		Line 75:
	// Setting the equation interior		// Setting the equation interior
	for (int i : 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) + (-dtnu) op.der2(i, 0, 0) = u[i];
	}		}
	</syntaxhighlight>		</syntaxhighlight>

VagifO: Boundary conditions for 2D

2024-04-23T14:38:37Z

Boundary conditions for 2D

← Older revision		Revision as of 14:38, 23 April 2024
Line 151:		Line 151:
	GeneralFill<Vec2d> fill; fill.seed(0);		GeneralFill<Vec2d> fill; fill.seed(0);
	discretization.fill(fill, step);		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()));
			discretization.findSupport(FindClosest(n).forNodes(discretization.boundary()).searchAmong(discretization.interior()).forceSelf(true));
	int domain_size = discretization.size();		int domain_size = discretization.size();
	</syntaxhighlight>		</syntaxhighlight>
Line 161:		Line 162:
	u(x,y,0) = sin(\pi x)cos(\pi y),\\		u(x,y,0) = sin(\pi x)cos(\pi y),\\
	v(x,y,0) = cos(\pi x)sin(\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,\\
			\frac{\partial u}{\partial y}\bigg\|_{y = 0} = \frac{\partial u}{\partial y}\bigg\|_{y = 1} = \frac{\partial v}{\partial x}\bigg\|_{x = 0} = \frac{\partial v}{\partial x}\bigg\|_{x = 1} = 0.
	\end{gather*}		\end{gather*}

VagifO: /* Convergence */

2024-03-18T15:26:46Z

Convergence

← Older revision		Revision as of 15:26, 18 March 2024
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$.		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]]		[[File:Error_step6.png\|400px\|center]]

VagifO at 15:25, 18 March 2024

2024-03-18T15:25:38Z

← Older revision		Revision as of 15:25, 18 March 2024
Line 19:		Line 19:
	\end{equation}		\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>		<syntaxhighlight lang="c++" line>
Line 42:		Line 42:
	</syntaxhighlight>		</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}		\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}
			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, t) = \frac{2\nu ake^{-\nu k^2t}\sin(kx)}{b + ae^{-\nu k^2t}\cos(kx)},
	\end{equation}		\end{equation}

Line 65:		Line 70:
	</syntaxhighlight>		</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>		<syntaxhighlight lang="c++" line>
Line 74:		Line 79:
	</syntaxhighlight>		</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>		<syntaxhighlight lang="c++" line>
	// Setting Dirichlet boundary conditions		// Setting Dirichlet boundary conditions
Line 106:		Line 111:

	=Results in 1D=		=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]]		[[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$.		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]]		[[File:Error_step6.png\|400px\|center]]
Line 203:		Line 208:

	=Results in 2D=		=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]]		[[File:Burgers2D_u.mp4\|360px]] [[File:Burgers2D_v.mp4\|360px]] [[File:Burgers2D_norm.mp4\|360px]]
Line 210:		Line 215:
	<b>Spatial discretization</b>		<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]]		[[File:Burgers2D_error_step.png\|400px\|\|center]]

VagifO at 15:09, 15 March 2024

2024-03-15T15:09:18Z

Show changes

VagifO: /* Convergence */

2024-03-07T15:16:50Z

Convergence

← Older revision		Revision as of 15:16, 7 March 2024
Line 113:		Line 113:
	<b>Spatial discretization</b>		<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.		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.

VagifO: /* Solution in 1D */

2024-03-07T15:15:01Z

Solution in 1D

← Older revision		Revision as of 15:15, 7 March 2024
Line 41:		Line 41:
	</syntaxhighlight>		</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}		\begin{equation}
Line 65:		Line 65:
	</syntaxhighlight>		</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>		<syntaxhighlight lang="c++" line>
Line 102:		Line 102:
	++t_save;		++t_save;
	hdf_out.writeDoubleArray("u", u);		hdf_out.writeDoubleArray("u", u);
	}		}
	</syntaxhighlight>		</syntaxhighlight>

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

← Older revision		Revision as of 17:15, 6 March 2024
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.		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>		<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.		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>		<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.		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=		=Solution in 2D=
Line 205:		Line 205:
	[[File:Burgers2D_u.mp4\|360px]] [[File:Burgers2D_v.mp4\|360px]] [[File:Burgers2D_norm.mp4\|360px]]		[[File:Burgers2D_u.mp4\|360px]] [[File:Burgers2D_v.mp4\|360px]] [[File:Burgers2D_norm.mp4\|360px]]

			==Convergence==
			<b>Spatial discretization</b>

			We use the same method as in the 1D case. We fix the time step to $dt = 10^{-6}$ and change the augmentation polinomial order and analize 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 the expected $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]]

			<b>Temporal discretization</b>

			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=		=References=