mirror of
https://github.com/smyalygames/FiniteVolumeGPU.git
synced 2025-11-27 22:16:14 +01:00
463 lines
16 KiB
Plaintext
463 lines
16 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"#Lets have matplotlib \"inline\"\n",
|
|
"%matplotlib inline\n",
|
|
"\n",
|
|
"# Add line profiler\n",
|
|
"%load_ext line_profiler\n",
|
|
"\n",
|
|
"#Import packages we need\n",
|
|
"import numpy as np\n",
|
|
"from matplotlib import animation, rc\n",
|
|
"from matplotlib import pyplot as plt\n",
|
|
"from matplotlib.colorbar import ColorbarBase\n",
|
|
"\n",
|
|
"from IPython.display import display, HTML\n",
|
|
"\n",
|
|
"from enum import Enum\n",
|
|
"import time\n",
|
|
"import os\n",
|
|
"import json\n",
|
|
"\n",
|
|
"try:\n",
|
|
" from StringIO import StringIO\n",
|
|
"except ImportError:\n",
|
|
" from io import StringIO\n",
|
|
"\n",
|
|
"from GPUSimulators.common import Timer, DataDumper, ProgressPrinter\n",
|
|
"from GPUSimulators.model import EE2DKP07Dimsplit\n",
|
|
"from GPUSimulators.helpers import initial_conditions, visualization"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": "%setup_logging --out test_schemes.log logger"
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"metadata": {},
|
|
"source": [
|
|
"%cuda_context_handler --no_autotuning my_context "
|
|
],
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"metadata": {},
|
|
"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)"
|
|
],
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"def plot_vars(rho, rho_u, rho_v, E):\n",
|
|
" plt.figure()\n",
|
|
" plt.subplot(1, 4, 1)\n",
|
|
" plt.imshow(rho, origin='lower')\n",
|
|
" plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n",
|
|
"\n",
|
|
" plt.subplot(1, 4, 2)\n",
|
|
" plt.imshow(rho_u, origin='lower')\n",
|
|
" plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n",
|
|
"\n",
|
|
" plt.subplot(1, 4, 3)\n",
|
|
" plt.imshow(rho_v, origin='lower')\n",
|
|
" plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)\n",
|
|
"\n",
|
|
" plt.subplot(1, 4, 4)\n",
|
|
" plt.imshow(E, origin='lower')\n",
|
|
" plt.colorbar(orientation='horizontal', pad=0.02, shrink=0.9)"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"def run_simulation(outfile, t_end, sim_args):\n",
|
|
" with Timer(\"construct\") as t:\n",
|
|
" sim = EE2DKP07Dimsplit(**sim_args)\n",
|
|
" print(\"Constructed in \" + str(t.secs) + \" seconds\")\n",
|
|
"\n",
|
|
" #Create a netcdf file and simulate\n",
|
|
" with DataDumper(outfile, mode='w', clobber=False) as outdata:\n",
|
|
" outdata.nc.createDimension('time', None)\n",
|
|
" outdata.nc.createDimension('x', sim.nx)\n",
|
|
" outdata.nc.createDimension('y', sim.ny)\n",
|
|
"\n",
|
|
" #Create variables\n",
|
|
" outdata.time_var = outdata.nc.createVariable('time', np.dtype('float32').char, 'time')\n",
|
|
" outdata.x_var = outdata.nc.createVariable('x', np.dtype('float32').char, 'x')\n",
|
|
" outdata.y_var = outdata.nc.createVariable('y', np.dtype('float32').char, 'y')\n",
|
|
" outdata.rho_var = outdata.nc.createVariable('rho', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True,\n",
|
|
" least_significant_digit=3)\n",
|
|
" outdata.rho_u_var = outdata.nc.createVariable('rho_u', np.dtype('float32').char, ('time', 'y', 'x'),\n",
|
|
" zlib=True, least_significant_digit=3)\n",
|
|
" outdata.rho_v_var = outdata.nc.createVariable('rho_v', np.dtype('float32').char, ('time', 'y', 'x'),\n",
|
|
" zlib=True, least_significant_digit=3)\n",
|
|
" outdata.E_var = outdata.nc.createVariable('E', np.dtype('float32').char, ('time', 'y', 'x'), zlib=True,\n",
|
|
" least_significant_digit=3)\n",
|
|
"\n",
|
|
" #Create attributes\n",
|
|
" def to_json(in_dict):\n",
|
|
" out_dict = in_dict.copy()\n",
|
|
"\n",
|
|
" for key in out_dict:\n",
|
|
" if isinstance(out_dict[key], np.ndarray):\n",
|
|
" out_dict[key] = out_dict[key].tolist()\n",
|
|
" else:\n",
|
|
" try:\n",
|
|
" json.dumps(out_dict[key])\n",
|
|
" except:\n",
|
|
" out_dict[key] = str(out_dict[key])\n",
|
|
"\n",
|
|
" return json.dumps(out_dict)\n",
|
|
"\n",
|
|
" outdata.nc.created = time.ctime(time.time())\n",
|
|
" outdata.nc.sim_args = to_json(sim_args)\n",
|
|
"\n",
|
|
" outdata.x_var[:] = np.linspace(0, sim.nx * sim.dx, sim.nx)\n",
|
|
" outdata.y_var[:] = np.linspace(0, sim.ny * sim.dy, sim.ny)\n",
|
|
"\n",
|
|
" progress_printer = ProgressPrinter(n_save, print_every=10)\n",
|
|
" print(f\"Simulating to t={t_end}. Estimated {int(t_end / sim.dt)} timesteps (dt={sim.dt})\")\n",
|
|
" for i in range(n_save + 1):\n",
|
|
" #Sanity check simulator\n",
|
|
" try:\n",
|
|
" sim.check()\n",
|
|
" except AssertionError as e:\n",
|
|
" print(f\"Error after {sim.sim_steps()} steps (t={sim.sim_time()}: {str(e)}\")\n",
|
|
" return outdata.filename\n",
|
|
"\n",
|
|
" #Simulate\n",
|
|
" if i > 0:\n",
|
|
" sim.simulate(t_end / n_save)\n",
|
|
"\n",
|
|
" #Save to file\n",
|
|
" #rho = sim.u0[0].download(sim.stream)\n",
|
|
" rho, rho_u, rho_v, E = sim.download()\n",
|
|
" outdata.time_var[i] = sim.sim_time()\n",
|
|
" outdata.rho_var[i, :] = rho\n",
|
|
" outdata.rho_u_var[i, :] = rho_u\n",
|
|
" outdata.rho_v_var[i, :] = rho_v\n",
|
|
" outdata.E_var[i, :] = E\n",
|
|
"\n",
|
|
" #Write progress to screen\n",
|
|
" print_string = progress_printer.getPrintString(i)\n",
|
|
" if print_string:\n",
|
|
" print(print_string)\n",
|
|
"\n",
|
|
" return outdata.filename"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"class VisType(Enum):\n",
|
|
" Schlieren = 0\n",
|
|
" Density = 1\n",
|
|
"\n",
|
|
"\n",
|
|
"def create_animation(infile, vis_type, vmin, vmax, save_anim=False, cmap=plt.cm.coolwarm, fig_args=None):\n",
|
|
" fig = plt.figure(**fig_args)\n",
|
|
"\n",
|
|
" with DataDumper(infile, 'r') as indata:\n",
|
|
" time = indata.nc.variables['time'][:]\n",
|
|
" x = indata.nc.variables['x'][:]\n",
|
|
" y = indata.nc.variables['y'][:]\n",
|
|
" rho = indata.nc.variables['rho'][0]\n",
|
|
" rho_u = indata.nc.variables['rho_u'][0]\n",
|
|
" rho_v = indata.nc.variables['rho_v'][0]\n",
|
|
"\n",
|
|
" created = indata.nc.created\n",
|
|
" sim_args = json.loads(indata.nc.sim_args)\n",
|
|
" for key in sim_args:\n",
|
|
" if isinstance(sim_args[key], list):\n",
|
|
" sim_args[key] = \"[...]\"\n",
|
|
" num_frames = len(time)\n",
|
|
" print(f\"{infile} created {created} contains {num_frames} timesteps\")\n",
|
|
" print(\"Simulator arguments: \\n\", sim_args)\n",
|
|
"\n",
|
|
" ax1 = fig.gca()\n",
|
|
" extents = [0, x.max(), 0, y.max()]\n",
|
|
"\n",
|
|
" if vis_type == VisType.Schlieren:\n",
|
|
" im = ax1.imshow(Visualization.gen_colors(rho, rho_u, rho_v, cmap, vmax, vmin), origin='lower',\n",
|
|
" extent=extents, cmap='gray', vmin=0.0, vmax=1.0)\n",
|
|
" fig.suptitle(\"Schlieren / vorticity at t={:.2f}\".format(time[0]))\n",
|
|
" elif vis_type == VisType.Density:\n",
|
|
" im = ax1.imshow(rho, origin='lower', extent=extents, cmap=cmap, vmin=vmin, vmax=vmax)\n",
|
|
" fig.suptitle(\"Density at t={:.2f}\".format(time[0]))\n",
|
|
" else:\n",
|
|
" assert False, \"Wrong vis_type\"\n",
|
|
"\n",
|
|
" #Create colorbar\n",
|
|
" from matplotlib.colors import Normalize\n",
|
|
" from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
|
|
" divider = make_axes_locatable(ax1)\n",
|
|
" ax2 = divider.append_axes(\"right\", size=0.1, pad=0.05)\n",
|
|
" cb1 = ColorbarBase(ax2, cmap=cmap,\n",
|
|
" norm=Normalize(vmin=vmin, vmax=vmax),\n",
|
|
" orientation='vertical')\n",
|
|
"\n",
|
|
" #Label colorbar\n",
|
|
" if vis_type == VisType.Schlieren:\n",
|
|
" cb1.set_label('Vorticity')\n",
|
|
" elif vis_type == VisType.Density:\n",
|
|
" cb1.set_label('Density')\n",
|
|
"\n",
|
|
" progress_printer = ProgressPrinter(num_frames, print_every=5)\n",
|
|
"\n",
|
|
" def animate(i):\n",
|
|
" rho = indata.nc.variables['rho'][i]\n",
|
|
" rho_u = indata.nc.variables['rho_u'][i]\n",
|
|
" rho_v = indata.nc.variables['rho_v'][i]\n",
|
|
"\n",
|
|
" if vis_type == VisType.Schlieren:\n",
|
|
" im.set_data(Visualization.gen_colors(rho, rho_u, rho_v, cmap, vmax, vmin))\n",
|
|
" fig.suptitle(\"Schlieren / vorticity at t={:.2f}\".format(time[i]))\n",
|
|
" elif vis_type == VisType.Density:\n",
|
|
" im.set_data(rho)\n",
|
|
" fig.suptitle(\"Density at t={:.2f}\".format(time[i]))\n",
|
|
"\n",
|
|
" #Write progress to screen\n",
|
|
" print_string = progress_printer.getPrintString(i)\n",
|
|
" if print_string:\n",
|
|
" print(print_string)\n",
|
|
"\n",
|
|
" anim = animation.FuncAnimation(fig, animate, interval=50, frames=range(num_frames))\n",
|
|
" plt.close()\n",
|
|
"\n",
|
|
" if (save_anim):\n",
|
|
" root, _ = os.path.splitext(infile)\n",
|
|
" movie_outpath = os.path.abspath(root + \".mp4\")\n",
|
|
" if (os.path.isfile(movie_outpath)):\n",
|
|
" 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 FFMpegWriter\n",
|
|
" writer = FFMpegWriter(fps=25)\n",
|
|
" anim.save(movie_outpath, writer=writer)\n",
|
|
" display(HTML(\"\"\"\n",
|
|
" <div align=\"middle\">\n",
|
|
" <video width=\"80%\" controls>\n",
|
|
" <source src=\"{:s}\" type=\"video/mp4\">\n",
|
|
" </video>\n",
|
|
" </div>\n",
|
|
" \"\"\".format(movie_outpath)))\n",
|
|
" else:\n",
|
|
" #plt.rcParams[\"animation.html\"] = \"html5\"\n",
|
|
" plt.rcParams[\"animation.html\"] = \"jshtml\"\n",
|
|
" display(anim)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Shock-bubble"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"nx = 400\n",
|
|
"arguments = InitialConditions.gen_shock_bubble(nx, nx // 4, 1.4)\n",
|
|
"plt.figure()\n",
|
|
"plt.imshow(Visualization.gen_schlieren(arguments['rho']), cmap='gray')\n",
|
|
"plt.colorbar(orientation='vertical', aspect=20, pad=0.02, shrink=0.3)"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"nx = 800 #1600\n",
|
|
"ny = nx // 4\n",
|
|
"g = 0.0\n",
|
|
"gamma = 1.4\n",
|
|
"t_end = 0.5 #3.0\n",
|
|
"n_save = 20 #500\n",
|
|
"outfile = \"data/euler_shock-bubble.nc\"\n",
|
|
"outdata = None\n",
|
|
"\n",
|
|
"arguments = InitialConditions.gen_shock_bubble(nx, ny, gamma)\n",
|
|
"arguments['context'] = my_context\n",
|
|
"outfile = run_simulation(outfile, t_end, arguments)"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"#outfile = 'data/euler_shock-bubble_0008.nc'\n",
|
|
"create_animation(outfile, vis_type=VisType.Schlieren, vmin=-0.2, vmax=0.2, cmap=plt.cm.RdBu, save_anim=False,\n",
|
|
" fig_args={'figsize': (16, 5)})"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Kelvin-Helmholtz"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"nx = 400\n",
|
|
"arguments = InitialConditions.gen_kelvin_helmholtz(nx, nx // 4, 1.4)\n",
|
|
"\n",
|
|
"plt.figure()\n",
|
|
"plt.imshow(arguments['rho'])\n",
|
|
"plt.colorbar(orientation='vertical', aspect=20, pad=0.02, shrink=0.3)"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"nx = 400\n",
|
|
"ny = nx // 2\n",
|
|
"roughness = 0.125\n",
|
|
"t_end = 10.0\n",
|
|
"n_save = 100 #1000\n",
|
|
"outfile = \"data/euler_kelvin_helmholtz.nc\"\n",
|
|
"outdata = None\n",
|
|
"\n",
|
|
"arguments = InitialConditions.gen_kelvin_helmholtz(nx, ny, gamma, roughness)\n",
|
|
"arguments['context'] = my_context\n",
|
|
"outfile = run_simulation(outfile, t_end, arguments)"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"#outfile='data/euler_kelvin_helmholtz_0012.nc'\n",
|
|
"create_animation(outfile, vis_type=VisType.Density, vmin=1.0, vmax=2.0, save_anim=False, fig_args={'figsize': (16, 9)})"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Rayleigh-Taylor\n",
|
|
"\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",
|
|
"* "
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"nx = 400\n",
|
|
"arguments = InitialConditions.gen_rayleigh_taylor(nx, nx * 3, 1.4, version=0)\n",
|
|
"plot_vars(arguments['rho'], arguments['rho_u'], arguments['rho_v'], arguments['E'])"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"nx = 151\n",
|
|
"ny = nx * 3\n",
|
|
"g = 0.1\n",
|
|
"gamma = 1.4\n",
|
|
"t_end = 30\n",
|
|
"n_save = 300\n",
|
|
"outfile = \"data/euler_rayleigh-taylor.nc\"\n",
|
|
"outdata = None\n",
|
|
"\n",
|
|
"arguments = InitialConditions.gen_rayleigh_taylor(nx, ny, gamma)\n",
|
|
"arguments['context'] = my_context\n",
|
|
"outfile = run_simulation(outfile, t_end, arguments)"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": [
|
|
"#outfile = 'data/euler_rayleigh-taylor_0007.nc'\n",
|
|
"create_animation(outfile, vis_type=VisType.Density, vmin=1, vmax=2, cmap=plt.cm.coolwarm, save_anim=False,\n",
|
|
" fig_args={'figsize': (3.4, 8)})"
|
|
]
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"outputs": [],
|
|
"execution_count": null,
|
|
"source": ""
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "ShallowWaterGPU",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.9.19"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|