Mesh Generation in 1-D

Question

I am using the method by Brandt in the FAS Multigrid algorithm to estimate the truncation error in a 1-D flow problem and then use that as the basis for generating a mesh(r-adaptation) in 1-D.

The process is standard, in that, transform the equations to a computational domain $\xi$, perform/run computations and map the solution $u$ and $x$ back to the physical domain $x$.

FLOW PROBLEM in computational domain:

$$ u_{\xi\xi} - \left(5 x_\xi + \frac{x_{\xi\xi}}{x_\xi}\right) u_\xi = 0 \quad\tag{1}$$

with boundary conditions $u(\xi=0)=0$ and $u(\xi=1) =1$

MESH EQUATION in computational domain:

$$ x_{\xi\xi} + \frac{2 e_\xi}{e} x_\xi = 0\qquad\tag{2}$$

with boundary conditions $x(\xi=0)=0$ and $x(\xi=1)=1$

where $e$ is the local truncation error from equation-1 and $e_\xi$ its derivative. I am using a simple second order central finite difference scheme for the linear coupled problem.

Things to know:

1. Using the FAS methodology, once an accurate truncation error estimate is made(which requires a reasonably fine grid for the flow problem, say N=65), it is then treated as an known input to equation-2.

2. The analytical solution to equation-1 is : $u(x) = \frac{exp(5 x(\xi)-1)}{exp(5)-1} $ and equation-2 treating the coefficient as known, is smooth that is: $x(\xi)= \frac{1- exp\left[-2\frac{ e_\xi}{e} \xi\right]}{1- exp\left[-2\frac{ e_\xi}{e}\right]}$

3.Truncation error from flow problem:

3.a Coefficient $\frac{e_\xi}{e}$:

4. I have written code in FORTRAN 77

Questions:

1. I have tried using the thomas algorithm for equation-2 and it never seems to give a reasonable solution other than $x(\xi)=0$ for all interior points. Why is this? Can it be that the system is weakly diagonally dominant?

2. I then switched the method to solve equation-2 to a simple Seidel iterative scheme, and that seems to work, but the solution to $x(\xi)$ is not smooth like the analytical counterpart. Why ?

3. Typically what happens is the mesh points would re-distribute accordingly and the flow problem is well resolved. In my case, the solution I obtain for $x(\xi)$ seems to completely change the behavior of the flow problem implying it is obviously not a valid solution.

Any suggestions relevant to any of the above or any insight will be much appreciated.

Answer 1

You have $$ \begin{gathered} u_{xx} - 5u_x = 0 \qquad \text{for} \quad x\in[0,1]\,,\\ u(0) = 0\,,\\ u(1) = 1\,, \end{gathered} $$ which has analytic solution $$u(x) = \frac{\exp(5x)-1}{\exp(5)-1}\,.$$

Using a computation variable $\xi$ such that $x(\xi)$, we can express derivatives of $u$ in terms of $\xi$ $$ \begin{gathered} u_\xi = u_x x_\xi\,,\\ u_{\xi\xi} = u_xx_{\xi\xi} + u_{xx}x_\xi^2\,, \end{gathered} $$ hence can transform the ODE into $$ \begin{gathered} u_{\xi\xi} - \left(5 x_\xi + \frac{x_{\xi\xi}}{x_\xi}\right)u_\xi = 0\,,\\ u(\xi = 0) = 0\,,\\ u(\xi = 1) = 1\,. \end{gathered} $$

The equidistribution principal for a monitor function $M(x)$ (Huang Russell 1997, Generating a non-uniform grid) can be expressed as $$\left(M(x) x_\xi\right)_\xi = 0\,,$$ which I think you've rearranged to give $$x_{\xi\xi} - \frac{M_\xi}{M}x_\xi = 0\,,$$ which will put a high density of points where $M$ is high.

Example

A finite difference solution on a uniform mesh is given by $$ \left(\frac{1}{\Delta x^2}\delta^2_x - \frac{5}{2\Delta x}\delta_{2x}\right)U = \mathbf{b} $$ where $$\mathbf{b} = \left(0,\dotsc,0,\frac{1}{\Delta x^2} - \frac{5}{2\Delta x}\right)^T\,,$$ and $\delta_{2x}$ and $\delta^2_x$ are the central difference operators for the first and second derivatives respectively.

If solving this ODE with finite differences, we could set the monitor function to be the truncation error of the method. Taking a uniform mesh on 20 points, we can plot the error $E_i = |u(x_i) - U|$ seen in the figure below.

(Using any monitor function which goes to zero should ring alarm bells, since we'll get zero mesh density, but I'll carry on for now adding $10^{-8}$.)

Taking $M_i = E_i$, equidistributing to the linear interpolant of $M_i$ and solving with non-uniform finite differences, I got the error shown below.

An example of a smoothing algorithm [Zegeling2007, Chapter 7 of Adaptive Computing:Theory and Algorithms by Tao Tang and Jinchao Xu (2007, Scientific Press)],[Verwer, Blom, Furzeland and Zegeling (1988)], $$ \tilde{M}_i = \sum_{j=0}^N M_j \left(\frac{\sigma}{\sigma+1}\right)^{|i-j|}\,, $$ guarantees that $$\frac{\sigma}{\sigma+1}\leq \frac{x_{i+1}-x_i}{x_i-x_{i-1}} \leq \frac{\sigma+1}{\sigma}\,.$$ Using this with $\sigma=2$, I get the errors shown below.

Not a great advert for adaptivity, but this was just a quick experiment. I did try to iterate this process, using the smoothed error on the non-uniform mesh to define a third mesh etc., which gave lower errors, but this process did not converge.