Skip to content

Schwarz Algorithms for Chebyshev Pseudospectral Method

Recently I bacame interested in Domain Decomposition algorithms for spectral methods, and in particular for Chebyshev collocation method.

There has been a lot of activity in this field since the second half of the 80′s. It will be hard to list all the papers, but maybe in some future post I will make a short list of some interesting ones.

Thanks to a funtastic scientific python blog I was able to get my hands on a masterfully written FD code implementing Multiplicative Schwarz method for Poisson equation on a square. With minor changes in the routine evaluating the 2D Laplacian, I was able to adapt the code for Chebyshev pseudospectral method. The results were wonderful!

I made two examples, and implemented them in another-chebpy. Both examples are in my google code git repository. One is solving Poisson eq. and other solves Biharmonic BVP (clamped BCs). I uploaded a small video on youtube that shows (rather slow) convergence for the case of biharmonic equation.

This can be accelerated by multigrid, which will be done in the future. Also parallelization is possible, but in this cae multicoloring (available here already) is required to decouple domains.

Surfaces defined by tensor product of Bernstein polynomials

Following equations will help us in developing a collocation method using Bernstein polynomials for two-dimensional problems.

Let the following equation represent a surface in \mathbb{R}^3 defined by a tensor product of Bernstein basis functions:

f(x,y)=\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}B_{i,n}(x)B_{j,m}(y)

We can define p-th order partial derivative in x -axis direction in a following way

\frac{\partial^pf(x,y)}{\partial x^p}=\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}\frac{\partial^p B_{i,n}(x)}{\partial x^p}B_{j,m}(y)

and an q-th order derivative in y -axis direction

\frac{\partial^qf(x,y)}{\partial y^q}=\sum_{i=0}^n \sum_{j=0}^m\beta_{i,j}B_{i,n}(x)\frac{\partial^qB_{j,m}(y)}{\partial y^q}

Finally mixed derivative of order p+q is defined as follows:

\frac{\partial^{p+q}f(x,y)}{\partial x^p y^q}=\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}\frac{\partial^pB_{i,n}(x)}{\partial x^p}\frac{\partial^qB_{j,m}(y)}{\partial y^q}.

Laplacian operator can be written like this

\Delta f(x,y)=\frac{\partial^2 f(x,y)}{\partial x^2 }+\frac{\partial^2 f(x,y)}{\partial y^2}

=\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}B_{j,m}(y)\frac{\partial^2B_{i,n}(x)}{\partial x^2}+\sum_{i=0}^n\sum_{j=0}^m\beta_{i,j}B_{i,n}(x)\frac{\partial^2 B_{j,m}(y)}{\partial y^2},

which can further be simplified to

\Delta f(x,y)=\sum_{i=0}^n \sum_{j=0}^m\beta_{i,j}\left[B_{j,m}(y)\frac{\partial^2 B_{i,n}(x)}{\partial x^2}+B_{i,n}(x)\frac{\partial^2 B_{j,m}(y)}{\partial y^2}\right].

Aside

What did Boris Grigorievich Galerkin say?

“I want you, not to be able to make residual with the same basis functions, you used to create the solution.”

(P.S. This story is completely not true, and I urge my readers to find a better explanation of (Bubnov-) Galerkin method, if it exists.)

Notes on implementation of Realizable k-epsilon model

FLUENT users are aware of the “realizable” k-\epsilon turbulence model. For more experienced modelers this means an implementation of a model presented in the paper by T.H. Shih, W.W. Liou,  A. Shabbir, Z. Yang and J. Zhu, “A New k-\epsilon Eddy Viscosity Model for High Reynolds Number Turbulent Flows”, Computers and Fluids, Vol. 24 No. 3, pp. 227-238, 1995.

In this model, realizability is achieved by letting  C_\mu coefficient appearing in the formulation of eddy viscosity, to be no longer a constant. Experiments also suggest it shouldn’t be a constant. While it has a value of 0.09 in inertial sublayer of BL flow, in homogeneous shear flow it is around 0,05. Therefore some way of changing the C_\mu coefficient depending on the flow situation is desirable.

Let’s first say what realizability means mathematically, and then we’ll go into details about  C_\mu.

The realizability constraints any model should meet are following inequalities regarding Reynolds stress tensor components:

\overline{u_\alpha^2}\textgreater0,\,(\alpha = 1,2,3)

\frac{\overline{u_\alpha u_\beta}^2}{\overline{u_\alpha^2} \, \overline{u_\beta ^2}} \leq 1, (\alpha = 1,2,3; \beta = 1,2,3)

Eq. (2) is know as Schwarz’s inequality. So we need positivity of normal stresses (Eq.(1)), and we need Reynolds stress tensor components to obey Schwartz’s inequality.  Based on these realizability constraints W. Reynolds as well as Shih, Zhu and Lumley  proposed following form of an equation for C_\mu:

C_\mu=\frac{1}{A_0+A_sU^{(*)}\frac{k}{\epsilon}}

Shih, Zhu, Lumley give this formulation for the constituents of the previous equation

U^{(*)}=\sqrt{S_{ij}S_{ij}+\tilde{\Omega}_{ij}\tilde{\Omega}_{ij}},

Where

\tilde{\Omega}_{ij}=\Omega_{ij}-2\epsilon_{ijk}\omega_k

\Omega_{ij}=\overline{\Omega_{ij}}-\epsilon_{ijk}\omega_k

where \overline{\Omega_{ij}} is the mean rate-of-rotation tensor viewed in a rotating reference frame with the angular velocity \omega_k. So this model can be applied in cases like rotating pipe etc.

A_s=\sqrt{6}cos\phi, \phi=\frac{1}{3} arccos(\sqrt{6}W), W=\frac{S_{ij}S_{jk}S_{ki}}{\tilde{S}^3}, \tilde{S}=\sqrt{S_{ij}S_{ij}}.

Only undetermined is coefficient A_0. As far as I’m aware all authors treat it as a constant, and calibrate it in a way that it returns desired values in some cases (like 0.09 in inertial sublayer in BL flow). In original reference authors came to the value of A_0=4.0, while in FLUENT this model is implemented with value A_0=4.04.  I look forward to some study where A_0 might not be a constant.

Now going to transport equations for k and \epsilon . The first one (k-equation) remains the same as in standard k-\epsilon model, while there are some changes with equation for dissipation. Here I’ll show them both for the sake of completeness

\frac{\partial}{\partial t} (\rho k) + \frac{\partial}{\partial x_j} (\rho k u_j) =

\frac{\partial}{\partial x_j} \left [ \left(\mu + \frac{\mu_t}{\sigma_k}\right) \frac{\partial k} {\partial x_j} \right ] + P_k + P_b - \rho \epsilon - Y_M + S_k

\frac{\partial}{\partial t} (\rho \epsilon) + \frac{\partial}{\partial x_j} (\rho \epsilon u_j) =

\frac{\partial}{\partial x_j} \left[ \left(\mu + \frac{\mu_t}{\sigma_{\epsilon}}\right) \frac{\partial \epsilon}{\partial x_j} \right ] + \rho \, C_1 S \epsilon - \rho \, C_2 \frac{{\epsilon}^2} {k + \sqrt{\nu \epsilon}} + C_{1 \epsilon}\frac{\epsilon}{k} C_{3 \epsilon} P_b + S_{\epsilon}

Where

C_1 = \max\left[0.43, \frac{\eta}{\eta + 5}\right] , \;\;\;\;\; \eta = S \frac{k}{\epsilon}, \;\;\;\;\; S =\sqrt{2 S_{ij} S_{ij}}

Production of turbulence kinetic energy due to velocity gradients P_k, as well as production of TKE due to buoyancy P_b, are treated in the same way as in standard k-\epsilon model.

When we find both k and \epsilon fields we move to computation of eddy viscosity in the standard form for k-\epsilon model

\mu_t=\rho C_\mu \frac{k^2}{\epsilon}

The purpose of this note is to point out a few things I found useful for stable implementation of the model.

Argument of the arccos function should be in interval [-1,1], therefore I suggest implementation of \phi in a following way:

\phi=\frac{1}{3}, \;\;arccos(max(-1, min(\sqrt{6}W, 1)))

Otherwise you might get  a few NaN’s in the C_\mu field which usually ruin the solution.

Experiments suggest there is no reason to let C_\mu be larger than 0.09, so we will limit it from above

C_\mu =min(C_\mu, 0.09)

It would be useful to see what are the lowest values in the C\mu field and check with some experimental data, just to make sure we are on the right track.

Finally constants of the model are:

C_{1 \epsilon} = 1.44, \;\; C_2 = 1.9, \;\; \sigma_k = 1.0, \;\; \sigma_{\epsilon} = 1.2

I wish you happy computing!

Image

Chebyshev pseudospectral solution to Poisson equation

Chebyshev pseudospectral solution to Poisson equation

I was interested to see how mlab handles plots from my Chebyshev pseudospectral code. The result was quite surprising. The light effects give nice impression of this object’s position in space. The colors you might chose from are the same as in Matplotlib, but here they change smoothly – the one thing I was missing there.
I also posted this figure to figshare here.

Follow

Get every new post delivered to your Inbox.