5
$\begingroup$

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: enter image description here

3.a Coefficient $\frac{e_\xi}{e}$: enter image description here

  1. 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 ?

The mesh points mapped to the original physical domain

  1. 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.

$\endgroup$
  • $\begingroup$ What is your truncation error $e$ looking like? The mesh looks like it's trying to put nearly all the points by $x=1$ $\endgroup$ – Steve Mar 17 '16 at 13:44
  • $\begingroup$ Also I think your analytic solution for $x(\xi)$ is only valid for $\frac{e_\xi}{e} = \text{const.}$. $\endgroup$ – Steve Mar 17 '16 at 13:57
  • $\begingroup$ Hi Mr. @Steve , thank you for your response. I have added in the plot for the truncation error as an edit. I would like to add that the flow problem is a boundary layer type problem, so the maximum truncation error occurs near the x=1 boundary and I expect the points to move towards that area. Also, yes the ratio $e_\xi/e$ will be a constant at every point. I believe I see the confusion now, the analytical solution should have been $(e_\xi/e)_i$, where $i$ $\forall$ $Z \in [1,65]$ at points $\xi_i=0...1$ in increments of $1/64$ for this case. $\endgroup$ – MSIngh Mar 17 '16 at 18:43
  • $\begingroup$ Only if $e_\xi/e$ the same constant over all $i$ (it might well be since the truncation error depends on derivatives of $u$ which is a exponential function). That last point in the truncation error to the right of the turning looks a bit suspicious, that's over one step? That large gradient would drive most of the points into that region, but I'm not sure if that accounts for what looks like all of the points being in $[0.9,1]$. $\endgroup$ – Steve Mar 18 '16 at 9:33
  • $\begingroup$ @Steve, I have recently discovered that the coefficient $\frac{e_\xi}{e}$ is very "zig-zag" , almost like a perfect oscillation, and I suspect that is contributing to the algorithm diverging. I have added the plot of this coefficient to the post. Any suggestions to overcome the smoothness issue? Also, since I am using Finite Differences and the boundary conditions for the error array are homogeneous Dirichlet at both boundaries, the derivative near the boundaries happens to be very high, and that is causing the ridiculous blow ups at the boundaries. No idea how to overcome this. $\endgroup$ – MSIngh Apr 7 '16 at 20:34
3
$\begingroup$

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.

Truncation error for uniform finite difference solution

(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.

Truncation error for non-uniform FD solution with M=E

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.

Truncation error for non-uniform FD solution with M=Etilde

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.

$\endgroup$
  • $\begingroup$ Sorry for the delayed response, in word, FINALS. First of all, I really appreciate your quick effort on this one. I do have a question, how did you get a minus sign when you rearranged the equidistribution principle by Huang & Russell. Secondly, do you think it is possible that using finite differences is causing the "rough" behavior? $\endgroup$ – MSIngh May 18 '16 at 22:31
  • $\begingroup$ In addition, I have been thinking of introducing global constraints on the mesh sizing. Like you said, the monitor function tending to zero is a RED FLAG. Any thoughts on this? I am trying to see how to introduce such a thing that would compensate for the infinite cell size behavior when there seems to be a blow up. $\endgroup$ – MSIngh May 18 '16 at 22:34
  • $\begingroup$ @MSIngh Regarding the minus sign, that looks like an error on my part. I think finite difference should be okay, provided you get a smooth enough mesh. Regarding the monitor function: Negative $M$ is straight up not allowed, check that's not happening. If $M$ is zero anywhere, the smoothing I referenced above should help a bit, but you can achieve a maximum grid-width of $2/N$ by using half the monitor function, and half a uniform mesh, e.g. $\tilde{M}(x) = \frac12 (M(x) + \theta)$ where $\theta = \int_0^1M(x) \mathrm{d}x$. $\endgroup$ – Steve May 19 '16 at 10:40
  • $\begingroup$ For my mesh equidistribution I used a piecewise linear monitor function and integrated it to find $x_i$ as in this answer which I linked above, since I already had the code. Not sure how finite difference is for solving the mesh equation $\endgroup$ – Steve May 19 '16 at 10:52
  • $\begingroup$ I have also read that mesh motion can be controlled by parabolizing elliptic grid equation(such as the one it this thread), which is more mathematical than using a smoothing operation separately. Have you ever read about this or have any idea about this technique? Also, is there an easier way to chat in this forum, it would make communicating back and forth much faster and efficient. $\endgroup$ – MSIngh May 25 '16 at 6:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.