Commit babe9686 authored by ph's avatar ph

last bug found ?

parent a535d83a
......@@ -31,7 +31,7 @@
#define _LAMBDA _lambda_
#ifndef M_PI
#define M_PI (3.14159265358979323846264338328_F)
#define M_PI (3.14159265358979323846264338328)
#endif
#define _VOL (_DX * _DY)
......@@ -63,7 +63,7 @@ __constant real ds[4] = {_DY, _DY, _DX, _DX};
void flux_phy(const real *W, const real *vn, real *flux) {
real gam = 1.6666666666_F;
real gam = 1.6666666666;
real un = (W[eU] * vn[eX] + W[eV] * vn[eY] + W[eW] * vn[eZ]) / W[eRho];
real bn = W[eBx] * vn[eX] + W[eBy] * vn[eY] + W[eBz] * vn[eZ];
......@@ -142,7 +142,7 @@ void f2w(const real *f, real *w) {
// mhd
void conservatives(real *y, real *w) {
real gam = 1.6666666666_F;
real gam = 1.6666666666;
w[eRho] = y[eRho]; // rho
w[eU] = y[eRho] * y[eU]; // rho u
......@@ -159,7 +159,7 @@ void conservatives(real *y, real *w) {
// orszag tang init data
void exact_sol(real *x, real t, real *w) {
real gam = 1.6666666666_F;
real gam = 1.6666666666;
real yL[_M];
yL[eRho] = gam * gam; // rho
......@@ -175,27 +175,32 @@ void exact_sol(real *x, real t, real *w) {
conservatives(yL, w);
}
real gauss(real r){
return r < 6 ? exp(-r*r/2) : 0;
}
void exact_smooth_vortex(real *x, real t, real *w) {
const real gam = 1.6666666666_F;
const real gam = 1.6666666666;
const real rref2 = 1;
const real sigmam2 = 0.5_F;
const real sigmam2 = 0.5;
const real kappa = 1;
const real uref = 0.2_F;
const real bref = 0.2_F;
const real uref = 0.2;
const real bref = 0.2;
const real mu = 1;
const real udrift[2] = {1.0_F, 1.0_F};
const real xstart[2] = {5._F, 5._F};
const real udrift[2] = {1.0, 1.0};
const real xstart[2] = {5., 5.};
//
real xrel[3] = {x[0] - uref * udrift[0] * t - xstart[0],
x[1] - uref * udrift[1] * t - xstart[1], 0.0_F};
x[1] - uref * udrift[1] * t - xstart[1], 0.0};
//
real r = sqrt(xrel[0] * xrel[0] + xrel[1] * xrel[1]);
real expf = exp(sigmam2 * (1.0_F - r * r));
real expf = exp(sigmam2) * gauss(r);
// rho
w[eRho] = 1.0_F;
w[eRho] = 1.0;
// qx = rho * ux
w[eU] = uref * (udrift[0] - xrel[1] * kappa * expf);
// qy = rho * uy
......@@ -205,19 +210,19 @@ void exact_smooth_vortex(real *x, real t, real *w) {
// By
w[eBy] = bref * xrel[0] * mu * expf;
// psi
w[ePsi] = 0.0_F;
w[ePsi] = 0.0;
//
// real T = 1.0;
// real e = T/(gam-1.0);
//
real p = 1.0_F - 0.5_F * bref * bref * mu * mu * r * r * expf * expf;
real qdotq = w[1] * w[1] + w[2] * w[2];
real BdotB = w[4] * w[4] + w[5] * w[5];
real p = 1.0 - 0.5 * bref * bref * mu * mu * r * r * expf * expf;
real qdotq = w[eU] * w[eU] + w[eV] * w[eV];
real BdotB = w[eBx] * w[eBx] + w[eBy] * w[eBy];
//
w[eE] = p / (gam - 1.0_F) + 0.5_F * (qdotq / w[0] + BdotB);
w[eE] = p / (gam - 1.0) + 0.5 * (qdotq / w[eRho] + BdotB);
w[eW] = 0._F;
w[eBz] = 0._F;
w[eW] = 0.;
w[eBz] = 0.;
}
......@@ -307,18 +312,18 @@ __kernel void time_step(__global const real *fn, __global real *fnp1) {
// first order relaxation
w2f(wnow + 0, fnext + 0);
real om = 2._F;
real om = 1.9;
// second order relaxation
for (int iv = 0; iv < _M; iv++) {
for (int d = 0; d < 4; d++) {
int ik = d * _M + iv;
// if (iv == ePsi) {
// fnext[ik] = 1.0_F * fnext[ik];
// fnext[ik] = 1.0 * fnext[ik];
fnext[ik] = om * fnext[ik] - (om - 1) * fnow[ik];
//} else {
// fnext[ik] = 1.9_F * fnext[ik] - 0.9_F * fnow[ik];
// fnext[ik] = 2.0_F * fnext[ik] - 1.0_F * fnow[ik];
// fnext[ik] = 1.0_F * fnext[ik] - 0.0_F * fnow[ik];
// fnext[ik] = 1.9 * fnext[ik] - 0.9 * fnow[ik];
// fnext[ik] = 2.0 * fnext[ik] - 1.0 * fnow[ik];
// fnext[ik] = 1.0 * fnext[ik] - 0.0 * fnow[ik];
//}
}
}
......@@ -386,8 +391,8 @@ __kernel void kinetic(const real t, __global real *fn, __global real *kin) {
real floc[_N];
// save
int imem = i + j * _NX;
kin[imem] = 0.0_F;
int imemk = i + j * _NX;
kin[imemk] = 0;
for (int ik = 0; ik < _N; ik++) {
int imem = i + j * _NX + ik * ngrid;
......@@ -402,5 +407,5 @@ __kernel void kinetic(const real t, __global real *fn, __global real *kin) {
real kinloc = (wn[eU]-wnow[eU]) * (wn[eU]-wnow[eU]);
//kinloc += wn[eV] * wn[eV];
//kinloc *= 0.5 / wn[eRho];
kin[imem] += kinloc * _DX * _DY;
kin[imemk] += kinloc * _DX * _DY;
}
......@@ -59,7 +59,7 @@ def solve_ocl(m=_m, n=_n, nx=_nx, ny=_ny, Lx=_Lx, Ly=_Ly, Tmax=_Tmax,
dy = Ly / ny
# lattice speed
vmax = 12
vmax = 20
# time stepping
cfl = 1
......@@ -104,8 +104,15 @@ def solve_ocl(m=_m, n=_n, nx=_nx, ny=_ny, Lx=_Lx, Ly=_Ly, Tmax=_Tmax,
ctx = cl.create_some_context()
mf = cl.mem_flags
#opt = "-cl-denorms-are-zero -cl-fast-relaxed-math"
opt = ""
if precision == "single":
opt = opt + " -cl-single-precision-constant"
print(opt)
# compile OpenCL C program
prg = cl.Program(ctx, source).build(options="")
prg = cl.Program(ctx, source).build(options= opt)
# create a queue (for submitting opencl operations)
queue = cl.CommandQueue(
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment