nsCouette-未完。。。

nsCouette介绍

The Taylor-Couette (TC) set-up is one of the most famous paradigmatic systems for wall-bounded shear flows in general and maybe the most important one if you are interested in rotating shear flows in particular. For understanding the following chapters of this documentation, it might help to acquaint yourself with a bit of crucial terminology and basic concepts of the TC system as well as with the notation we will use throughout this documentation and in the source code of nsCouette.

• 代码来源：

John F. Gibson提出Plane Couette算法 http://channelflow.org/dokuwiki/doku.php?id=docs

https://github.com/epfl-ecps/channelflow

• 参考AVILA课题组：

NONLINEAR DYNAMICS OF MODE COMPETITION IN ANNULAR FLOWS- 2008-MARC AVILA-phd-thesis

Direct Numerical Simulation of Transition to Turbulence and Turbulence Control in Pipe Flow- Baofang Song-phd-thesis

Governing equation and geometry¹

Start from N-S function²

The naiver-Stoke model are very similar to the corresponding equation of the Euler equation. The different is additional terms that take into count the effect of internal friction and heat conduction in the fluid.

The momentum equation in this model is of the form of the Cauchy momentum equation:

$${D \vec{u}\over Dt} = {1\over \rho} {\partial \sigma_{\alpha \beta} \over \partial x_\beta} + \vec{g}$$

$\sigma$ is the stress tensor, g contains all the body forces per unit mass.

Write in terms of equation in cylindrical coordinates:

$$\vec{u}(r,\theta, z, t)= u_r (r,\theta,z,t) \vec{e_r} + u_\theta(r,\theta,z,t)\vec{e_\theta} + u_z(r,\theta,z,t)\vec{e_z}$$

where $vec{e}$ is a right-handed triad of unit vectors.

Material Derivative

$${D \vec{u}\over Dt} = \partial_t \vec{u} + \vec{u} \cdot \nabla \vec{u}$$

Cauchy stress tensor³

$$\sigma = -p \vec{I} +\sigma'\\\\ \sigma_{\alpha\beta}'=\mu({\partial u_\beta \over \partial x_\alpha}+{\partial u_\alpha\over x_\beta}-{2\over 3} \sigma_{\alpha\beta}{\partial u_\gamma \over \partial x_\gamma}) + \mu_B \sigma_{\alpha\beta}{\partial u_\gamma \over \partial x_\gamma}$$ $$\nu={\mu \over \rho}$$

其中，$\sigma$ is the tensor.

$$\nabla \vec{u}=[{\partial u_\beta \over \partial x_\alpha}]=\vec{u} \otimes \nabla=(u_r \vec{e_r}+u_\theta \vec{e_\theta}+ u_z \vec{e_z})\otimes (\vec{e_r} {\partial \over \partial r} + \vec{e_\theta}{1\over r}{\partial \over \partial \theta} + \vec{e_\theta}{\partial \over \partial z})$$

one term for example:

$${\partial \vec{e_r}\over \partial r}= {\partial \vec{e_\theta}\over \partial r} ={\partial \vec{e_z}\over \partial r}=0\\\\ {\partial \vec{e_r} \over \partial \theta} =\vec{e_\theta}\\\\ {\partial \vec{e_\theta} \over \partial \theta} =-\vec{e_r}$$ $${1\over r}{\partial \over \partial \theta} (u_r \vec{e_r}) \otimes \vec{e_\theta} = {1\over r}{\partial u_r \over \partial \theta} \vec{e_r} \otimes {\vec{e_\theta}} + {u_r \over r} {\partial \vec{e_r}\over \partial \theta} \otimes \vec{e_\theta}= {1\over r}{\partial u_r \over \partial \theta} \vec{e_r} \otimes {\vec{e_\theta}} + {u_r \over r} \vec{e_\theta} \otimes \vec{e_\theta}$$

The final result, expressed as a tensor, is

$$\nabla \vec{u}=\begin{bmatrix} {\partial u_r\over \partial r} & {1\over r}{\partial u_r \over \partial \theta}-{u_\theta\over r} & {\partial u_r \over \partial z} \\ {\partial u_\theta \over \partial r} & {1\over r}{\partial u_\theta \over \partial \theta}+{u_r\over r} & {\partial u_\theta \over \partial z} \\ {\partial u_z\over \partial r} & {1\over r}{\partial u_\theta \over \partial \theta} & {\partial u_z\over \partial z} \\ \end{bmatrix}$$

Thus,

$$\nabla \vec{u}+\nabla \vec{u}^T=[{\partial u_\beta \over \partial x_\alpha}+{\partial u_\alpha\over x_\beta}]_{tensor}=\begin{bmatrix} 2{\partial u_r\over \partial r} & {1\over r}{\partial u_r \over \partial \theta}-{u_\theta\over r}+{\partial u_\theta \over \partial r} & {\partial u_r \over \partial z} + {\partial u_z\over \partial r}\\ {1\over r}{\partial u_r \over \partial \theta}-{u_\theta\over r}+{\partial u_\theta \over \partial r} & 2({1\over r}{\partial u_\theta \over \partial \theta}+{u_r\over r}) & {\partial u_\theta \over \partial z} + {1\over r}{\partial u_\theta \over \partial \theta} \\ {\partial u_r \over \partial z} + {\partial u_z\over \partial r} & {\partial u_\theta \over \partial z} + {1\over r}{\partial u_\theta \over \partial \theta}& 2{\partial u_z\over \partial z} \\ \end{bmatrix}$$

最后， 可以简化stress term$\nabla \cdot \sigma$ in N-S equation:³

$$\nabla \cdot \sigma_{tensor}= ({\sigma_{rr}\over r}+{\partial \sigma_{rr}\over\partial r}+ {1\over r}{\partial \sigma_{r\theta}\over \partial \theta}+ {\partial \sigma_{rz}\over \partial z}-{\sigma_{\theta\theta}\over r})\vec{e_r}\\\\ +({2 \sigma_{r\theta}\over r}+{\partial \sigma_{r\theta}\over \partial r}+ {1\over r}{\partial \sigma_{\theta\theta}\over \partial\theta}+{\partial \sigma_{\theta z}\over\partial z})\vec{e_\theta} \\\\ + ({\partial \sigma_{rz}\over r}+ {\partial \sigma_{r z}\over \partial r}+{1\over r}{\partial \sigma_{\theta z}\over \partial \theta}+{\partial \sigma_{zz}\over \partial z})\vec{e_z}$$

带入上面的公式，

第一项-这里有一大堆整理，参考³

$$({\sigma_{rr}\over r}+{\partial \sigma_{rr}\over\partial r}+ {1\over r}{\partial \sigma_{r\theta}\over \partial \theta}+ {\partial \sigma_{rz}\over \partial z}-{\sigma_{\theta\theta}\over r})\vec{e_r} \\\\ = -{\partial p \over \partial r} + \mu （-{u_r \over r^2}+{\partial \over \partial r }[{1\over r}{\partial (r u_r)\over \partial r}+{1\over r}{\partial u_\theta \over \partial \theta}+{\partial u_z \over \partial z}]+{1\over r}{\partial \over \partial r}(r{\partial u_r \over \partial r})+{1\over r^2}{\partial ^2 u_r \over \partial \theta^2}+ {\partial^2 u_r \over\partial z^2}-{2\over r^2}{\partial u_\theta \over \partial \theta}）\vec{e_r}$$

For an incompressible fluid, the term in the square brackets is equal to zero because of continuity equation.

同理，第二项，第三项

Incompressible N-S equation

In may subfields of fluid mechanics, the fluid is often considered to be incompressible, meaning that the density $\rho$ is constant. This assumption satisfied the mass conservation equation to:

$$\nabla \cdot \vec{u}=0$$

Consequently, the bulk viscosity in N-S equation becomes irrelevant, as its terms is zero. For this reason, bulk viscosity is often neglected in fluid mechanics. Still, it relevant in acoustics and high-velocity compressible flow.

同时，2/3那一项也忽略掉。

建立不可压N-S方程, 这里假设$\rho$ is constant

$$\partial_t \vec{u} + \vec{u} \cdot \nabla \vec{u} = -{1\over \rho}\nabla p^h +\nu \nabla \cdot({\partial u_\beta \over \partial x_\alpha}+{\partial u_\alpha \over \partial x_\beta}）$$

其中，$\vec{u}(\vec{r},t)$ 是速度场， $p^h$ 是hydrodynamic pressure. 雷诺数$Re= \Omega r d/ \nu$, $d= r_0-r_i$, $\eta= r_i / r_o$, the length-to-gap aspect ratio $\Gamma = L_z/d$, $L_z$ 为柱子的高度。

—》由上面的详细推导过程，阐述了 $\nabla \cdot \sigma$的详细推导过程：

因此，公式最终可以整理成不可压的形式为（一些细节和约化参见附录1）：

$${\partial u_r \over \partial t} + u_r {\partial u_r \over \partial r} +{u_\theta \over r}{\partial u_r \over \partial \theta} - {u_\theta^2 \over r}+u_z{\partial u_r \over \partial z}=\\\\ -{1\over\rho}{\partial p \over \partial r} + {\mu\over \rho}(-{u_r \over r^2}+ {1\over r}{\partial \over \partial r}(r{\partial u_r \over \partial r})+{1\over r^2}{\partial^2 u_r \over \partial \theta^2}+ {\partial^2 u_r \over \partial z^2} - {2\over r^2}{\partial u_\theta \over \partial \theta}) + g_r$$ $${\partial u_\theta \over \partial t} + u_r {\partial u_\theta \over \partial r} +{u_\theta \over r}{\partial u_\theta \over \partial \theta} + {u_r u_\theta \over r}+u_z{\partial u_\theta \over \partial z}=\\\\ -{1\over r \rho}{\partial p \over \partial \theta} + {\mu\over \rho}(-{u_\theta\over r^2}+ {1\over r}{\partial \over \partial r}(r{\partial u_\theta \over \partial r})+{1\over r^2}{\partial^2 u_\theta \over \partial \theta^2}+{\partial^2 u_\theta \over \partial z^2}+ {2\over r^2}{\partial u_r \over \partial \theta}) + g_\theta$$ $${\partial u_z \over \partial t} + u_r {\partial u_z \over \partial r} +{u_\theta \over r}{\partial u_z \over \partial \theta} +u_z{\partial u_\theta \over \partial z}=\\\\ -{1\over\rho}{\partial p \over \partial z} + {\mu\over \rho}({1\over r }{\partial \over \partial r}(r {\partial u_z \over \partial r})+{1\over r^2}{\partial^2 u_z \over \partial \theta^2}+{\partial^2 u_z \over \partial z^2}) + g_z$$

进一步，我们对其作dimensionless处理，通过d约化，$\tau=d^2/\nu$, 和 $\nu^2 / d^2$ 作为空间，时间的计量单位，压力通过$p=p^h / \rho$. 这样上述方程转化为，即为文献¹中的公式（3）：

$$\partial_t \vec{u} + (\vec{u} \cdot \nabla) \vec{u} =-\nabla p + \nabla \cdot \sigma \\\\ \nabla \cdot \vec{u} =0$$ $${\partial u_r \over \partial t} + u_r {\partial u_r \over \partial r} +{u_\theta \over r}{\partial u_r \over \partial \theta} - {u_\theta^2 \over r}+u_z{\partial u_r \over \partial z}=\\\\ -{\partial p \over \partial r} + (-{u_r \over r^2}+ {1\over r}{\partial \over \partial r}(r{\partial u_r \over \partial r})+{1\over r^2}{\partial^2 u_r \over \partial \theta^2}+ {\partial^2 u_r \over \partial z^2} - {2\over r^2}{\partial u_\theta \over \partial \theta}) + g_r$$ $${\partial u_\theta \over \partial t} + u_r {\partial u_\theta \over \partial r} +{u_\theta \over r}{\partial u_\theta \over \partial \theta} + {u_r u_\theta \over r}+u_z{\partial u_\theta \over \partial z}=\\\\ -{1\over r}{\partial p \over \partial \theta} + (-{u_\theta\over r^2}+ {1\over r}{\partial \over \partial r}(r{\partial u_\theta \over \partial r})+{1\over r^2}{\partial^2 u_\theta \over \partial \theta^2}+{\partial^2 u_\theta \over \partial z^2} +{2\over r^2}{\partial u_r \over \partial \theta}) + g_\theta$$ $${\partial u_z \over \partial t} + u_r {\partial u_z \over \partial r} +{u_\theta \over r}{\partial u_z \over \partial \theta} +u_z{\partial u_\theta \over \partial z}=\\\\ -{\partial p \over \partial z} + ({1\over r }{\partial \over \partial r}(r {\partial u_z \over \partial r})+{1\over r^2}{\partial^2 u_z \over \partial \theta^2}+{\partial^2 u_z \over \partial z^2}) + g_z$$

• 备注：后续可进一步简化为轴对称，去掉$\theta$ direction. 不再进行赘述。

Numerical simulation

空间离散

1. Fourier-Galerkin method: 相当于每个时间步下，都要求解一次方程，算出不同l、n模态或者波数下所对应的r方向的chebshev系数。
   • 备注：这里面也许可以做一些文章。稳定性问题判定，通过eigenvalue。

$$\vec{u}(r,\theta, z)=\sum_{l=-L}^L \sum_{n=-N}^N \hat{u}(r)e^{i (lk_z z+n k_\theta \theta)}\\\\ \vec{p}(r,\theta, z)=\sum_{l=-L}^L \sum_{n=-N}^N \hat{p}(r)e^{i (lk_z z+n k_\theta \theta)}\\\\$$

其中， $kz$ 是最小波数， $k\theta$ 是简化的周期数。

$$\Gamma=2\pi/k_z\\\\ L_0=2\pi / k_\theta$$

将上面的离散形式带入到最终的推导公式，转化到每个(l,n)：

$$\partial_t \hat{u}_t + \hat{N}_r = -\partial_r \hat{p} + \hat{\triangle} \hat{u}_r - \hat{u}_r/r^2 - 2inlk_\theta\hat{u}_\theta/r^2\\\\ \partial_t \hat{u}_\theta + \hat{N}_\theta = -ink_\theta \hat{p}/r +\hat{\triangle}\hat{u}_\theta -\hat{u}_\theta/r^2 - 2ink_\theta \hat{u}_r/r^2 \\\\ \partial_t \hat{u}_z +\hat{N}= -ilk_z \hat{p} + \hat{\triangle} \hat{u}_z$$

其中，非线性项包含$u\cdot u$和$ u^2$项。

$$\hat{\triangle}=\partial_r/r +\partial_{rr}-n^2k_\theta^2/r^2-l^2k_z^2$$

接下来做的是解耦径向和周向速度：

带入：

$$\hat{u}_+=\hat{u}_r + i\hat{u}_\theta\\\\ \hat{u}_+=\hat{u}_r - i\hat{u}_\theta$$

得到：

$$\partial_t \hat{u}_+(r)+\hat{N}_+(r)= -\partial_r \hat{p}(r)+ nk_\theta \hat{p}(r)/r +(\hat{\triangle}-1/r^2 -2nk_\theta/r^2)\hat{u}_+ \\\\ \partial_t \hat{u}_-(r)+\hat{N}_-(r)= -\partial_r \hat{p}(r)- nk_\theta \hat{p}(r)/r +(\hat{\triangle}-1/r^2 +2nk_\theta/r^2)\hat{u}_- \\\\ \partial_t \hat{u}_z +\hat{N}= -ilk_z \hat{p} + \hat{\triangle} \hat{u}_z$$ $$\hat{N}_{\pm}= \hat{N}_r \pm i \hat{N}_\theta$$

1. finite-difference method in r方向的Chebyshev离散⁵

Use a standard high-order, central finite-difference method to approximate the radial derivatives.

$$r_j={1+\eta\over 2(1-\eta)}+{sin^{-1}(-\alpha cos(\pi j/M))\over 2 sin^{-1}\alpha}, j=0,...,M$$

1. 时间离散格式-Adams-Bashforth back-ward-difference

   $${3u^{i+1}-4u^i + u^{i-1}\over 2\triangle t} + 2 N^i(u) - N^{i-1}(u) = -\nabla p^{i+1}+\triangle u^{i+1}$$

   • 计算频域的非线性项系数$\hat{N}^i(u)$, 通过3/2-dealiasing rule:

     • Do matrix-vector multiplication to calculate $\partial_r\hat{u}^i$ (FD method) - chebfun: diff
     • ??Compute dot product in Fourier space to calculate $\partial_\theta \hat{u}^i$ and $\partial_z \hat{u}^i$
     • Perform Fourier transform of $\partial_{r,\theta,z}\hat{u}^i$ and $\hat{u}^i$ to obtain the velocity field and all its derivatives in physical space.
     • Calculate N^i(u)=u^i\cdot \triangle u^i$.
     • Perform inverse Fourier transform to obtain the spectral coefficients $\hat{N}^i(u)$.

   • Obtain the pressure prediction, $\hat{p}*$: solve the Possion equation, with consistent Neumann boundary condition.

     $$\triangle \hat{p}*= \nabla \cdot [-2 \hat{N}^i(u)+\hat{N}^{i-1}(u)]$$

   • Obtain the velocity prediction, $\hat{u}*$: solve the three three Helmholtz equations

     $${3\hat{u}*-4\hat{u}^i+ \hat{u}^{i-1}\over 2\triangle t}+ 2\hat{N}^i(u) -\hat{N}^{i-1}(u)=-\nabla\hat{p}* + \triangle \hat{u}*$$

     with Dirichlet boundary conditions.

   • Correct via an intermediate variable $\phi=2\triangle t (\hat{p}^{i+1}-\hat{p}*)/3$. The incompressibility condition $\nabla\cdot \hat{u}^{i+1}=0$ leads to a Possion equation for $\phi$ with homogeneous Neumann boundary conditions:

     $$\triangle\phi =\nabla \cdot \hat{u}*\\\\ \partial_r \phi_{r_{i,o}}=0$$

   • Compute pressure and velocity correction, $\hat{p}^{i+1}$ and $\hat{u}^{i+1}$:

     $$\hat{p}^{i+1}=\hat{p}* + 3\phi /(2\triangle t)\\\\ \hat{u}^{i+1}=\hat{u}*-\nabla\phi$$

   • Go back to step 1

   The Navier-Stokes equations are thus advanced in time by solving five systems of linear euqation, of Possion or Helmholtz type, for each Fourier mode.

附录1 ：

在柱坐标系下：

$$h_r=1, h_\theta=r; h_z=1, \\\\ \vec{u}=[u_r; u_\theta; u_z;]$$

令：

$$\nabla \cdot \vec{u}={1\over h_1 h_2 h_3}[{\partial(v_1 h_2 h_3)\over \partial w_1}+{\partial(v_2 h_3 h_1)\over \partial w_2}+{\partial(v_3 h_1 h_2)\over \partial w_3}] \\\\ = {1\over r}[{\partial u_z r \over \partial z} + {\partial (u_r r) \over \partial r} + {\partial u_\theta\over \partial \theta}]={\partial u_z \over \partial z} + {\partial u_r \over \partial r} + {u_r \over r}+ {1\over r}{\partial u_\theta\over \partial \theta}\\\\$$ $$D={v_1\over{h_1}}{\partial \over {\partial \xi_1}}+{v_2\over{h_2}}{\partial \over {\partial \xi_2}}+{v_3\over{h_3}}{\partial \over {\partial \xi_3}}=u_r {\partial \over {\partial r}}+{u_\theta \over r}{\partial \over {\partial \theta}}+u_z{\partial \over {\partial z}}$$

因:

$$(\vec{u}\cdot\nabla)\vec{u}=\begin{cases} term1:D v_1+{v_2\over{h_2h_1}}(v_1{\partial h_1\over \partial \xi_2}-v_2{\partial h_2\over \partial \xi_1})+{v_3\over{h_3h_1}}(v_1{\partial h_1\over \partial \xi_3}-v_2{\partial h_3\over \partial \xi_1})\\ term2:D v_2+{v_3\over h_3h_2}(v_2{\partial h_2 \over \partial \xi_3}-v_3{\partial h_3 \over \partial \xi_2})+{v_1\over h_1h_2}(v_2{\partial h_2 \over \partial \xi_1}-v_1{\partial h_1 \over \partial \xi_2})\\ term3:D v_3+{v_1\over h_1h_3 }(v_3{\partial h_3 \over \partial \xi_1}-v_1{ \partial h_1 \over \partial \xi_3})+{v_2\over h_2h_3}(v_3{\partial h_3\over \partial \xi_2}-v_3{\partial h_2 \over \partial \xi_3})\\ \end{cases} \\\\$$

得：

$$(\vec{u}\cdot\nabla)\vec{u}=\begin{cases} term1: u_r {\partial u_r \over {\partial r}}+{u_\theta \over r}{\partial u_r \over {\partial \theta}}+u_z{\partial u_r \over {\partial z}} - u_\theta^2 / r\\ term2:u_r {\partial u_\theta \over {\partial r}}+{u_\theta \over r}{\partial u_\theta \over {\partial \theta}}+u_z{\partial u_\theta \over {\partial z}} + u_ru_\theta / r \\ term3:u_r {\partial u_z \over {\partial r}}+{u_\theta \over r}{\partial u_z \over {\partial \theta}}+u_z{\partial u_z \over {\partial z}} +0 \\ \end{cases} \\\\$$

最终得到：

$$(\partial_t + \vec{u}\cdot \nabla )u_r - u_\theta ^2 / r=-\partial_r p + \triangle u_r - u_r / r^2 -2 \partial_\theta u_\theta / r^2 \\\\ (\partial_t + \vec{u}\cdot \nabla) u_\theta + u_\theta u_r /r = -\partial_\theta p/ r +\triangle u_\theta - u_\theta /r^2 + 2 \partial_\theta u_r /r^2 \\\\ (\partial_t + \vec{u}\cdot \nabla)u_z = -\partial_z p + \triangle u_z \\\\ u_r/r +\partial_r u_r +\partial_\theta u_r +\partial_zu_z =0$$ $$\nabla p = \partial_r p \vec{e_r} + {1\over r}\partial_\theta p \vec{e_\theta} +\partial_z p \vec{e_z}$$

Boundary condition

$$u_\theta (r_{i,o},\theta,z)=Re_{i,o}\\\\ u_{r,z}(r_{i,o},\theta,z)=0\\\\ \vec{u}(r,\theta,z)=\vec{u}(r,\theta+L_\theta,z)\\\\ \vec{u}(r,\theta,z)=\vec{u}(r,\theta,z+\Gamma)\\\\$$

压力项可以用Possion方程进行近似处理：

$$\triangle p =- \nabla\cdot N(u), where N(u)=\vec{u} \cdot \nabla \vec{u}$$

需要满足一定的边界条件：⁴

$$\partial_n p_{r_{i,o}}=\vec{n}\cdot [-\partial_t \vec{u}-N(\vec{u})+\nabla\cdot \sigma]$$

Install

环境：ubuntu 20.04

软件有三个branch，其中nsCouette-gpu 通过GPU计算，支持nivida的cuda-c。

具体操作是，主要需要按照库fftw，hdf5，libblas3, liblapack-dev

$ git clone http://github.com/dfeldmann/nsCouette.git
$ git clone -b nsPipe-1.0 https://github.com/dfeldmann/nsCouette.git
$ git clone -b nsCouette-gpu https://github.com/dfeldmann/nsCouette.git
$ cd nsCouette
$ sudo apt install libblas3 liblapack-dev
$ sudo apt install libhdf5-openmpi-dev
$ sudo apt install libfftw3-dev
# 需要进入ARCH文件夹中，修改make.arch.gcc-linux中hdf5的路径，
# LIBHDF5 = -L/usr/lib/x86_64-linux-gnu/hdf5/openmpi -lhdf5_fortran -lhdf5 -lz //-L表示lib库
# INCLUDES = -I/usr/include/hdf5/openmpi //-I表示include
$ make ARCH=gcc-linux
$ cd gcc-linux && ./nsCouette.x

nsCouette-gpu环境配置(类似C语言)：

# 首先得有Nivida的显卡，安装好cuda环境
# 修改Makefile
# all: $(GPU_OBJECTS)
#        $(NVCC) -o taylorC $(GPU_OBJECTS) $(LIBS) -L/usr/lib/x86_64-linux-gnu/hdf5/openmpi 
#$ (GPU_OBJECTS): %.o: %.cu
#        $(NVCC) -c   $< -o $@  -I/usr/include/hdf5/openmpi -I/usr/lib/x86_64-linux-gnu/openmpi/include
# 该修改根据自己openmpi和hdf5的安装路径进行选择，如果按照上面库的安装，该路径如上所示。
# 本机只有一个独立显卡，因此还需要修改main.cu 函数的”dev=0“。
$ make clean && make 

代码解析

// main.cu 代码
#include"TayC.h"

int main( int argc, const char* argv[]){
  
  printf("\n+++++++++++++++++++++++++++");
  printf("\nStarting GPU Taylor-Couette");
  printf("\n+++++++++++++++++++++++++++");
  
  int dev=0; //根据自己主机的GPU数目来定
  double time;

  printf("\n");
  printf("\nSetting device: %d",dev);
  CHECK_CUDART(cudaSetDevice(dev));
  
  //Timing variables
  clock_t start, end;
  double gpu_time_used;

  //Set up
  size_p sizes;
    
  sizes.Nr=NR; //R方向的模态数
  sizes.Nt=NT; //theta方向的模态数
  sizes.Nz=NZ; //z方向的模态数
    
  //Initialize the mesh  
  double* grid_mesh=(double*)malloc(NR*sizeof(double));

  double r_i = 1.0;         //ETA/(1.0-ETA);
  double r_o = 1.0/ETA;     //1.0/(1.0-ETA); 
  
  
  //Chebyshev mesh in r
  
  for(int j=0;j<NR;j++){
      //grid_mesh[j]=2.0*(double)j/(double)(NR-1)+1.0;
        grid_mesh[j]= (r_i+r_o)/2.0 - cos(PI2/2.0*j/(NR-1))/2.0;
    
  }
  
  //Write mesh

  //Set modules
  setImplic_7_exp(grid_mesh,sizes);
  setDeriv_9_exp(grid_mesh,sizes);
  setBoundary(grid_mesh);
  setNonlinear(grid_mesh);
  setFft(sizes);
  setCublas();
  setIntegrator(sizes);
  setLinear(grid_mesh);
  setStatistics(grid_mesh);
  
  //Allocate memory buffers
  vfield u, uw, rhs, rhsw;
  
  size_t size_p=NR*NT*NZ*sizeof(double2);

  CHECK_CUDART(cudaMalloc(&u.r,size_p));
  CHECK_CUDART(cudaMalloc(&u.t,size_p));
  CHECK_CUDART(cudaMalloc(&u.z,size_p));
  
  CHECK_CUDART(cudaMalloc(&uw.r,size_p));
  CHECK_CUDART(cudaMalloc(&uw.t,size_p));
  CHECK_CUDART(cudaMalloc(&uw.z,size_p));
  
  CHECK_CUDART(cudaMalloc(&rhs.r,size_p));
  CHECK_CUDART(cudaMalloc(&rhs.t,size_p));
  CHECK_CUDART(cudaMalloc(&rhs.z,size_p));
  
  CHECK_CUDART(cudaMalloc(&rhsw.r,size_p));
  CHECK_CUDART(cudaMalloc(&rhsw.t,size_p));
  CHECK_CUDART(cudaMalloc(&rhsw.z,size_p));

  
  //Start initial field
  initField(u, grid_mesh);
  
  //Or read check point
  //readCheckpoint(u,grid_mesh,&time,"./checkpoint.h5");

  //Time-step
  printf("\nParameter");
  printf("\nr_o,r_i=%e,%e",r_o,r_i);
  printf("\nRe_i,Re_o=%e,%e",REYNOLDS_OUTER,REYNOLDS_INNER);
  
  double U_i=REYNOLDS_OUTER*r_o;
  double U_o=REYNOLDS_INNER*r_i;
  
  double t_i=PI2*r_i/NT/U_i;
  double t_o=PI2*r_o/NT/U_o;
 
  printf("\nt_i,t_o=%e,%e",t_o,t_i);
  
  double dt=0.1*min(abs(t_i),abs(t_o));
  
  if(VARIABLE_DT){
    printf("\nRunning with fixed Courant number=%e",COURANT);
  }else{
    printf("\nRunning with fixed Dt=%e",dt);
  }
  

// TayC.h 代码块
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cublas_v2.h>
#include <cufft>
#include <time.h>
#include <sys/time.h>
#include <hdf5.h>
#include <hdf5_hl.h>


#define CHECK_CUDART(x) do { \
  cudaError_t res = (x); \
  if(res != cudaSuccess) { \
    fprintf(stderr, " CUDART: %s = %d (%s) at (%s:%d)\n",  #x, res, cudaGetErrorString(res),__FILE__,__LINE__); \
    exit(1); \
  } \
} while(0)

#define NZ 65 //Number of grid points in real space is 2*(NZ-1)
#define NT 64 //Number of grid points in real space is NT
#define NR 64 //Number of grid points in real space is NR

#define ETA (0.5) //eta=inner_radius/outer_radius, eta<1

#define NZP (3*(NZ-1)/2+1)
#define NTP (3*NT/2)

#define REYNOLDS_INNER (-700.0)
#define REYNOLDS_OUTER ( 700.0)

#define PI2 6.283185307179586
#define LT PI2
#define LZ (PI2/2)

//Number of steps to perform non-linear convolution: save memory
#define NSTEPS_CONV 1 

//Size of the influence matrix
#define MAT 8

//Parameter of the implicit time-integrator
#define C_IMPLICIT (0.51)

//Variable dt 
#define VARIABLE_DT 1
#define MAXDT 0.01
#define COURANT 0.25
#define TOLERANCE_DTERR 5e-5

typedef struct { double2* r;double2* t;double2* z;} vfield;
typedef struct { int Nr; int Nt; int Nz; double Lt ; double Lz;} size_p;

//Globals
extern double2* AUX_T;
extern double Cfl_r, Cfl_z, Cfl_th;

//Finite differences
void fd_weigths(double x0,double* x,int n,double* L1,double*L2);

//Cublas
void setCublas(void);
void transpose_A(double2* u_2,double2* u_1);
void transpose_B(double2* u_2,double2* u_1);
void transpose_infbackward(double* u_2,double* u_1,int mr);
void transpose_infforward(double* u_2,double* u_1,int mr);
void invert_infmatrix(double* src);

//Derivatives compact 3-3
void setDeriv_3_comp(double* mesh,size_p sizes);
void deriv_R_3C(double2* u,double2* aux);
void deriv_RR_3C(double2* u, double2* aux);

//Derivatives explicit 9
void setDeriv_9_exp(double* mesh,size_p sizes);
void deriv_R_9E(double2* u);
void deriv_RR_9E(double2* u);

//Implicit steps explicit 7
void setImplic_7_exp(double* mesh,size_p sizes);
void implicit_step(double2* u, double2* aux1,double* aux2, double* aux3, double* aux4, 
                   double d_implicit,double dt,int flag);
void calcPressure(double2* p, vfield u, double2* dr,
                   double2* aux1,double* aux2, double* aux3, double* aux4);
void deriv_R_7E(double2* u);

//IO
void writeBufferBinary(double* w,const char* file,size_t elsize, size_t elements);
void writeBufferBinary(double* w,const char* file,size_t elsize, size_t elements);
void writeBufferBinaryGPU(double* w,const char* file,size_t elsize,size_t elements);
void writeCheckpoint(vfield u, double* grid_mesh,double* time,char* fileName);
void readCheckpoint(vfield u, double* grid_mesh,double* time,char* fileName);
void writeFieldVis(vfield u, double* grid_mesh,double* time,char* fileName);

//fft
void setFft(size_p sizes);
void fftDestroy(void);
void fftForward(double2* buffer);
void fftBackward(double2* buffer);
void fftForward_reduced(double2* buffer);
void fftBackward_reduced(double2* buffer);


//Utils
void copyVfield(vfield u1, vfield u2);
void onesGPU(double2* u1, size_t elements);
void normalize(double2* u, double norm, size_t elements);
void decoupleForward(vfield u);
void decoupleBackward(vfield u);
void copyBuffer(double2* u1, double2* u2);
void updateNonlinear(vfield n1, vfield n2, double d_implicit);
void calcnorm(double2* u_out, double2* u_in);
void setZeromode(double2* u);
void initField(vfield u, double* mesh);
void writeBufferBinaryGPUreal(double* w,const char* file,
                              size_t elsize,size_t elements);
void zerosGPU(double2* u1, size_t elements);
void diffAbs(double2* u_diff, double2* u);
void fieldAbs(double2* u_diff, double2* u);


//Non-linear
void setNonlinear(double* mesh);
void nonLinear(vfield u, vfield dur,int flag, double time);


//Check routines
void checkDerivatives(void);
void checkNonlinear(void);
void checkImplicit(void);
void checkPressure(void);
void checkLaplacian(void);
void checkTranspose(void);
void checkBoundaries(void);

//Linear 
void setLinear(double* mesh);
void pressureProject(vfield rhs);
void implictStep(vfield u, vfield rhs, double d_implicit, double dt);
void calcLaplacian(vfield rhs, vfield u, double d_implicit, double dt);
void calcPressure_hom(double2* u,double2* aux1,double* aux2, double* aux3, double* aux4
                      ,double2 Bci, double2 Bco);
void implicit_step_hom(double2* u, double2* aux1,double* aux2, double* aux3, 
                        double* aux4, double d_implicit,double dt,int flag,double2 Bci, double2 Bco);
void calcDivergence(double2* div, vfield u);
void calcVorticity(vfield vor, vfield u);
void forcing(double2* uz);

//Boundary conditions
void setBoundary(double* mesh);
void calcHsol(double* M, vfield hsolu_i, vfield hsolu_o, vfield hsolp_i, vfield hsolp_o,
              double d_implicit, double dt);
void imposeBoundarycond(vfield u, vfield hsolu_i, vfield hsolu_o, vfield hsolp_i,
                        vfield hsolp_o, double* M);

//Time integrator
void integrate(vfield u, vfield u_w, vfield rhs, vfield rhs_w,
               int nsteps, double dt,double* time);
void setIntegrator(size_p sizes);
void check_convergence(int* iter, double dt);
void new_step(double* dt);
void measurecorr(vfield up, vfield um);

//Statistics
void setStatistics(double* mesh);
void readTorq(double2* ut,double2* utdr,double time);

//Non-linear
void padForward(double2* aux,double2* u);
void padBackward(double2* u,double2* aux);
void calcDut(vfield u, vfield du);
void calcDuz(vfield u, vfield du);

//Check cublas
void cublasCheck(cublasStatus_t error, const char* function);                                                               

¹. A hybrid MPI-OpenMP parallel implementation for pseudo-spectral simulations with application to Taylor-Couette flow ↩

². The lattice Boltzmann method: Fundamentals and acoustcis ↩

³. https://github.com/demichie/NS_cylindrical/blob/master/NS_.cylindrical.pdf ↩

⁴. ON PRESSURE BOUNDARY CONDITIONS FOR THE INCOMPRESSIBLE NAVIER-STOKES EQUATIONS* ↩

⁵. A Practical Guide to Pseudospectral Methods; Cambridge ↩