Commit 38105574 authored by ph's avatar ph

up

parent 30220022
......@@ -66,7 +66,7 @@ def solve_ocl(m=_m, nx=_nx, ny=_ny, Tmax=_Tmax, Lx=_Lx, Ly=_Ly, animate=True, iv
mf = cl.mem_flags
# compile OpenCL C program
prg = cl.Program(ctx, source).build(options="-cl-fast-relaxed-math")
prg = cl.Program(ctx, source).build(options="-cl-fast-relaxed-math -cl-single-precision-constant")
dev = ctx.devices[0];
......
......@@ -5,6 +5,9 @@
#define _DT _dt_
#define _M _m_
#define double float
#define _LAMBDA _lambda_
#define _G (9.81_F)
......@@ -13,9 +16,7 @@
__constant int dir[4][2] = {{1, 0}, {-1, 0}, {0, 1}, {0, -1}};
__constant float ds[4] = {_DY, _DY, _DX, _DX};
#define double float
__constant double ds[4] = {_DY, _DY, _DX, _DX};
void riem_stvenant(double *wL, double *wR, double xi, double *w);
......@@ -25,22 +26,22 @@ double dZ(double h1, double h2);
void flux_riem_2d(double *wL, double *wR, double *vn, double *flux);
void fluxphy(float *w, float *n, float *flux) {
void fluxphy(double *w, double *n, double *flux) {
float h = w[0];
float u = w[1] / h;
float v = w[2] / h;
double h = w[0];
double u = w[1] / h;
double v = w[2] / h;
float un = u * n[0] + v * n[1];
double un = u * n[0] + v * n[1];
flux[0] = h * un;
flux[1] = h * u * un + _G * h * h / 2 * n[0];
flux[2] = h * v * un + _G * h * h / 2 * n[1];
}
void fluxnum(float *wL, float *wR, float *vnorm, float *flux) {
float fL[_M];
float fR[_M];
void fluxnum(double *wL, double *wR, double *vnorm, double *flux) {
double fL[_M];
double fR[_M];
fluxphy(wL, vnorm, fL);
fluxphy(wR, vnorm, fR);
......@@ -189,7 +190,7 @@ double Z(double hs, double h) {
double dZ(double hs, double h) {
double t0 =
-sqrt(g) / pow(sqrt(hs) + sqrt(h), 2.0f) * Heaviside(h - hs) / sqrt(hs) -
-sqrt(g) / pow(sqrt(hs) + sqrt(h), 2.0) * Heaviside(h - hs) / sqrt(hs) -
2.0 * sqrt(g) / (sqrt(hs) + sqrt(h)) * Dirac(-h + hs) +
sqrt(2.0) * sqrt(g) / sqrt(h + hs) / sqrt(h * hs) / 4.0 -
sqrt(2.0) * sqrt(g) * sqrt(h + hs) / sqrt(h * h * h * hs * hs * hs) * h /
......@@ -203,10 +204,10 @@ double dZ(double hs, double h) {
return t0;
}
void exact_sol(float *xy, float t, float *w) {
float h = 1;
float x = xy[0]; // - 0.5f;
float y = xy[1]; // - 0.5f;
void exact_sol(double *xy, double t, double *w) {
double h = 1;
double x = xy[0]; // - 0.5f;
double y = xy[1]; // - 0.5f;
if (x < 0.75f && x > 0.25f && y < 0.75f && y > 0.25f)
h = 2;
......@@ -216,7 +217,7 @@ void exact_sol(float *xy, float t, float *w) {
w[2] = 0;
}
__kernel void init_sol(__global float *wn) {
__kernel void init_sol(__global double *wn) {
int id = get_global_id(0);
......@@ -225,10 +226,10 @@ __kernel void init_sol(__global float *wn) {
int ngrid = _NX * _NY;
float wnow[_M];
double wnow[_M];
float t = 0;
float xy[2] = {i * _DX + _DX / 2, j * _DY + _DY / 2};
double t = 0;
double xy[2] = {i * _DX + _DX / 2, j * _DY + _DY / 2};
exact_sol(xy, t, wnow);
// printf("x=%f, y=%f \n",xy[0],xy[1]);
......@@ -242,7 +243,34 @@ __kernel void init_sol(__global float *wn) {
}
}
__kernel void time_step(__global float *wn, __global float *wnp1) {
double minmod(double a, double b, double c)
{
double res = 0.0;
if ((a > 0) && (b > 0) && (c > 0))
res = a < b ? (a < c ? a : c) : (b < c ? b : c);
else if ((a < 0) && (b < 0) && (c < 0))
res = a > b ? (a > c ? a : c) : (b > c ? b : c);
return res;
}
__kernel void muscl(__global const double *wn, __global double *dxwn,
__global double *dywn, __global double *dtwn) {
int id = get_global_id(0);
int i = id % _NX;
int j = id / _NX;
int ngrid = _NX * _NY;
// slope exists only for internal cells
if (i > 0 && i < _NX - 1 && j > 0 && j < _NY - 1) {
// etc...
}
}
__kernel void time_step(__global const double *wn, __global double *wnp1) {
int id = get_global_id(0);
......@@ -251,19 +279,19 @@ __kernel void time_step(__global float *wn, __global float *wnp1) {
int ngrid = _NX * _NY;
float wnow[_M];
float wnext[_M];
double wnow[_M];
double wnext[_M];
// load middle value
for (int iv = 0; iv < _M; iv++) {
int imem = i + j * _NX + iv * ngrid;
wnow[iv] = wn[imem]; // moralement: wn[iv][j][i]
wnow[iv] = wn[imem]; // moralement: wn[iv][j][i]
wnext[iv] = wnow[iv];
}
// only compute internal cells
if (i > 0 && i < _NX - 1 && j > 0 && j < _NY - 1) {
float flux[_M];
double flux[_M];
// loop on all directions:
// idir = 0 (east)
......@@ -271,25 +299,25 @@ __kernel void time_step(__global float *wn, __global float *wnp1) {
// idir = 2 (north)
// idir = 3 (south)
for (int idir = 0; idir < 4; idir++) {
float vn[2];
double vn[2];
vn[0] = dir[idir][0];
vn[1] = dir[idir][1];
int iR = i + dir[idir][0];
int jR = j + dir[idir][1];
// load neighbour values
float wR[_M];
double wR[_M];
for (int iv = 0; iv < _M; iv++) {
int imem = iv * ngrid + iR + jR * _NX;
wR[iv] = wn[imem];
}
if (iR == 0 || iR == _NX - 1 || jR == 0 || jR == _NY - 1) {
wR[0] = wnow[0];
float qn = wnow[1] * vn[0] + wnow[2] * vn[1];
double qn = wnow[1] * vn[0] + wnow[2] * vn[1];
wR[1] = wnow[1] - 2 * qn * vn[0];
wR[2] = wnow[2] - 2 * qn * vn[1];
}
flux_riem_2d(wnow, wR, vn, flux);
//fluxnum(wnow, wR, vn, flux);
// flux_riem_2d(wnow, wR, vn, flux);
fluxnum(wnow, wR, vn, flux);
// time evolution
for (int iv = 0; iv < _M; iv++) {
......
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