Revision as of 13:49, 16 January 2024

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 current velocity as an approximant of the velociy in the next time step. Using the implicit Euler method we arrive at the equation

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

where the subscript $i$ is the index of the domain node, superscripts $n$ and $n+1$ denote the current and the next 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+1}_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+1}_i + [u^{n}_i\frac{\partial u^{n+1}}{\partial x} - \nu\frac{\partial^2 u^{n+1}}{\partial x^2}]\Delta t = u^{n}_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)

1 // Build 1D domain
2 BoxShape<Vec1d> domain((-l), (l));
3 auto discretization = domain.discretizeWithStep(dx);
4 int domain_size = discretization.size();
5 // Find support nodes
6 FindClosest find_support(5);
7 discretization.findSupport(find_support);

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.

1 Polyharmonic<double, 3> p;
2 Monomials<Vec1d> mon(4);
3 RBFFD<Polyharmonic<double, 3>, Vec1d, ScaleToClosest>
4         approx(p, mon);
5 auto storage = discretization.computeShapes(approx); 
6 VectorXd E1 = VectorXd::Zero(domain_size);  // vector for containing the current state

We will use Polyharmonic Radial Basis Functions(RBFs) of order 3 augmented with monomials up to 3rd 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$.

1     //Set initial condition 
2     for (int i : discretization.all()){ 
3         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(- 
4         visc*pow(k,2)*t_0)*std::cos(k*discretization.pos(i,0)));

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.

1 for (tt = 1; tt <= t_steps; ++tt) {
2     SparseMatrix<double> M(domain_size, domain_size); // construct matrix of apropreate size
3     VectorXd rhs = VectorXd::Zero(domain_size);  // set empty vector for rhs
4     M.reserve(Range<int>(domain_size, n));
5     auto op = storage.implicitOperators(M, rhs);  // compute the implicit operators

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$ - result of the previous step.

1     // Setting the equation 
2     for (int i : interior) {
3         op.value(i) + (E1[i] * dt) * op.der1(i, 0) + (-dt*visc) * op.der2(i, 0, 0) = E1[i];  
4 }

Boundary conditions are handled in the same way. The rows of the boundary nodes are modified, here with Dirichlet B.C.

1     // Setting boundary conditions
2     for (int i : boundary) {
3         op.value(i) = 0;  
4 }

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.

 1     M.makeCompressed(); // Optimization
 2     BiCGSTAB<SparseMatrix<double, RowMajor>> solver; // Creating solver object
 3     solver.compute(M); // Initialize the iterative solver
 4 
 5     // Solve the linear system and update the state
 6     VectorXd E2 = solver.solve(rhs);
 7     E1 = E2;
 8 
 9     if (tt%(t_factor) == 0) {
10         // Saving state in output file
11         hdf_out.openGroup("/step" + std::to_string(t_save));
12         hdf_out.writeSparseMatrix("M", M);
13         hdf_out.writeDoubleAttribute("TimeStep", t_save);
14         std::cout<< "Time step " << t_save << " of " << t_steps/t_factor << "." << std::endl;
15         hdf_out.writeDoubleAttribute("time",(t_0 + tt * dt));
16         ++t_save;
17         hdf_out.writeDoubleArray("E", E1);
18         }

Results

For our simulation we used the following parameters: \begin{align*} \nu &= 1\\ l &= 1\\ a &= 4\\ b &= 4\\ k &= \frac{4\pi}{2l} \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 order 3.