座標変換係数を用いて、直角座標以外の計算をやってみよう。

座標変換係数を用いると下のような格子も作れる。

参考：流体計算と差分法、桑原・河村著書

x, y座標系から理想的なξ, η座標系に移して計算を行う。ξ, ηの座標系では、Δξ, Δηは1であり、直四角形である。図3.20のように実際のx, y座標系で、(a)のように作って、座標変換係数を用いることで、計算内では(b)のような格子での計算になる。変換の仕方は簡単である。

2次元座標系の1次偏微分方程式の場合は上の式のように変換してあげるといい。ここで、

とは普通に差分して求めればいい。 Δξ, Δηは1。しかし、は求める必要がある。求め方は別紙で配る。（3次元だと9つ）

2次偏微分方程式の場合以下のように変換するが、

のようにまとめて整理したほうが、後で楽になる。まとめたときの、式は別紙で配る。

ときの求めなければならない変数としては、

の2つである。（3次元の場合は3つ）この変数を求め方も別紙で配る。

注意！！必ず自ら紙でまとめてみること。

後は、Navier-Stokes式を座標変換して、それを今までのように差分で計算すればいい。

2次元の場合、求めないといけない変数6つを先に求める。（nループの前）

以前計算した管内流れのCコードから修正すればいい。以下は修正したコードである。

#include <iostream>

#include <fstream>

#include <cmath>

#include <time.h>

#include <cstdlib>

#include <cstring>

#include <stdio.h>

#include <omp.h>

using namespace std;

#define IMAX 100

#define JMAX 20

#define NMAX 500

#define dt 0.01

#define dx 0.1

#define dy 0.1

#define Re 100.0

int n = 0;

void initialization_values(double u0[][IMAX + 1], double u[][IMAX + 1], double v0[][IMAX + 1], double v[][IMAX + 1], double p0[][IMAX + 1], double rhs[][IMAX + 1],

double x0[][IMAX + 1], double y0[][IMAX + 1], double us[][IMAX + 1], double vs[][IMAX + 1], double i0x[][IMAX + 1], double i0y[][IMAX + 1],

double j0x[][IMAX + 1], double j0y[][IMAX + 1], double d2i[][IMAX + 1], double d2j[][IMAX + 1]) {

int i, j;

#pragma omp parallel for private(i)

for (j = 0; j < JMAX + 1; j++) {

for (i = 0; i < IMAX + 1; i++) {

u0[j][i] = 0.0;

u[j][i] = 0.0;

v0[j][i] = 0.0;

v[j][i] = 0.0;

p0[j][i] = 0.0;

rhs[j][i] = 0.0;

x0[j][i] = 0.0;

y0[j][i] = 0.0;

us[j][i] = 0.0;

vs[j][i] = 0.0;

i0x[j][i] = 0.0;

i0y[j][i] = 0.0;

j0x[j][i] = 0.0;

j0y[j][i] = 0.0;

d2i[j][i] = 0.0;

d2j[j][i] = 0.0;

}

}

}

void domain(double x0[][IMAX + 1], double y0[][IMAX + 1]) {

int i, j;

#pragma omp parallel for private(i)

for (j = 0; j < JMAX + 1; j++) {

for (i = 0; i < IMAX + 1; i++) {

x0[j][i] = dx * (double)i;

y0[j][i] = dy * (double)j;

}

}

}

void change_domain_calculation(double x0[][IMAX + 1], double y0[][IMAX + 1], double i0x[][IMAX + 1], double i0y[][IMAX + 1],

double j0x[][IMAX + 1], double j0y[][IMAX + 1], double d2i[][IMAX + 1], double d2j[][IMAX + 1]) {

double J, kvx, kvy, x0i, y0i, x0j, y0j;

int i, j;

#pragma omp parallel for private(i, J, kvx, kvy, x0i, y0i, x0j, y0j)

for (j = 1; j < JMAX; j++) {

for (i = 1; i < IMAX; i++) {

x0i = (-x0[j][i - 1] + x0[j][i + 1])*0.5;

y0i = (-y0[j][i - 1] + y0[j][i + 1])*0.5;

x0j = (-x0[j - 1][i] + x0[j + 1][i])*0.5;

y0j = (-y0[j - 1][i] + y0[j + 1][i])*0.5;

J = x0i * y0j - x0j * y0i;

i0x[j][i] = y0j / J;

i0y[j][i] = -x0j / J;

j0x[j][i] = -y0i / J;

j0y[j][i] = x0i / J;

kvx = (i0x[j][i] * i0x[j][i] + i0y[j][i] * i0y[j][i])*(x0[j][i - 1] + x0[j][i + 1] - 2.0 * x0[j][i])

+ (j0x[j][i] * j0x[j][i] + j0y[j][i] * j0y[j][i])*(x0[j - 1][i] + x0[j + 1][i] - 2.0 * x0[j][i])

+ 2.0*(i0x[j][i] * j0x[j][i] + i0y[j][i] * j0y[j][i])*(x0[j - 1][i - 1] - x0[j - 1][i + 1] - x0[j + 1][i - 1] + x0[j + 1][i + 1])*0.25;

kvy = (i0x[j][i] * i0x[j][i] + i0y[j][i] * i0y[j][i])*(y0[j][i - 1] + y0[j][i + 1] - 2.0 * y0[j][i])

+ (j0x[j][i] * j0x[j][i] + j0y[j][i] * j0y[j][i])*(y0[j - 1][i] + y0[j + 1][i] - 2.0 * y0[j][i])

+ 2.0*(i0x[j][i] * j0x[j][i] + i0y[j][i] * j0y[j][i])*(y0[j - 1][i - 1] - y0[j - 1][i + 1] - y0[j + 1][i - 1] + y0[j + 1][i + 1])*0.25;

d2i[j][i] = (-kvx * y0j + kvy * x0j) / J;

d2j[j][i] = (kvx * y0i - kvy * x0i) / J;

}

}

}

void initial_condition(double u0[][IMAX + 1], double v0[][IMAX + 1]) {

int i, j;

#pragma omp parallel for private(i)

for (j = 0; j < JMAX + 1; j++) {

for (i = 0; i < IMAX + 1; i++) {

u0[j][i] = 1.0;

v0[j][i] = 0.0;

}

}

}

void vel_boundary_condition(double u0[][IMAX + 1], double v0[][IMAX + 1]) {

int i, j;

#pragma omp parallel for

for (i = 0; i < IMAX + 1; i++) {

u0[0][i] = 0.0;

v0[0][i] = 0.0;

u0[JMAX][i] = 0.0;

v0[JMAX][i] = 0.0;

}

#pragma omp parallel for

for (j = 0; j < JMAX + 1; j++) {

u0[j][0] = 1.0;

v0[j][0] = 0.0;

u0[j][IMAX] = u0[j][IMAX - 1];

v0[j][IMAX] = v0[j][IMAX - 1];

}

}

void pre_boundary_condition(double p0[][IMAX + 1]) {

int i, j;

#pragma omp parallel for

for (i = 0; i < IMAX + 1; i++) {

p0[0][i] = p0[1][i];

p0[JMAX][i] = p0[JMAX - 1][i];

}

#pragma omp parallel for

for (j = 0; j < JMAX + 1; j++) {

p0[j][0] = 0.0;

p0[j][IMAX] = p0[j][IMAX - 1];

}

}

void rhs_calculation(double u0[][IMAX + 1], double v0[][IMAX + 1], double rhs[][IMAX + 1], double i0x[][IMAX + 1], double i0y[][IMAX + 1],

double j0x[][IMAX + 1], double j0y[][IMAX + 1]) {

double umi1, vmi1, upi1, vpi1, umj1, vmj1, upj1, vpj1, ix, iy, jx, jy;

int i, j;

#pragma omp parallel for private(i, umi1, vmi1, upi1, vpi1, umj1, vmj1, upj1, vpj1, ix, iy, jx, jy)

for (j = 1; j < JMAX; j++) {

for (i = 1; i < IMAX; i++) {

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

umi1 = u0[j][i - 1];

vmi1 = v0[j][i - 1];

upi1 = u0[j][i + 1];

vpi1 = v0[j][i + 1];

umj1 = u0[j - 1][i];

vmj1 = v0[j - 1][i];

upj1 = u0[j + 1][i];

vpj1 = v0[j + 1][i];

rhs[j][i] = ((-umi1 + upi1)*ix + (-umj1 + upj1)*jx + (-vmi1 + vpi1)*iy + (-vmj1 + vpj1)*jy)*0.5 / dt

- (((-umi1 + upi1)*ix + (-umj1 + upj1)*jx) * ((-umi1 + upi1)*ix + (-umj1 + upj1)*jx))*0.25

- (((-vmi1 + vpi1)*iy + (-vmj1 + vpj1)*jy) * ((-vmi1 + vpi1)*iy + (-vmj1 + vpj1)*jy))*0.25

- 0.5 * (((-vmi1 + vpi1)*ix + (-vmj1 + vpj1)*jx) * ((-umi1 + upi1)*iy + (-umj1 + upj1)*jy));

}

}

}

void pre_calculation(double p0[][IMAX + 1], double rhs[][IMAX + 1], double i0x[][IMAX + 1], double i0y[][IMAX + 1],

double j0x[][IMAX + 1], double j0y[][IMAX + 1], double d2i[][IMAX + 1], double d2j[][IMAX + 1]) {

double p000, pmi1, ppi1, pmj1, ppj1, pmm1, ppp1, pmp1, ppm1, ix, iy, jx, jy, di, dj, SGM = 1.05;

int i, j;

#pragma omp parallel for private(i,p000, pmi1, ppi1, pmj1, ppj1, pmm1, ppp1, pmp1, ppm1, ix, iy, jx, jy, di, dj) firstprivate(SGM)

for (j = 1; j < JMAX; j += 2) {

for (i = 1; i < IMAX; i++) {

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

p000 = p0[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

pmm1 = p0[j - 1][i - 1];

ppp1 = p0[j + 1][i + 1];

pmp1 = p0[j - 1][i + 1];

ppm1 = p0[j + 1][i - 1];

p0[j][i] = (

(ix * ix + iy * iy)*(pmi1 + ppi1) + (jx * jx + jy * jy)*(pmj1 + ppj1)

+ 2.0*(ix * jx + iy * jy)*(pmm1 - pmp1 - ppm1 + ppp1)

+ di * (-pmi1 + ppi1)*0.5 + dj * (-pmj1 + ppj1)*0.5

- rhs[j][i]) / (2.0*(ix * ix + iy * iy + jx * jx + jy * jy));

p0[j][i] = p000 + SGM * (p0[j][i] - p000);

}

}

#pragma omp parallel for private(i,p000, pmi1, ppi1, pmj1, ppj1, pmm1, ppp1, pmp1, ppm1, ix, iy, jx, jy, di, dj) firstprivate(SGM)

for (j = 2; j < JMAX; j += 2) {

for (i = 1; i < IMAX; i++) {

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

p000 = p0[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

pmm1 = p0[j - 1][i - 1];

ppp1 = p0[j + 1][i + 1];

pmp1 = p0[j - 1][i + 1];

ppm1 = p0[j + 1][i - 1];

p0[j][i] = (

(ix * ix + iy * iy)*(pmi1 + ppi1) + (jx * jx + jy * jy)*(pmj1 + ppj1)

+ 2.0*(ix * jx + iy * jy)*(pmm1 - pmp1 - ppm1 + ppp1)

+ di * (-pmi1 + ppi1)*0.5 + dj * (-pmj1 + ppj1)*0.5

- rhs[j][i]) / (2.0*(ix * ix + iy * iy + jx * jx + jy * jy));

p0[j][i] = p000 + SGM * (p0[j][i] - p000);

}

}

}

void vel_explicit_calculation(double u0[][IMAX + 1], double u[][IMAX + 1], double v0[][IMAX + 1], double v[][IMAX + 1], double p0[][IMAX + 1], double us[][IMAX + 1], double vs[][IMAX + 1],

double i0x[][IMAX + 1], double i0y[][IMAX + 1], double j0x[][IMAX + 1], double j0y[][IMAX + 1], double d2i[][IMAX + 1], double d2j[][IMAX + 1]) {

double u000, umi1, upi1, umj1, upj1, umi2, upi2, umj2, upj2,

v000, vmi1, vpi1, vmj1, vpj1, vmi2, vpi2, vmj2, vpj2,

rhs_u, rhs_v, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1,

pmi1, ppi1, pmj1, ppj1;

int i, j;

#pragma omp parallel for private(i, u000, umi1, upi1, umj1, upj1, umi2, upi2, umj2, upj2, v000, vmi1, vpi1, vmj1, vpj1, vmi2, vpi2, vmj2, vpj2, rhs_u, rhs_v, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1)

for (j = 2; j < JMAX - 1; j++) {

for (i = 2; i < IMAX - 1; i++) {

u000 = u[j][i];

v000 = v[j][i];

umi1 = u[j][i - 1];

vmi1 = v[j][i - 1];

upi1 = u[j][i + 1];

vpi1 = v[j][i + 1];

umj1 = u[j - 1][i];

vmj1 = v[j - 1][i];

upj1 = u[j + 1][i];

vpj1 = v[j + 1][i];

umi2 = u[j][i - 2];

vmi2 = v[j][i - 2];

upi2 = u[j][i + 2];

vpi2 = v[j][i + 2];

umj2 = u[j - 2][i];

vmj2 = v[j - 2][i];

upj2 = u[j + 2][i];

vpj2 = v[j + 2][i];

umm1 = u[j - 1][i - 1];

upp1 = u[j + 1][i + 1];

ump1 = u[j - 1][i + 1];

upm1 = u[j + 1][i - 1];

vmm1 = v[j - 1][i - 1];

vpp1 = v[j + 1][i + 1];

vmp1 = v[j - 1][i + 1];

vpm1 = v[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1 - 2.0*u000) + (jx * jx + jy * jy)*(umj1 + upj1 - 2.0*u000)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(umi2 + 8.0*(-umi1 + upi1) - upi2) + fabs(u000*ix + v000 * iy)*(umi2 - 4.0*(umi1 + upi1) + upi2 + 6.0*u000)) / 12.0

+ ((u000*jx + v000 * jy)*(umj2 + 8.0*(-umj1 + upj1) - upj2) + fabs(u000*jx + v000 * jy)*(umj2 - 4.0*(umj1 + upj1) + upj2 + 6.0*u000)) / 12.0

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1 - 2.0*v000) + (jx * jx + jy * jy)*(vmj1 + vpj1 - 2.0*v000)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(vmi2 + 8.0*(-vmi1 + vpi1) - vpi2) + fabs(u000*ix + v000 * iy)*(vmi2 - 4.0*(vmi1 + vpi1) + vpi2 + 6.0*v000)) / 12.0

+ ((u000*jx + v000 * jy)*(vmj2 + 8.0*(-vmj1 + vpj1) - vpj2) + fabs(u000*jx + v000 * jy)*(vmj2 - 4.0*(vmj1 + vpj1) + vpj2 + 6.0*v000)) / 12.0

);

u0[j][i] = u000 + dt * 0.5*(3.0*rhs_u - us[j][i]);

v0[j][i] = v000 + dt * 0.5*(3.0*rhs_v - vs[j][i]);

us[j][i] = rhs_u;

vs[j][i] = rhs_v;

}

}

for (j = 1; j < JMAX; j += JMAX - 2) {

#pragma omp parallel for private(u000, umi1, upi1, umj1, upj1, v000, vmi1, vpi1, vmj1, vpj1, rhs_u, rhs_v, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1)

for (i = 1; i < IMAX; i++) {

u000 = u[j][i];

v000 = v[j][i];

umi1 = u[j][i - 1];

vmi1 = v[j][i - 1];

upi1 = u[j][i + 1];

vpi1 = v[j][i + 1];

umj1 = u[j - 1][i];

vmj1 = v[j - 1][i];

upj1 = u[j + 1][i];

vpj1 = v[j + 1][i];

umm1 = u[j - 1][i - 1];

upp1 = u[j + 1][i + 1];

ump1 = u[j - 1][i + 1];

upm1 = u[j + 1][i - 1];

vmm1 = v[j - 1][i - 1];

vpp1 = v[j + 1][i + 1];

vmp1 = v[j - 1][i + 1];

vpm1 = v[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1 - 2.0*u000) + (jx * jx + jy * jy)*(umj1 + upj1 - 2.0*u000)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(-umi1 + upi1) + fabs(u000*ix + v000 * iy)*(-umi1 - upi1 + 2.0*u000)) *0.5

+ ((u000*jx + v000 * jy)*(-umj1 + upj1) + fabs(u000*jx + v000 * jy)*(-umj1 - upj1 + 2.0*u000)) *0.5

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1 - 2.0*v000) + (jx * jx + jy * jy)*(vmj1 + vpj1 - 2.0*v000)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(-vmi1 + vpi1) + fabs(u000*ix + v000 * iy)*(-vmi1 - vpi1 + 2.0*v000)) *0.5

+ ((u000*jx + v000 * jy)*(-vmj1 + vpj1) + fabs(u000*jx + v000 * jy)*(-vmj1 - vpj1 + 2.0*v000)) *0.5

);

u0[j][i] = u000 + dt * 0.5*(3.0*rhs_u - us[j][i]);

v0[j][i] = v000 + dt * 0.5*(3.0*rhs_v - vs[j][i]);

us[j][i] = rhs_u;

vs[j][i] = rhs_v;

}

}

#pragma omp parallel for private(i, u000, umi1, upi1, umj1, upj1, v000, vmi1, vpi1, vmj1, vpj1, rhs_u, rhs_v, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1)

for (j = 2; j < JMAX - 1; j++) {

for (i = 1; i < IMAX; i += IMAX - 2) {

u000 = u[j][i];

v000 = v[j][i];

umi1 = u[j][i - 1];

vmi1 = v[j][i - 1];

upi1 = u[j][i + 1];

vpi1 = v[j][i + 1];

umj1 = u[j - 1][i];

vmj1 = v[j - 1][i];

upj1 = u[j + 1][i];

vpj1 = v[j + 1][i];

umm1 = u[j - 1][i - 1];

upp1 = u[j + 1][i + 1];

ump1 = u[j - 1][i + 1];

upm1 = u[j + 1][i - 1];

vmm1 = v[j - 1][i - 1];

vpp1 = v[j + 1][i + 1];

vmp1 = v[j - 1][i + 1];

vpm1 = v[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1 - 2.0*u000) + (jx * jx + jy * jy)*(umj1 + upj1 - 2.0*u000)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(-umi1 + upi1) + fabs(u000*ix + v000 * iy)*(-umi1 - upi1 + 2.0*u000)) *0.5

+ ((u000*jx + v000 * jy)*(-umj1 + upj1) + fabs(u000*jx + v000 * jy)*(-umj1 - upj1 + 2.0*u000)) *0.5

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1 - 2.0*v000) + (jx * jx + jy * jy)*(vmj1 + vpj1 - 2.0*v000)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(-vmi1 + vpi1) + fabs(u000*ix + v000 * iy)*(-vmi1 - vpi1 + 2.0*v000)) *0.5

+ ((u000*jx + v000 * jy)*(-vmj1 + vpj1) + fabs(u000*jx + v000 * jy)*(-vmj1 - vpj1 + 2.0*v000)) *0.5

);

u0[j][i] = u000 + dt * 0.5*(3.0*rhs_u - us[j][i]);

v0[j][i] = v000 + dt * 0.5*(3.0*rhs_v - vs[j][i]);

us[j][i] = rhs_u;

vs[j][i] = rhs_v;

}

}

}

void vel_implicit_calculation(double u0[][IMAX + 1], double u[][IMAX + 1], double v0[][IMAX + 1], double v[][IMAX + 1], double p0[][IMAX + 1], double us[][IMAX + 1], double vs[][IMAX + 1],

double i0x[][IMAX + 1], double i0y[][IMAX + 1], double j0x[][IMAX + 1], double j0y[][IMAX + 1], double d2i[][IMAX + 1], double d2j[][IMAX + 1]) {

double u000, umi1, upi1, umj1, upj1, umi2, upi2, umj2, upj2,

v000, vmi1, vpi1, vmj1, vpj1, vmi2, vpi2, vmj2, vpj2,

rhs_u, rhs_v, leftvel, SGM2 = 1.05, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1,

pmi1, ppi1, pmj1, ppj1;

int i, j;

#pragma omp parallel for private(i, u000, umi1, upi1, umj1, upj1, umi2, upi2, umj2, upj2, v000, vmi1, vpi1, vmj1, vpj1, vmi2, vpi2, vmj2, vpj2, rhs_u, rhs_v, leftvel, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1) firstprivate(SGM2)

for (j = 2; j < JMAX - 1; j += 3) {

for (i = 2; i < IMAX - 1; i++) {

u000 = u0[j][i];

v000 = v0[j][i];

umi1 = u0[j][i - 1];

vmi1 = v0[j][i - 1];

upi1 = u0[j][i + 1];

vpi1 = v0[j][i + 1];

umj1 = u0[j - 1][i];

vmj1 = v0[j - 1][i];

upj1 = u0[j + 1][i];

vpj1 = v0[j + 1][i];

umi2 = u0[j][i - 2];

vmi2 = v0[j][i - 2];

upi2 = u0[j][i + 2];

vpi2 = v0[j][i + 2];

umj2 = u0[j - 2][i];

vmj2 = v0[j - 2][i];

upj2 = u0[j + 2][i];

vpj2 = v0[j + 2][i];

umm1 = u0[j - 1][i - 1];

upp1 = u0[j + 1][i + 1];

ump1 = u0[j - 1][i + 1];

upm1 = u0[j + 1][i - 1];

vmm1 = v0[j - 1][i - 1];

vpp1 = v0[j + 1][i + 1];

vmp1 = v0[j - 1][i + 1];

vpm1 = v0[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

leftvel = 2.0*(ix * ix + iy * iy + jx * jx + jy * jy) / Re + (fabs(u000*ix + v000 * iy) + fabs(u000*jx + v000 * jy))*0.5;

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1) + (jx * jx + jy * jy)*(umj1 + upj1)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(umi2 + 8.0*(-umi1 + upi1) - upi2) + fabs(u000*ix + v000 * iy)*(umi2 - 4.0*(umi1 + upi1) + upi2)) / 12.0

+ ((u000*jx + v000 * jy)*(umj2 + 8.0*(-umj1 + upj1) - upj2) + fabs(u000*jx + v000 * jy)*(umj2 - 4.0*(umj1 + upj1) + upj2)) / 12.0

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1) + (jx * jx + jy * jy)*(vmj1 + vpj1)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(vmi2 + 8.0*(-vmi1 + vpi1) - vpi2) + fabs(u000*ix + v000 * iy)*(vmi2 - 4.0*(vmi1 + vpi1) + vpi2)) / 12.0

+ ((u000*jx + v000 * jy)*(vmj2 + 8.0*(-vmj1 + vpj1) - vpj2) + fabs(u000*jx + v000 * jy)*(vmj2 - 4.0*(vmj1 + vpj1) + vpj2)) / 12.0

);

u0[j][i] = (u[j][i] + dt * 0.5 * (rhs_u + us[j][i])) / (1.0 + dt * 0.5 * leftvel);

v0[j][i] = (v[j][i] + dt * 0.5 * (rhs_v + vs[j][i])) / (1.0 + dt * 0.5 * leftvel);

u0[j][i] = u000 + SGM2 * (u0[j][i] - u000);

v0[j][i] = v000 + SGM2 * (v0[j][i] - v000);

}

}

#pragma omp parallel for private(i, u000, umi1, upi1, umj1, upj1, umi2, upi2, umj2, upj2, v000, vmi1, vpi1, vmj1, vpj1, vmi2, vpi2, vmj2, vpj2, rhs_u, rhs_v, leftvel, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1) firstprivate(SGM2)

for (j = 3; j < JMAX - 1; j += 3) {

for (i = 2; i < IMAX - 1; i++) {

u000 = u0[j][i];

v000 = v0[j][i];

umi1 = u0[j][i - 1];

vmi1 = v0[j][i - 1];

upi1 = u0[j][i + 1];

vpi1 = v0[j][i + 1];

umj1 = u0[j - 1][i];

vmj1 = v0[j - 1][i];

upj1 = u0[j + 1][i];

vpj1 = v0[j + 1][i];

umi2 = u0[j][i - 2];

vmi2 = v0[j][i - 2];

upi2 = u0[j][i + 2];

vpi2 = v0[j][i + 2];

umj2 = u0[j - 2][i];

vmj2 = v0[j - 2][i];

upj2 = u0[j + 2][i];

vpj2 = v0[j + 2][i];

umm1 = u0[j - 1][i - 1];

upp1 = u0[j + 1][i + 1];

ump1 = u0[j - 1][i + 1];

upm1 = u0[j + 1][i - 1];

vmm1 = v0[j - 1][i - 1];

vpp1 = v0[j + 1][i + 1];

vmp1 = v0[j - 1][i + 1];

vpm1 = v0[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

leftvel = 2.0*(ix * ix + iy * iy + jx * jx + jy * jy) / Re + (fabs(u000*ix + v000 * iy) + fabs(u000*jx + v000 * jy))*0.5;

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1) + (jx * jx + jy * jy)*(umj1 + upj1)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(umi2 + 8.0*(-umi1 + upi1) - upi2) + fabs(u000*ix + v000 * iy)*(umi2 - 4.0*(umi1 + upi1) + upi2)) / 12.0

+ ((u000*jx + v000 * jy)*(umj2 + 8.0*(-umj1 + upj1) - upj2) + fabs(u000*jx + v000 * jy)*(umj2 - 4.0*(umj1 + upj1) + upj2)) / 12.0

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1) + (jx * jx + jy * jy)*(vmj1 + vpj1)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(vmi2 + 8.0*(-vmi1 + vpi1) - vpi2) + fabs(u000*ix + v000 * iy)*(vmi2 - 4.0*(vmi1 + vpi1) + vpi2)) / 12.0

+ ((u000*jx + v000 * jy)*(vmj2 + 8.0*(-vmj1 + vpj1) - vpj2) + fabs(u000*jx + v000 * jy)*(vmj2 - 4.0*(vmj1 + vpj1) + vpj2)) / 12.0

);

u0[j][i] = (u[j][i] + dt * 0.5 * (rhs_u + us[j][i])) / (1.0 + dt * 0.5 * leftvel);

v0[j][i] = (v[j][i] + dt * 0.5 * (rhs_v + vs[j][i])) / (1.0 + dt * 0.5 * leftvel);

u0[j][i] = u000 + SGM2 * (u0[j][i] - u000);

v0[j][i] = v000 + SGM2 * (v0[j][i] - v000);

}

}

#pragma omp parallel for private(i, u000, umi1, upi1, umj1, upj1, umi2, upi2, umj2, upj2, v000, vmi1, vpi1, vmj1, vpj1, vmi2, vpi2, vmj2, vpj2, rhs_u, rhs_v, leftvel, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1) firstprivate(SGM2)

for (j = 4; j < JMAX - 1; j += 3) {

for (i = 2; i < IMAX - 1; i++) {

u000 = u0[j][i];

v000 = v0[j][i];

umi1 = u0[j][i - 1];

vmi1 = v0[j][i - 1];

upi1 = u0[j][i + 1];

vpi1 = v0[j][i + 1];

umj1 = u0[j - 1][i];

vmj1 = v0[j - 1][i];

upj1 = u0[j + 1][i];

vpj1 = v0[j + 1][i];

umi2 = u0[j][i - 2];

vmi2 = v0[j][i - 2];

upi2 = u0[j][i + 2];

vpi2 = v0[j][i + 2];

umj2 = u0[j - 2][i];

vmj2 = v0[j - 2][i];

upj2 = u0[j + 2][i];

vpj2 = v0[j + 2][i];

umm1 = u0[j - 1][i - 1];

upp1 = u0[j + 1][i + 1];

ump1 = u0[j - 1][i + 1];

upm1 = u0[j + 1][i - 1];

vmm1 = v0[j - 1][i - 1];

vpp1 = v0[j + 1][i + 1];

vmp1 = v0[j - 1][i + 1];

vpm1 = v0[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

leftvel = 2.0*(ix * ix + iy * iy + jx * jx + jy * jy) / Re + (fabs(u000*ix + v000 * iy) + fabs(u000*jx + v000 * jy))*0.5;

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1) + (jx * jx + jy * jy)*(umj1 + upj1)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(umi2 + 8.0*(-umi1 + upi1) - upi2) + fabs(u000*ix + v000 * iy)*(umi2 - 4.0*(umi1 + upi1) + upi2)) / 12.0

+ ((u000*jx + v000 * jy)*(umj2 + 8.0*(-umj1 + upj1) - upj2) + fabs(u000*jx + v000 * jy)*(umj2 - 4.0*(umj1 + upj1) + upj2)) / 12.0

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1) + (jx * jx + jy * jy)*(vmj1 + vpj1)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(vmi2 + 8.0*(-vmi1 + vpi1) - vpi2) + fabs(u000*ix + v000 * iy)*(vmi2 - 4.0*(vmi1 + vpi1) + vpi2)) / 12.0

+ ((u000*jx + v000 * jy)*(vmj2 + 8.0*(-vmj1 + vpj1) - vpj2) + fabs(u000*jx + v000 * jy)*(vmj2 - 4.0*(vmj1 + vpj1) + vpj2)) / 12.0

);

u0[j][i] = (u[j][i] + dt * 0.5 * (rhs_u + us[j][i])) / (1.0 + dt * 0.5 * leftvel);

v0[j][i] = (v[j][i] + dt * 0.5 * (rhs_v + vs[j][i])) / (1.0 + dt * 0.5 * leftvel);

u0[j][i] = u000 + SGM2 * (u0[j][i] - u000);

v0[j][i] = v000 + SGM2 * (v0[j][i] - v000);

}

}

for (j = 1; j < JMAX; j += JMAX - 2) {

#pragma omp parallel for private(u000, umi1, upi1, umj1, upj1, v000, vmi1, vpi1, vmj1, vpj1, rhs_u, rhs_v, leftvel, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1) firstprivate(SGM2)

for (i = 1; i < IMAX; i += 2) {

u000 = u0[j][i];

v000 = v0[j][i];

umi1 = u0[j][i - 1];

vmi1 = v0[j][i - 1];

upi1 = u0[j][i + 1];

vpi1 = v0[j][i + 1];

umj1 = u0[j - 1][i];

vmj1 = v0[j - 1][i];

upj1 = u0[j + 1][i];

vpj1 = v0[j + 1][i];

umm1 = u0[j - 1][i - 1];

upp1 = u0[j + 1][i + 1];

ump1 = u0[j - 1][i + 1];

upm1 = u0[j + 1][i - 1];

vmm1 = v0[j - 1][i - 1];

vpp1 = v0[j + 1][i + 1];

vmp1 = v0[j - 1][i + 1];

vpm1 = v0[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

leftvel = 2.0*(ix * ix + iy * iy + jx * jx + jy * jy) / Re + fabs(u000*ix + v000 * iy) + fabs(u000*jx + v000 * jy);

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1) + (jx * jx + jy * jy)*(umj1 + upj1)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(-umi1 + upi1) + fabs(u000*ix + v000 * iy)*(-umi1 - upi1)) *0.5

+ ((u000*jx + v000 * jy)*(-umj1 + upj1) + fabs(u000*jx + v000 * jy)*(-umj1 - upj1)) *0.5

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1) + (jx * jx + jy * jy)*(vmj1 + vpj1)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(-vmi1 + vpi1) + fabs(u000*ix + v000 * iy)*(-vmi1 - vpi1)) *0.5

+ ((u000*jx + v000 * jy)*(-vmj1 + vpj1) + fabs(u000*jx + v000 * jy)*(-vmj1 - vpj1)) *0.5

);

u0[j][i] = (u[j][i] + dt * 0.5 * (rhs_u + us[j][i])) / (1.0 + dt * 0.5 * leftvel);

v0[j][i] = (v[j][i] + dt * 0.5 * (rhs_v + vs[j][i])) / (1.0 + dt * 0.5 * leftvel);

u0[j][i] = u000 + SGM2 * (u0[j][i] - u000);

v0[j][i] = v000 + SGM2 * (v0[j][i] - v000);

}

#pragma omp parallel for private(u000, umi1, upi1, umj1, upj1, v000, vmi1, vpi1, vmj1, vpj1, rhs_u, rhs_v, leftvel, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1) firstprivate(SGM2)

for (i = 2; i < IMAX; i += 2) {

u000 = u0[j][i];

v000 = v0[j][i];

umi1 = u0[j][i - 1];

vmi1 = v0[j][i - 1];

upi1 = u0[j][i + 1];

vpi1 = v0[j][i + 1];

umj1 = u0[j - 1][i];

vmj1 = v0[j - 1][i];

upj1 = u0[j + 1][i];

vpj1 = v0[j + 1][i];

umm1 = u0[j - 1][i - 1];

upp1 = u0[j + 1][i + 1];

ump1 = u0[j - 1][i + 1];

upm1 = u0[j + 1][i - 1];

vmm1 = v0[j - 1][i - 1];

vpp1 = v0[j + 1][i + 1];

vmp1 = v0[j - 1][i + 1];

vpm1 = v0[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

leftvel = 2.0*(ix * ix + iy * iy + jx * jx + jy * jy) / Re + fabs(u000*ix + v000 * iy) + fabs(u000*jx + v000 * jy);

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1) + (jx * jx + jy * jy)*(umj1 + upj1)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(-umi1 + upi1) + fabs(u000*ix + v000 * iy)*(-umi1 - upi1)) *0.5

+ ((u000*jx + v000 * jy)*(-umj1 + upj1) + fabs(u000*jx + v000 * jy)*(-umj1 - upj1)) *0.5

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1) + (jx * jx + jy * jy)*(vmj1 + vpj1)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(-vmi1 + vpi1) + fabs(u000*ix + v000 * iy)*(-vmi1 - vpi1)) *0.5

+ ((u000*jx + v000 * jy)*(-vmj1 + vpj1) + fabs(u000*jx + v000 * jy)*(-vmj1 - vpj1)) *0.5

);

u0[j][i] = (u[j][i] + dt * 0.5 * (rhs_u + us[j][i])) / (1.0 + dt * 0.5 * leftvel);

v0[j][i] = (v[j][i] + dt * 0.5 * (rhs_v + vs[j][i])) / (1.0 + dt * 0.5 * leftvel);

u0[j][i] = u000 + SGM2 * (u0[j][i] - u000);

v0[j][i] = v000 + SGM2 * (v0[j][i] - v000);

}

}

#pragma omp parallel for private(i, u000, umi1, upi1, umj1, upj1, v000, vmi1, vpi1, vmj1, vpj1, rhs_u, rhs_v, leftvel, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1) firstprivate(SGM2)

for (j = 2; j < JMAX - 1; j += 2) {

for (i = 1; i < IMAX; i += IMAX - 2) {

u000 = u0[j][i];

v000 = v0[j][i];

umi1 = u0[j][i - 1];

vmi1 = v0[j][i - 1];

upi1 = u0[j][i + 1];

vpi1 = v0[j][i + 1];

umj1 = u0[j - 1][i];

vmj1 = v0[j - 1][i];

upj1 = u0[j + 1][i];

vpj1 = v0[j + 1][i];

umm1 = u0[j - 1][i - 1];

upp1 = u0[j + 1][i + 1];

ump1 = u0[j - 1][i + 1];

upm1 = u0[j + 1][i - 1];

vmm1 = v0[j - 1][i - 1];

vpp1 = v0[j + 1][i + 1];

vmp1 = v0[j - 1][i + 1];

vpm1 = v0[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

leftvel = 2.0*(ix * ix + iy * iy + jx * jx + jy * jy) / Re + fabs(u000*ix + v000 * iy) + fabs(u000*jx + v000 * jy);

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1) + (jx * jx + jy * jy)*(umj1 + upj1)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(-umi1 + upi1) + fabs(u000*ix + v000 * iy)*(-umi1 - upi1)) *0.5

+ ((u000*jx + v000 * jy)*(-umj1 + upj1) + fabs(u000*jx + v000 * jy)*(-umj1 - upj1)) *0.5

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1) + (jx * jx + jy * jy)*(vmj1 + vpj1)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(-vmi1 + vpi1) + fabs(u000*ix + v000 * iy)*(-vmi1 - vpi1)) *0.5

+ ((u000*jx + v000 * jy)*(-vmj1 + vpj1) + fabs(u000*jx + v000 * jy)*(-vmj1 - vpj1)) *0.5

);

u0[j][i] = (u[j][i] + dt * 0.5 * (rhs_u + us[j][i])) / (1.0 + dt * 0.5 * leftvel);

v0[j][i] = (v[j][i] + dt * 0.5 * (rhs_v + vs[j][i])) / (1.0 + dt * 0.5 * leftvel);

u0[j][i] = u000 + SGM2 * (u0[j][i] - u000);

v0[j][i] = v000 + SGM2 * (v0[j][i] - v000);

}

}

#pragma omp parallel for private(i, u000, umi1, upi1, umj1, upj1, v000, vmi1, vpi1, vmj1, vpj1, rhs_u, rhs_v, leftvel, ix, iy, jx, jy, di, dj, umm1, upp1, ump1, upm1, vmm1, vpp1, vmp1, vpm1, pmi1, ppi1, pmj1, ppj1) firstprivate(SGM2)

for (j = 3; j < JMAX - 1; j += 2) {

for (i = 1; i < IMAX; i += IMAX - 2) {

u000 = u0[j][i];

v000 = v0[j][i];

umi1 = u0[j][i - 1];

vmi1 = v0[j][i - 1];

upi1 = u0[j][i + 1];

vpi1 = v0[j][i + 1];

umj1 = u0[j - 1][i];

vmj1 = v0[j - 1][i];

upj1 = u0[j + 1][i];

vpj1 = v0[j + 1][i];

umm1 = u0[j - 1][i - 1];

upp1 = u0[j + 1][i + 1];

ump1 = u0[j - 1][i + 1];

upm1 = u0[j + 1][i - 1];

vmm1 = v0[j - 1][i - 1];

vpp1 = v0[j + 1][i + 1];

vmp1 = v0[j - 1][i + 1];

vpm1 = v0[j + 1][i - 1];

ix = i0x[j][i];

jx = j0x[j][i];

iy = i0y[j][i];

jy = j0y[j][i];

di = d2i[j][i];

dj = d2j[j][i];

pmi1 = p0[j][i - 1];

ppi1 = p0[j][i + 1];

pmj1 = p0[j - 1][i];

ppj1 = p0[j + 1][i];

leftvel = 2.0*(ix * ix + iy * iy + jx * jx + jy * jy) / Re + fabs(u000*ix + v000 * iy) + fabs(u000*jx + v000 * jy);

rhs_u = ((ix * ix + iy * iy)*(umi1 + upi1) + (jx * jx + jy * jy)*(umj1 + upj1)

+ 2.0*(ix * jx + iy * jy)*(umm1 - ump1 - upm1 + upp1)

+ di * (-umi1 + upi1)*0.5 + dj * (-umj1 + upj1)*0.5) / Re

- (

((-pmi1 + ppi1) *ix + (-pmj1 + ppj1)*jx)*0.5

+ ((u000*ix + v000 * iy)*(-umi1 + upi1) + fabs(u000*ix + v000 * iy)*(-umi1 - upi1)) *0.5

+ ((u000*jx + v000 * jy)*(-umj1 + upj1) + fabs(u000*jx + v000 * jy)*(-umj1 - upj1)) *0.5

);

rhs_v = ((ix * ix + iy * iy)*(vmi1 + vpi1) + (jx * jx + jy * jy)*(vmj1 + vpj1)

+ 2.0*(ix * jx + iy * jy)*(vmm1 - vmp1 - vpm1 + vpp1)

+ di * (-vmi1 + vpi1)*0.5 + dj * (-vmj1 + vpj1)*0.5) / Re

- (

((-pmi1 + ppi1) *iy + (-pmj1 + ppj1)*jy)*0.5

+ ((u000*ix + v000 * iy)*(-vmi1 + vpi1) + fabs(u000*ix + v000 * iy)*(-vmi1 - vpi1)) *0.5

+ ((u000*jx + v000 * jy)*(-vmj1 + vpj1) + fabs(u000*jx + v000 * jy)*(-vmj1 - vpj1)) *0.5

);

u0[j][i] = (u[j][i] + dt * 0.5 * (rhs_u + us[j][i])) / (1.0 + dt * 0.5 * leftvel);

v0[j][i] = (v[j][i] + dt * 0.5 * (rhs_v + vs[j][i])) / (1.0 + dt * 0.5 * leftvel);

u0[j][i] = u000 + SGM2 * (u0[j][i] - u000);

v0[j][i] = v000 + SGM2 * (v0[j][i] - v000);

}

}

}

void write_data_vt(double u0[][IMAX + 1], double x[][IMAX + 1], double y[][IMAX + 1]) {

ofstream rdi;

rdi.precision(18);

char str[20];

sprintf_s(str, "data%d.dat", n);// Convert i to a char

string datafile = "./data/";

datafile.append(str);

rdi.open(datafile.c_str());

for (int j = 0; j < JMAX + 1; j++) {

for (int i = 0; i < IMAX + 1; i++) {

rdi << x[j][i] << ' ' << y[j][i] << ' ' << u0[j][i] << endl;

}

rdi << endl;

}

rdi.close();

}

int main() {

cout.precision(15);

double check0 = 0.0, check1 = 0.0, check2 = 0.0, check3 = 0.0;

typedef double IJ[IMAX + 1];

IJ *u0, *u, *v0, *v, *p0, *rhs, *x0, *y0, *us, *vs, *i0x, *i0y, *j0x, *j0y, *d2i, *d2j;

u0 = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

u = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

v0 = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

v = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

p0 = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

rhs = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

x0 = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

y0 = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

us = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

vs = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

i0x = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

i0y = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

j0x = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

j0y = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

d2i = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

d2j = (IJ *)malloc((IMAX + 1) * (JMAX + 1) * sizeof(double));

initialization_values(u0, u, v0, v, p0, rhs, x0, y0, us, vs, i0x, i0y, j0x, j0y, d2i, d2j);

domain(x0, y0);

change_domain_calculation(x0, y0, i0x, i0y, j0x, j0y, d2i, d2j);

initial_condition(u0, v0);

vel_boundary_condition(u0, v0);

pre_boundary_condition(p0);

write_data_vt(u0, x0, y0);

for (n = 1; n < NMAX + 1; n++) {

rhs_calculation(u0, v0, rhs, i0x, i0y, j0x, j0y);

for (int loop = 1; loop < 51; loop++) {

pre_calculation(p0, rhs, i0x, i0y, j0x, j0y, d2i, d2j);

}

check0 = p0[5][5];

pre_calculation(p0, rhs, i0x, i0y, j0x, j0y, d2i, d2j);

check1 = p0[5][5];

pre_boundary_condition(p0);

memcpy(u, u0, ((IMAX + 1)*(JMAX + 1) * sizeof(double)));

memcpy(v, v0, ((IMAX + 1)*(JMAX + 1) * sizeof(double)));

vel_explicit_calculation(u0, u, v0, v, p0, us, vs, i0x, i0y, j0x, j0y, d2i, d2j);

for (int loop = 1; loop < 6; loop++) {

vel_implicit_calculation(u0, u, v0, v, p0, us, vs, i0x, i0y, j0x, j0y, d2i, d2j);

}

check2 = u0[5][5];

vel_implicit_calculation(u0, u, v0, v, p0, us, vs, i0x, i0y, j0x, j0y, d2i, d2j);

check3 = u0[5][5];

if (n % 100 == 0) {

cout << n << ' ' << fabs(check1 - check0) << ' ' << fabs(check2 - check3) << endl;

}

vel_boundary_condition(u0, v0);

if (n % 100 == 0) {

write_data_vt(u0, x0, y0);

}

}

free(u0);

free(u);

free(v0);

free(v);

free(p0);

free(rhs);

free(x0);

free(y0);

free(us);

free(vs);

free(i0x);

free(i0y);

free(j0x);

free(j0y);

free(d2i);

free(d2j);

return 0;

}

実際座標変換係数を入れて計算をしてみて、以前と同じ結果が出ればいい。修正するときは、rhs_calculationの部分を先に修正して、そのまま計算を回して、正しい結果が出るか確認する。できたら、pre_calculationの部分を修正し確認し、正しくなっているかを確認する。その後は陽解法の速度、最終的には陰解法の速度を修正し確認すればいい。

しかし、今の格子だと格子変換係数を使う意味がない。そのため、y方向格子間隔を非等間隔に作って計算をしてみよう。

現在のy方向格子の式は

y0[j][i] = dy * (double)j;

である。

これを、等比数列の和の式を使って、yが0付近の格子を細かくし、2に近くなるほど粗くなる格子を作ってみよう。

等比数列の和の式は、

ここで、aは初項、rは公比である。実際、y方向の格子を以下のように修正し計算する。

y0[j][i] = 0.05*(pow(1.0677448,(double)(j))-1.0)/(0.0677448);

等間隔の格子と非等間隔の格子の値を比較すると以下のような表になる。

j	0	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20
等間隔	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	1.1	1.2	1.3	1.4	1.5	1.6	1.7	1.8	1.9	2
非等間隔	0	0.05	0.103387	0.160391	0.221257	0.286246	0.355638	0.42973	0.508842	0.593314	0.683507	0.779811	0.88264	0.992434	1.109666	1.23484	1.368494	1.511203	1.663579	1.826278	1.999998

計算結果は以下になる。

yが0の付近は格子が細かくて、2の付近は粗くなっていることがわかる。自分で確認すること。

ここで、座標変換係数について、もう一回考えてみよう。

2次元の場合、求めないといけない変数が6つ存在する。しかし、上のような格子だと、y方向に非等間隔であり座標変換空間のη方向しか関与しない。このように1つの方向だけを変える場合に必要な座標変換係数は2つだけでいい。上の格子条件だと、

になる項はΔxを使う一般の差分でいい。

yだけη方向に非等間隔になっているため、yに関する偏微分だけ座標変換係数を使えば計算ができる。メリットとしては直角座標系の計算ことから座標変換係数を導入するとき、修正が必要な部分が少なくなることから、バクが劇的に減る。また、計算上、参照する変数の数が減ることから、計算速度も速くなる。この時の座標変換係数は以下のように求められる。

取り合えず、基本形から考えてみると、

になる。理由としては、yはξと直交しているからである。結局、

となり、を求めたら、の差分ができる。の求め方は既存の求め方で求めればいい。

自ら求めてみたら、

のようになる。2次偏微分方程式は、

となり、を求めると、の差分もできる。も自ら求めてみると、

である。必ず自ら求めてみること。

6つの変数が全部必要な格子は以下のように、

r0 = (D0*0.5 + 0.004*(pow(1.0221, (double)(i)) - 1.0) / 0.0221);

x0[j][i] = r0 * cos((double)(j)*1.0 / rde);

y0[j][i] = r0 * sin((double)(j)*1.0 / rde);

円筒座標を作ったときである。

xはiとj、つまり、ξとηに関わっている。yもξとηに関わっている。この場合だと、6つの座標変換係数が必要である。

また、

x0[j][i] = 0.05*(pow(1.0677448,(double)(i))-1.0)/(0.0677448);

y0[j][i] = 0.05*(pow(1.0677448,(double)(j))-1.0)/(0.0677448);

このように、xはiと、yはjと関係していて、非等間隔であれは、4つの変数が必要である。