smyalygames/FiniteVolumeGPU
1
0
Fork 0
mirror of https://github.com/smyalygames/FiniteVolumeGPU.git synced 2025-07-05 07:51:00 +02:00
Code Issues Actions Packages Projects Releases Wiki Activity

Implemented RT instability

Browse Source
This commit is contained in:
André R. Brodtkorb 2018-11-07 11:06:45 +01:00
parent 0f68c7867b
commit ae668a40d3
3 changed files with 170 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

172
EulerTesting.ipynb
View File

@ -37,11 +37,6 @@
" from StringIO import StringIO\n",
"except ImportError:\n",
" from io import StringIO\n",
" \n",
"#Set large figure sizes as default\n",
"rc('figure', figsize=(16.0, 12.0))\n",
"rc('animation', html='html5')\n",
"rc('animation', bitrate=1800)\n",
"\n",
"from GPUSimulators import Common, IPythonMagic"
]
@ -64,6 +59,18 @@
"%cuda_context_handler --no_autotuning my_context "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#Set large figure sizes as default\n",
"rc('figure', figsize=(16.0, 12.0))\n",
"rc('animation', html='html5')\n",
"rc('animation', bitrate=1800)"
]
},
{
"cell_type": "code",
"execution_count": null,
@ -127,7 +134,7 @@
"metadata": {},
"outputs": [],
"source": [
"def genKelvinHelmholz(nx, ny, gamma, roughness=0.5):\n",
"def genKelvinHelmholtz(nx, ny, gamma, roughness=0.5):\n",
" \"\"\"\n",
" Roughness parameter in (0, 1.0] determines how squiggly the interface betweeen the zones is\n",
" \"\"\"\n",
@ -211,12 +218,76 @@
" return rho, rho*u, rho*v, E, dx, dy, dt\n",
"\n",
"nx = 400\n",
"rho0, rho_u0, rho_v0, E0, dx, dy, dt = genKelvinHelmholz(nx, nx//4, 1.4, 0.25)\n",
"rho0, rho_u0, rho_v0, E0, dx, dy, dt = genKelvinHelmholtz(nx, nx//4, 1.4, 0.25)\n",
"plt.figure()\n",
"plt.imshow(rho0)\n",
"plt.colorbar(orientation='vertical', aspect=20, pad=0.02, shrink=0.3)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def genRayleighTaylor(nx, ny, gamma):\n",
" width = 0.5\n",
" height = 1.5\n",
" dx = width / float(nx)\n",
" dy = height / float(ny)\n",
" g = 0.1\n",
"\n",
" rho = np.zeros((ny, nx), dtype=np.float32)\n",
" u = np.zeros((ny, nx), dtype=np.float32)\n",
" v = np.zeros((ny, nx), dtype=np.float32)\n",
" p = np.zeros((ny, nx), dtype=np.float32)\n",
" \n",
" x = dx*np.arange(0, nx, dtype=np.float32)-width*0.5\n",
" y = dy*np.arange(0, ny, dtype=np.float32)-height*0.5\n",
" xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')\n",
" \n",
" #y_threshold = -0.01*np.cos(2*np.pi*np.abs(x)/0.5)\n",
" #rho = np.where(yv <= y_threshold, 1.0, rho)\n",
" #rho = np.where(yv > y_threshold, 2.0, rho)\n",
" \n",
" rho = np.where(yv <= 0.0, 1.0, rho)\n",
" rho = np.where(yv > 0.0, 2.0, rho)\n",
" \n",
" v = 0.01*(1.0 + np.cos(2*np.pi*xv/0.5))/4\n",
" \n",
" p = 2.5 - rho*g*yv\n",
" E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)\n",
" \n",
" #Estimate dt\n",
" scale = 0.9\n",
" max_rho_estimate = 3.0\n",
" max_u_estimate = 1.0\n",
" dx = width/nx\n",
" dy = height/ny\n",
" dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate))\n",
"\n",
" return rho, rho*u, rho*v, E, dx, dy, dt\n",
"\n",
"nx = 400\n",
"rho0, rho_u0, rho_v0, E0, dx, dy, dt = genRayleighTaylor(nx, nx*3, 1.4)\n",
"plt.figure()\n",
"plt.subplot(1,4,1)\n",
"plt.imshow(rho0, origin='bottom')\n",
"plt.colorbar(orientation='horizontal', aspect=20, pad=0.02, shrink=0.9)\n",
"\n",
"plt.subplot(1,4,2)\n",
"plt.imshow(rho_u0, origin='bottom')\n",
"plt.colorbar(orientation='horizontal', aspect=20, pad=0.02, shrink=0.9)\n",
"\n",
"plt.subplot(1,4,3)\n",
"plt.imshow(rho_v0, origin='bottom')\n",
"plt.colorbar(orientation='horizontal', aspect=20, pad=0.02, shrink=0.9)\n",
"\n",
"plt.subplot(1,4,4)\n",
"plt.imshow(E0, origin='bottom')\n",
"plt.colorbar(orientation='horizontal', aspect=20, pad=0.02, shrink=0.9)"
]
},
{
"cell_type": "code",
"execution_count": null,
@ -330,8 +401,8 @@
" Schlieren = 0\n",
" Density = 1\n",
"\n",
"def makeAnimation(infile, vis_type, vmin, vmax, save_anim=False, cmap=plt.cm.coolwarm):\n",
" fig = plt.figure(figsize=(16, 8))#, tight_layout=True)\n",
"def makeAnimation(infile, vis_type, vmin, vmax, save_anim=False, cmap=plt.cm.coolwarm, fig_args=None):\n",
" fig = plt.figure(**fig_args)\n",
" \n",
" with Common.DataDumper(infile, 'r') as indata:\n",
" time = indata.ncfile.variables['time'][:]\n",
@ -397,8 +468,10 @@
" print(\"Reusing previously created movie \" + movie_outpath)\n",
" else:\n",
" print(\"Creating movie \" + movie_outpath)\n",
" from matplotlib.animation import FFMpegFileWriter\n",
" writer = FFMpegFileWriter(fps=25)\n",
" #from matplotlib.animation import FFMpegFileWriter\n",
" #writer = FFMpegFileWriter(fps=25)\n",
" from matplotlib.animation import FFMpegWriter\n",
" writer = FFMpegWriter(fps=25)\n",
" anim.save(movie_outpath, writer=writer)\n",
" display(HTML(\"\"\"\n",
" <div align=\"middle\">\n",
@ -413,6 +486,44 @@
" display(anim)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"* Liska and Wendroff, COMPARISON OF SEVERAL DIFFERENCE SCHEMES ON 1D AND 2D TEST PROBLEMS FOR THE EULER EQUATIONS, http://www-troja.fjfi.cvut.cz/~liska/CompareEuler/compare8.pdf\n",
"* https://www.astro.princeton.edu/~jstone/Athena/tests/rt/rt.html\n",
"* "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nx = 301\n",
"ny = nx*3\n",
"gamma = 1.4\n",
"t_end = 30\n",
"n_save = 300\n",
"outfile = \"data/euler_rayleigh-taylor.nc\"\n",
"outdata = None\n",
"\n",
"rho0, rho_u0, rho_v0, E0, dx, dy, dt = genRayleighTaylor(nx, ny, gamma)\n",
"outfile = makeSimulation(outfile, t_end, rho0, rho_u0, rho_v0, E0, nx, ny, dx, dy, dt, gamma)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#outfile = 'data/euler_rayleigh-taylor_0007.nc'\n",
"makeAnimation(outfile, vis_type=VisType.Density, vmin=1, vmax=2, cmap=plt.cm.RdBu, save_anim=True, fig_args={'figsize':(3.4, 8), 'tight_layout':True})"
]
},
{
"cell_type": "code",
"execution_count": null,
@ -422,8 +533,8 @@
"nx = 1600\n",
"ny = nx//4\n",
"gamma = 1.4\n",
"t_end = 3.0\n",
"n_save = 200\n",
"t_end = 0.5#3.0\n",
"n_save = 20#500\n",
"outfile = \"data/euler_shock-bubble.nc\"\n",
"outdata = None\n",
"\n",
@ -437,8 +548,8 @@
"metadata": {},
"outputs": [],
"source": [
"#makeAnimation('data/euler_simdata_0000.nc')\n",
"makeAnimation(outfile, vis_type=VisType.Schlieren, vmin=-0.2, vmax=0.2, cmap=plt.cm.RdBu)"
"#outfile = 'data/euler_shock-bubble_0011.nc'\n",
"makeAnimation(outfile, vis_type=VisType.Schlieren, vmin=-0.2, vmax=0.2, cmap=plt.cm.RdBu, save_anim=True, fig_args={'figsize':(16, 5), 'tight_layout':True})"
]
},
{
@ -451,12 +562,12 @@
"ny = nx//2\n",
"gamma = 1.4\n",
"roughness = 0.125\n",
"t_end = 2.0#10.0\n",
"n_save = 100#1000\n",
"outfile = \"data/euler_kelvin_helmholz.nc\"\n",
"t_end = 10.0\n",
"n_save = 1000\n",
"outfile = \"data/euler_kelvin_helmholtz.nc\"\n",
"outdata = None\n",
"\n",
"rho0, rho_u0, rho_v0, E0, dx, dy, dt = genKelvinHelmholz(nx, ny, gamma, roughness)\n",
"rho0, rho_u0, rho_v0, E0, dx, dy, dt = genKelvinHelmholtz(nx, ny, gamma, roughness)\n",
"outfile = makeSimulation(outfile, t_end, rho0, rho_u0, rho_v0, E0, nx, ny, dx, dy, dt, gamma)"
]
},
@ -466,27 +577,10 @@
"metadata": {},
"outputs": [],
"source": [
"makeAnimation(outfile, vis_type=VisType.Density, vmin=1.0, vmax=2.0)\n",
"#makeAnimation('data/euler_kelvin_helmholz_0012.nc', cmap=plt.cm.coolwarm, save_anim=True, vmin=1.1, vmax=1.9)\n",
"#makeAnimation('data/euler_kelvin_helmholz_0011.nc', save_anim=True)"
"#outfile='data/euler_kelvin_helmholtz_0012.nc'\n",
"makeAnimation(outfile, vis_type=VisType.Density, vmin=1.0, vmax=2.0, save_anim=True, fig_args={'figsize':(16, 9), 'tight_layout':True})"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"makeAnimation(outfile, vis_type=VisType.Schlieren, vmin=-0.2, vmax=0.2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
@ -511,7 +605,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
"version": "3.6.7"
}
},
"nbformat": 4,

30
GPUSimulators/cuda/EE2D_KP07_dimsplit.cu
View File

@ -59,6 +59,7 @@ void computeFluxF(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
// Compute flux based on prediction
const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
//const float4 flux = HLL_flux(Q_l_bar, Q_r_bar, gamma_);
//Write to shared memory
F[0][j][i] = flux.x;
@ -103,6 +104,7 @@ void computeFluxG(float Q[4][BLOCK_HEIGHT+4][BLOCK_WIDTH+4],
// Compute flux based on prediction
const float4 flux = CentralUpwindFlux(Q_l_bar, Q_r_bar, gamma_);
//const float4 flux = HLL_flux(Q_l_bar, Q_r_bar, gamma_);
//Write to shared memory
//Note that we here swap hu and hv back to the original
@ -162,10 +164,12 @@ __global__ void KP07DimsplitKernel(
//Fix boundary conditions
noFlowBoundary<w, h, gc, 1, 1>(Q[0], nx_, ny_);
noFlowBoundary<w, h, gc, 1, 1>(Q[1], nx_, ny_);
noFlowBoundary<w, h, gc, -1, 1>(Q[1], nx_, ny_);
noFlowBoundary<w, h, gc, 1, -1>(Q[2], nx_, ny_);
noFlowBoundary<w, h, gc, 1, 1>(Q[3], nx_, ny_);
__syncthreads();
const float g = 0.1f;
//Step 0 => evolve x first, then y
@ -182,7 +186,7 @@ __global__ void KP07DimsplitKernel(
//Set boundary conditions
noFlowBoundary<w, h, gc, 1, 1>(Q[0], nx_, ny_);
noFlowBoundary<w, h, gc, 1, 1>(Q[1], nx_, ny_);
noFlowBoundary<w, h, gc, -1, 1>(Q[1], nx_, ny_);
noFlowBoundary<w, h, gc, 1, -1>(Q[2], nx_, ny_);
noFlowBoundary<w, h, gc, 1, 1>(Q[3], nx_, ny_);
__syncthreads();
@ -196,6 +200,16 @@ __global__ void KP07DimsplitKernel(
evolveG<w, h, gc, vars>(Q, F, dy_, dt_);
__syncthreads();
//Gravity source term
{
const int i = threadIdx.x + gc;
const int j = threadIdx.y + gc;
const float rho_v = Q[2][j][i];
Q[2][j][i] -= g*Q[0][j][i]*dt_;
Q[3][j][i] -= g*rho_v*dt_;
}
__syncthreads();
}
//Step 1 => evolve y first, then x
@ -212,7 +226,7 @@ __global__ void KP07DimsplitKernel(
//Set boundary conditions
noFlowBoundary<w, h, gc, 1, 1>(Q[0], nx_, ny_);
noFlowBoundary<w, h, gc, 1, 1>(Q[1], nx_, ny_);
noFlowBoundary<w, h, gc, -1, 1>(Q[1], nx_, ny_);
noFlowBoundary<w, h, gc, 1, -1>(Q[2], nx_, ny_);
noFlowBoundary<w, h, gc, 1, 1>(Q[3], nx_, ny_);
__syncthreads();
@ -227,6 +241,16 @@ __global__ void KP07DimsplitKernel(
evolveF<w, h, gc, vars>(Q, F, dx_, dt_);
__syncthreads();
//Gravity source term
{
const int i = threadIdx.x + gc;
const int j = threadIdx.y + gc;
const float rho_v = Q[2][j][i];
Q[2][j][i] -= g*Q[0][j][i]*dt_;
Q[3][j][i] -= g*rho_v*dt_;
}
__syncthreads();
//This is the RK2-part
const int tx = threadIdx.x + gc;
const int ty = threadIdx.y + gc;

14
GPUSimulators/cuda/common.h
View File

@ -102,19 +102,25 @@ inline __device__ void readBlock(float* ptr_, int pitch_,
//Read into shared memory
//Loop over all variables
for (int j=threadIdx.y; j<block_height+2*ghost_cells; j+=block_height) {
//const int l = min(by + j, ny_+2*ghost_cells-1);
const int l = min(by + j, ny_+2*ghost_cells-1);
/*
const int y = by + j;
const int y_offset = ( (int) (y < gc_pad) - (int) (y >= ny_+2*ghost_cells-gc_pad) ) * (ny_+2*ghost_cells - 2*gc_pad);
const int l = y + y_offset;
const int l = min(y + y_offset, ny_+2*ghost_cells-1);
*/
float* row = (float*) ((char*) ptr_ + pitch_*l);
for (int i=threadIdx.x; i<block_width+2*ghost_cells; i+=block_width) {
//const int k = min(bx + i, nx_+2*ghost_cells-1);
const int k = min(bx + i, nx_+2*ghost_cells-1);
/*
const int x = bx + i;
const int gc_pad = 4;
const int x_offset = ( (int) (x < gc_pad) - (int) (x >= nx_+2*ghost_cells-gc_pad) ) * (nx_+2*ghost_cells - 2*gc_pad);
const int k = x + x_offset;
const int k = min(x + x_offset, nx_+2*ghost_cells-1);
*/
shmem[j][i] = row[k];
}