Commit a9cd0e12 authored by Matthieu Boileau's avatar Matthieu Boileau

fix vortex_lbm

parent 032043ed
......@@ -65,12 +65,12 @@ typedef bessel_j_scalar_type T;
// {{{ bessel_j0
__constant const bessel_j_scalar_type bessel_j0_P1[] = {
__constant bessel_j_scalar_type bessel_j0_P1[] = {
-4.1298668500990866786e+11, 2.7282507878605942706e+10,
-6.2140700423540120665e+08, 6.6302997904833794242e+06,
-3.6629814655107086448e+04, 1.0344222815443188943e+02,
-1.2117036164593528341e-01};
__constant const bessel_j_scalar_type bessel_j0_Q1[] = {
__constant bessel_j_scalar_type bessel_j0_Q1[] = {
2.3883787996332290397e+12,
2.6328198300859648632e+10,
1.3985097372263433271e+08,
......@@ -78,29 +78,29 @@ __constant const bessel_j_scalar_type bessel_j0_Q1[] = {
9.3614022392337710626e+02,
1.0,
0.0};
__constant const bessel_j_scalar_type bessel_j0_P2[] = {
__constant bessel_j_scalar_type bessel_j0_P2[] = {
-1.8319397969392084011e+03, -1.2254078161378989535e+04,
-7.2879702464464618998e+03, 1.0341910641583726701e+04,
1.1725046279757103576e+04, 4.4176707025325087628e+03,
7.4321196680624245801e+02, 4.8591703355916499363e+01};
__constant const bessel_j_scalar_type bessel_j0_Q2[] = {
__constant bessel_j_scalar_type bessel_j0_Q2[] = {
-3.5783478026152301072e+05, 2.4599102262586308984e+05,
-8.4055062591169562211e+04, 1.8680990008359188352e+04,
-2.9458766545509337327e+03, 3.3307310774649071172e+02,
-2.5258076240801555057e+01, 1.0};
__constant const bessel_j_scalar_type bessel_j0_PC[] = {
__constant bessel_j_scalar_type bessel_j0_PC[] = {
2.2779090197304684302e+04, 4.1345386639580765797e+04,
2.1170523380864944322e+04, 3.4806486443249270347e+03,
1.5376201909008354296e+02, 8.8961548424210455236e-01};
__constant const bessel_j_scalar_type bessel_j0_QC[] = {
__constant bessel_j_scalar_type bessel_j0_QC[] = {
2.2779090197304684318e+04, 4.1370412495510416640e+04,
2.1215350561880115730e+04, 3.5028735138235608207e+03,
1.5711159858080893649e+02, 1.0};
__constant const bessel_j_scalar_type bessel_j0_PS[] = {
__constant bessel_j_scalar_type bessel_j0_PS[] = {
-8.9226600200800094098e+01, -1.8591953644342993800e+02,
-1.1183429920482737611e+02, -2.2300261666214198472e+01,
-1.2441026745835638459e+00, -8.8033303048680751817e-03};
__constant const bessel_j_scalar_type bessel_j0_QS[] = {
__constant bessel_j_scalar_type bessel_j0_QS[] = {
5.7105024128512061905e+03, 1.1951131543434613647e+04,
7.2642780169211018836e+03, 1.4887231232283756582e+03,
9.0593769594993125859e+01, 1.0};
......@@ -149,12 +149,12 @@ bessel_j_scalar_type bessel_j0(bessel_j_scalar_type x) {
// {{{ bessel_j1
__constant const bessel_j_scalar_type bessel_j1_P1[] = {
__constant bessel_j_scalar_type bessel_j1_P1[] = {
-1.4258509801366645672e+11, 6.6781041261492395835e+09,
-1.1548696764841276794e+08, 9.8062904098958257677e+05,
-4.4615792982775076130e+03, 1.0650724020080236441e+01,
-1.0767857011487300348e-02};
__constant const bessel_j_scalar_type bessel_j1_Q1[] = {
__constant bessel_j_scalar_type bessel_j1_Q1[] = {
4.1868604460820175290e+12,
4.2091902282580133541e+10,
2.0228375140097033958e+08,
......@@ -162,17 +162,17 @@ __constant const bessel_j_scalar_type bessel_j1_Q1[] = {
1.0742272239517380498e+03,
1.0,
0.0};
__constant const bessel_j_scalar_type bessel_j1_P2[] = {
__constant bessel_j_scalar_type bessel_j1_P2[] = {
-1.7527881995806511112e+16, 1.6608531731299018674e+15,
-3.6658018905416665164e+13, 3.5580665670910619166e+11,
-1.8113931269860667829e+09, 5.0793266148011179143e+06,
-7.5023342220781607561e+03, 4.6179191852758252278e+00};
__constant const bessel_j_scalar_type bessel_j1_Q2[] = {
__constant bessel_j_scalar_type bessel_j1_Q2[] = {
1.7253905888447681194e+18, 1.7128800897135812012e+16,
8.4899346165481429307e+13, 2.7622777286244082666e+11,
6.4872502899596389593e+08, 1.1267125065029138050e+06,
1.3886978985861357615e+03, 1.0};
__constant const bessel_j_scalar_type bessel_j1_PC[] = {
__constant bessel_j_scalar_type bessel_j1_PC[] = {
-4.4357578167941278571e+06,
-9.9422465050776411957e+06,
-6.6033732483649391093e+06,
......@@ -180,7 +180,7 @@ __constant const bessel_j_scalar_type bessel_j1_PC[] = {
-1.0982405543459346727e+05,
-1.6116166443246101165e+03,
0.0};
__constant const bessel_j_scalar_type bessel_j1_QC[] = {
__constant bessel_j_scalar_type bessel_j1_QC[] = {
-4.4357578167941278568e+06,
-9.9341243899345856590e+06,
-6.5853394797230870728e+06,
......@@ -188,7 +188,7 @@ __constant const bessel_j_scalar_type bessel_j1_QC[] = {
-1.0726385991103820119e+05,
-1.4550094401904961825e+03,
1.0};
__constant const bessel_j_scalar_type bessel_j1_PS[] = {
__constant bessel_j_scalar_type bessel_j1_PS[] = {
3.3220913409857223519e+04,
8.5145160675335701966e+04,
6.6178836581270835179e+04,
......@@ -196,7 +196,7 @@ __constant const bessel_j_scalar_type bessel_j1_PS[] = {
1.7063754290207680021e+03,
3.5265133846636032186e+01,
0.0};
__constant const bessel_j_scalar_type bessel_j1_QS[] = {
__constant bessel_j_scalar_type bessel_j1_QS[] = {
7.0871281941028743574e+05,
1.8194580422439972989e+06,
1.4194606696037208929e+06,
......
......@@ -15,9 +15,7 @@
//#pragma OPENCL EXTENSION cl_khr_fp64 : enable
#include "vortex_kernels.h"
#define real _real_
#define _NX _nx_
#define _NY _ny_
#define _DX _dx_
......@@ -56,6 +54,8 @@
//#define dirichlet_updown
//#define dirichlet_leftright
void flux_phy(const real *W, const real *vn, real *flux);
__constant int dir[4][2] = {{-1, 0}, {1, 0}, {0, -1}, {0, 1}};
__constant real ds[4] = {_DY, _DY, _DX, _DX};
......@@ -176,11 +176,7 @@ void exact_sol(real *x, real t, real *w) {
conservatives(yL, w);
}
real gauss(real r){
return exp(-r*r/2);
}
real gauss(real r) { return exp(-r * r / 2); }
void exact_smooth_vortex(real *x, real t, real *w) {
......@@ -224,11 +220,8 @@ void exact_smooth_vortex(real *x, real t, real *w) {
w[eW] = 0.;
w[eBz] = 0.;
}
// initial condition on the macro data
__kernel void init_sol(__global real *fn) {
......@@ -244,10 +237,9 @@ __kernel void init_sol(__global real *fn) {
real t = 0;
real xy[2] = {i * _DX + _DX / 2, j * _DY + _DY / 2};
//exact_sol(xy, t, wnow);
// exact_sol(xy, t, wnow);
exact_smooth_vortex(xy, t, wnow);
real fnow[_N];
w2f(wnow, fnow);
......@@ -260,9 +252,6 @@ __kernel void init_sol(__global real *fn) {
}
}
// one time step of the LBM scheme
__kernel void time_step(__global const real *fn, __global real *fnp1) {
......@@ -283,25 +272,25 @@ __kernel void time_step(__global const real *fn, __global real *fnp1) {
for (int iv = 0; iv < _M; iv++) {
int ik = d * _M + iv;
int imem = iR + jR * _NX + ik * ngrid;
// #ifdef dirichlet_updown
// // dirichlet condition on up and down borders
// // (values of border cells are unchanged)
// if ((j == 0) || (j == _NY - 1)) {
// imem = i + j * _NX + ik * ngrid;
// }
// #elif defined dirichlet_leftright
// // dirichlet condition on left and right borders
// // (values of border cells are unchanged)
// if ((i == 0) || (i == _NX - 1)) {
// imem = i + j * _NX + ik * ngrid;
// }
// #elif defined dirichlet
// // dirichlet condition on all borders
// // (values of border cells are unchanged)
// if ((i == 0) || (i == _NX - 1) || (j == 0) || (j == _NY - 1)) {
// imem = i + j * _NX + ik * ngrid;
// }
// #endif
// #ifdef dirichlet_updown
// // dirichlet condition on up and down borders
// // (values of border cells are unchanged)
// if ((j == 0) || (j == _NY - 1)) {
// imem = i + j * _NX + ik * ngrid;
// }
// #elif defined dirichlet_leftright
// // dirichlet condition on left and right borders
// // (values of border cells are unchanged)
// if ((i == 0) || (i == _NX - 1)) {
// imem = i + j * _NX + ik * ngrid;
// }
// #elif defined dirichlet
// // dirichlet condition on all borders
// // (values of border cells are unchanged)
// if ((i == 0) || (i == _NX - 1) || (j == 0) || (j == _NY - 1)) {
// imem = i + j * _NX + ik * ngrid;
// }
// #endif
fnow[ik] = fn[imem];
}
}
......@@ -405,9 +394,9 @@ __kernel void kinetic(const real t, __global real *fn, __global real *kin) {
real xy[2] = {i * _DX + _DX / 2, j * _DY + _DY / 2};
real wnow[_M];
exact_smooth_vortex(xy, t, wnow);
//real kinloc = (wn[eU]-wnow[eU]) * (wn[eU]-wnow[eU]);
real kinloc = fabs(wn[eU]-wnow[eU]);
//kinloc += wn[eV] * wn[eV];
//kinloc *= 0.5 / wn[eRho];
// real kinloc = (wn[eU]-wnow[eU]) * (wn[eU]-wnow[eU]);
real kinloc = fabs(wn[eU] - wnow[eU]);
// kinloc += wn[eV] * wn[eV];
// kinloc *= 0.5 / wn[eRho];
kin[imemk] += kinloc * _DX * _DY;
}
void flux_phy(const real *W, const real *vn, real *flux);
......@@ -106,14 +106,14 @@ def solve_ocl(m=_m, n=_n, nx=_nx, ny=_ny, Lx=_Lx, Ly=_Ly, Tmax=_Tmax,
mf = cl.mem_flags
#opt = "-cl-denorms-are-zero -cl-fast-relaxed-math"
opt = ""
opt = []
if precision == "single":
opt = opt + " -cl-single-precision-constant"
opt.append("-cl-single-precision-constant")
print(opt)
# compile OpenCL C program
prg = cl.Program(ctx, source).build(options= opt)
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