Hi!

For a project in the Computational Physics Practicum I have to generate orthogonal curvilinear coordinates. As the base of my work I'm using the article: Orthogonal Mapping by G. Ryskin and L. G. Leal, J. Comp. Phys. 50:71-100 (1983), especially their so called weak constraint method (prescribed (x,y) at the boundary). This method can be seen as a variant of Thompson's method with special P and Q. I tried to reproduce their examples in Figs. 7 and 8, but my results won't become orthogonal over the whole grid. At a first glance they look orthogonal, but especially in the upper left and right corners the orthogonality is very poor (If I plot g_{12} vs. \xi,\eta, g_{12} has extrema in these regions). I wonder if someone has tried this method too and could give me some hints.

Rudi

I am not familiar with the sited reference (but can guess what they are probably doing). Specifying both ends of a grid line for orthogonal grids is usually a problem and wants avoiding if possible (but can be done if essential). Relying on source terms to get orthogonality also want avoiding since numerical error in the "bendy" regions can degrade the quality of the grid.

The area, in which I want to generate orthogonal coordinates, is a 2x2 square with the upper line replaced by a cos-function (y=1.6-0.4*cos(pi*x), x=-1..1). Interestingly enough, at this upper boundary the orthogonality is very good. The same is true for the lower boundary. The problems arise at the upper halfes of the side boundaries.

The method described in the mentioned article uses a so called weak constraint. This means that the aspect ratio of the image of a small square, f(\xi,\eta)=\sqrt{ \frac{x_\eta+y_\eta}{x_\xi+y_\xi}} is not constant, but is iterativly calculated. The principle algorithm is 1. Guess initial values for x(\xi,\eta), y{\xi,\eta} (I'm just interpolating between the boundary values of x,y) 2. Calculate f at boundary from x, y 3. Interpolate f 4. Calc new x,y from the PDE (I'm using one step of a SOR scheme) 5. Go to 2.

Rudi

The geometry is fairly simple and seems suitable for a conformal method. If you are able to use a different scheme, I would suggest the following:

- solve the f.d. form of the governing equations directly

- apply the average scale factor ratio between each iteration

- construct orthogonality at the boundary

Methods based on mappings are likely to be a bit cumbersome because of the curved upper boundary.

I am not crazy about the orthogonal curvilinear coordinates. When you have a corner, that is a slope discontinuity, lines normal to the surface from both side of the corner will intersect or diverge depending upon the angle between two surfaces. For a line generated from a surface located away from the corner, it will have more space to turn away from ( or toward ) the line generated from the corner point than the point which is located closer to the corner point. A line generated from a point near the corner point will have to make a very quick turn to avoid collision ( or divergence). In order to see all of these turning at different rates, the mesh size also must change accordingly, that is smaller sizes closer to the corner. As a result, nodal point distributions on the surface also must be non-uniform and concentrated near the corner point. The final results dictate this, regardless of the methods you use. It is easier to see this requirment when you replace a corner by a round arc of small radius. Sometimes, it is very useful to draw the mesh or coordinates by hands.

In the two mentioned examples in the article by Ryskin and Leal, the (x,y)-Nodes were equally distributed along the boundary. They used at least 40x40-Nodes (I counted them in the figures). If I increase this to say 60x60, the orthogonality gets better.

But I would like to reproduce their results. Somehow they must have made it, too.

Assuming that you have done everything right, there are still two possibilities: (1). the computation was carried out on finer mesh, and not every mesh line was printed so that the readers can see the picture more clearly, (2). the mesh plot has a higher order fit. This is sometimes used in the contour plot. At the same time, double check your method step-by-step with a friend. ( It is very difficult to find your own bug because most of time one just keep running his own program in his brain all the time. A friend can spot the problem much easier because his way of thinking is different. A friend by definition is a person who can find your problem and thus you can fix your problem.)