mirror of
https://github.com/smyalygames/FiniteVolumeGPU.git
synced 2025-10-31 20:17:41 +01:00
refactor(helpers): follow PEP8 formatting standard
This commit is contained in:
Anthony Berg
parent
d4f2ffc493
commit
1469ac1128
@ -62,7 +62,7 @@
|
||||
"#Finally, import our simulator\n",
|
||||
"from GPUSimulators.model import Force, HLL, HLL2, KP07, LxF, WAF, KP07Dimsplit\n",
|
||||
"from GPUSimulators.common import Timer\n",
|
||||
"from GPUSimulators.helpers import InitialConditions"
|
||||
"from GPUSimulators.helpers import initial_conditions"
|
||||
]
|
||||
},
|
||||
{
|
||||
|
||||
@ -62,7 +62,7 @@
|
||||
"#Finally, import our simulator\n",
|
||||
"from GPUSimulators.model import Force, HLL, HLL2, KP07, LxF, WAF, KP07Dimsplit\n",
|
||||
"from GPUSimulators.common import Timer\n",
|
||||
"from GPUSimulators.helpers import InitialConditions"
|
||||
"from GPUSimulators.helpers import initial_conditions"
|
||||
]
|
||||
},
|
||||
{
|
||||
|
||||
@ -32,7 +32,7 @@
|
||||
"\n",
|
||||
"from GPUSimulators.common import Timer, DataDumper, ProgressPrinter\n",
|
||||
"from GPUSimulators.model import EE2DKP07Dimsplit\n",
|
||||
"from GPUSimulators.helpers import InitialConditions, Visualization"
|
||||
"from GPUSimulators.helpers import initial_conditions, visualization"
|
||||
]
|
||||
},
|
||||
{
|
||||
@ -205,7 +205,7 @@
|
||||
" extents = [0, x.max(), 0, y.max()]\n",
|
||||
"\n",
|
||||
" if vis_type == VisType.Schlieren:\n",
|
||||
" im = ax1.imshow(Visualization.genColors(rho, rho_u, rho_v, cmap, vmax, vmin), origin='lower',\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",
|
||||
@ -237,7 +237,7 @@
|
||||
" rho_v = indata.ncfile.variables['rho_v'][i]\n",
|
||||
"\n",
|
||||
" if vis_type == VisType.Schlieren:\n",
|
||||
" im.set_data(Visualization.genColors(rho, rho_u, rho_v, cmap, vmax, vmin))\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",
|
||||
@ -290,9 +290,9 @@
|
||||
"execution_count": null,
|
||||
"source": [
|
||||
"nx = 400\n",
|
||||
"arguments = InitialConditions.genShockBubble(nx, nx // 4, 1.4)\n",
|
||||
"arguments = InitialConditions.gen_shock_bubble(nx, nx // 4, 1.4)\n",
|
||||
"plt.figure()\n",
|
||||
"plt.imshow(Visualization.genSchlieren(arguments['rho']), cmap='gray')\n",
|
||||
"plt.imshow(Visualization.gen_schlieren(arguments['rho']), cmap='gray')\n",
|
||||
"plt.colorbar(orientation='vertical', aspect=20, pad=0.02, shrink=0.3)"
|
||||
]
|
||||
},
|
||||
@ -311,7 +311,7 @@
|
||||
"outfile = \"data/euler_shock-bubble.nc\"\n",
|
||||
"outdata = None\n",
|
||||
"\n",
|
||||
"arguments = InitialConditions.genShockBubble(nx, ny, gamma)\n",
|
||||
"arguments = InitialConditions.gen_shock_bubble(nx, ny, gamma)\n",
|
||||
"arguments['context'] = my_context\n",
|
||||
"outfile = run_simulation(outfile, t_end, arguments)"
|
||||
]
|
||||
@ -341,7 +341,7 @@
|
||||
"execution_count": null,
|
||||
"source": [
|
||||
"nx = 400\n",
|
||||
"arguments = InitialConditions.genKelvinHelmholtz(nx, nx // 4, 1.4)\n",
|
||||
"arguments = InitialConditions.gen_kelvin_helmholtz(nx, nx // 4, 1.4)\n",
|
||||
"\n",
|
||||
"plt.figure()\n",
|
||||
"plt.imshow(arguments['rho'])\n",
|
||||
@ -362,7 +362,7 @@
|
||||
"outfile = \"data/euler_kelvin_helmholtz.nc\"\n",
|
||||
"outdata = None\n",
|
||||
"\n",
|
||||
"arguments = InitialConditions.genKelvinHelmholtz(nx, ny, gamma, roughness)\n",
|
||||
"arguments = InitialConditions.gen_kelvin_helmholtz(nx, ny, gamma, roughness)\n",
|
||||
"arguments['context'] = my_context\n",
|
||||
"outfile = run_simulation(outfile, t_end, arguments)"
|
||||
]
|
||||
@ -395,7 +395,7 @@
|
||||
"execution_count": null,
|
||||
"source": [
|
||||
"nx = 400\n",
|
||||
"arguments = InitialConditions.genRayleighTaylor(nx, nx * 3, 1.4, version=0)\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'])"
|
||||
]
|
||||
},
|
||||
@ -414,7 +414,7 @@
|
||||
"outfile = \"data/euler_rayleigh-taylor.nc\"\n",
|
||||
"outdata = None\n",
|
||||
"\n",
|
||||
"arguments = InitialConditions.genRayleighTaylor(nx, ny, gamma)\n",
|
||||
"arguments = InitialConditions.gen_rayleigh_taylor(nx, ny, gamma)\n",
|
||||
"arguments['context'] = my_context\n",
|
||||
"outfile = run_simulation(outfile, t_end, arguments)"
|
||||
]
|
||||
|
||||
@ -19,12 +19,14 @@ You should have received a copy of the GNU General Public License
|
||||
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
"""
|
||||
|
||||
from GPUSimulators.Simulator import BoundaryCondition
|
||||
import numpy as np
|
||||
import gc
|
||||
|
||||
import numpy as np
|
||||
|
||||
def getExtent(width, height, nx, ny, grid, index=None):
|
||||
from GPUSimulators.Simulator import BoundaryCondition
|
||||
|
||||
|
||||
def get_extent(width, height, nx, ny, grid, index=None):
|
||||
if grid is not None:
|
||||
gx = grid.grid[0]
|
||||
gy = grid.grid[1]
|
||||
@ -32,38 +34,38 @@ def getExtent(width, height, nx, ny, grid, index=None):
|
||||
i, j = grid.get_coordinate(index)
|
||||
else:
|
||||
i, j = grid.get_coordinate()
|
||||
|
||||
|
||||
dx = (width / gx) / nx
|
||||
dy = (height / gy) / ny
|
||||
|
||||
x0 = width*i/gx + 0.5*dx
|
||||
y0 = height*j/gy + 0.5*dy
|
||||
x1 = width*(i+1)/gx - 0.5*dx
|
||||
y1 = height*(j+1)/gy - 0.5*dx
|
||||
|
||||
|
||||
x0 = width * i / gx + 0.5 * dx
|
||||
y0 = height * j / gy + 0.5 * dy
|
||||
x1 = width * (i + 1) / gx - 0.5 * dx
|
||||
y1 = height * (j + 1) / gy - 0.5 * dx
|
||||
|
||||
else:
|
||||
dx = width / nx
|
||||
dy = height / ny
|
||||
|
||||
x0 = 0.5*dx
|
||||
y0 = 0.5*dy
|
||||
x1 = width-0.5*dx
|
||||
y1 = height-0.5*dy
|
||||
|
||||
|
||||
x0 = 0.5 * dx
|
||||
y0 = 0.5 * dy
|
||||
x1 = width - 0.5 * dx
|
||||
y1 = height - 0.5 * dy
|
||||
|
||||
return [x0, x1, y0, y1, dx, dy]
|
||||
|
||||
|
||||
|
||||
def downsample(highres_solution, x_factor, y_factor=None):
|
||||
if (y_factor == None):
|
||||
if y_factor is None:
|
||||
y_factor = x_factor
|
||||
|
||||
assert(highres_solution.shape[1] % x_factor == 0)
|
||||
assert(highres_solution.shape[0] % y_factor == 0)
|
||||
|
||||
if (x_factor*y_factor == 1):
|
||||
assert (highres_solution.shape[1] % x_factor == 0)
|
||||
assert (highres_solution.shape[0] % y_factor == 0)
|
||||
|
||||
if x_factor * y_factor == 1:
|
||||
return highres_solution
|
||||
|
||||
if (len(highres_solution.shape) == 1):
|
||||
|
||||
if len(highres_solution.shape) == 1:
|
||||
highres_solution = highres_solution.reshape((1, highres_solution.size))
|
||||
|
||||
nx = highres_solution.shape[1] / x_factor
|
||||
@ -72,100 +74,96 @@ def downsample(highres_solution, x_factor, y_factor=None):
|
||||
return highres_solution.reshape([int(ny), int(y_factor), int(nx), int(x_factor)]).mean(3).mean(1)
|
||||
|
||||
|
||||
def bump(nx: int, ny: int, width: int, height: int,
|
||||
bump_size=None,
|
||||
ref_nx=None, ref_ny=None,
|
||||
x_center=0.5, y_center=0.5,
|
||||
h_ref=0.5, h_amp=0.1, u_ref=0.0, u_amp=0.1, v_ref=0.0, v_amp=0.1):
|
||||
|
||||
if (ref_nx == None):
|
||||
def bump(nx: int, ny: int, width: int, height: int,
|
||||
bump_size=None,
|
||||
ref_nx=None, ref_ny=None,
|
||||
x_center=0.5, y_center=0.5,
|
||||
h_ref=0.5, h_amp=0.1, u_ref=0.0, u_amp=0.1, v_ref=0.0, v_amp=0.1):
|
||||
if ref_nx is None:
|
||||
ref_nx = nx
|
||||
assert(ref_nx >= nx)
|
||||
|
||||
if (ref_ny == None):
|
||||
assert (ref_nx >= nx)
|
||||
|
||||
if ref_ny is None:
|
||||
ref_ny = ny
|
||||
assert(ref_ny >= ny)
|
||||
|
||||
if (bump_size == None):
|
||||
bump_size = width/5.0
|
||||
|
||||
assert (ref_ny >= ny)
|
||||
|
||||
if bump_size is None:
|
||||
bump_size = width / 5.0
|
||||
|
||||
ref_dx = width / float(ref_nx)
|
||||
ref_dy = height / float(ref_ny)
|
||||
|
||||
x_center = ref_dx*ref_nx*x_center
|
||||
y_center = ref_dy*ref_ny*y_center
|
||||
|
||||
x = ref_dx*(np.arange(0, ref_nx, dtype=np.float32)+0.5) - x_center
|
||||
y = ref_dy*(np.arange(0, ref_ny, dtype=np.float32)+0.5) - y_center
|
||||
x_center = ref_dx * ref_nx * x_center
|
||||
y_center = ref_dy * ref_ny * y_center
|
||||
|
||||
x = ref_dx * (np.arange(0, ref_nx, dtype=np.float32) + 0.5) - x_center
|
||||
y = ref_dy * (np.arange(0, ref_ny, dtype=np.float32) + 0.5) - y_center
|
||||
xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
|
||||
r = np.sqrt(xv**2 + yv**2)
|
||||
r = np.sqrt(xv ** 2 + yv ** 2)
|
||||
xv = None
|
||||
yv = None
|
||||
gc.collect()
|
||||
|
||||
#Generate highres then downsample
|
||||
#h_highres = 0.5 + 0.1*np.exp(-(xv**2/size + yv**2/size))
|
||||
h_highres = h_ref + h_amp*0.5*(1.0 + np.cos(np.pi*r/bump_size)) * (r < bump_size)
|
||||
h = downsample(h_highres, ref_nx/nx, ref_ny/ny)
|
||||
|
||||
# Generate highres then downsample
|
||||
# h_highres = 0.5 + 0.1*np.exp(-(xv**2/size + yv**2/size))
|
||||
h_highres = h_ref + h_amp * 0.5 * (1.0 + np.cos(np.pi * r / bump_size)) * (r < bump_size)
|
||||
h = downsample(h_highres, ref_nx / nx, ref_ny / ny)
|
||||
h_highres = None
|
||||
gc.collect()
|
||||
|
||||
#hu_highres = 0.1*np.exp(-(xv**2/size + yv**2/size))
|
||||
u_highres = u_ref + u_amp*0.5*(1.0 + np.cos(np.pi*r/bump_size)) * (r < bump_size)
|
||||
hu = downsample(u_highres, ref_nx/nx, ref_ny/ny)*h
|
||||
|
||||
# hu_highres = 0.1*np.exp(-(xv**2/size + yv**2/size))
|
||||
u_highres = u_ref + u_amp * 0.5 * (1.0 + np.cos(np.pi * r / bump_size)) * (r < bump_size)
|
||||
hu = downsample(u_highres, ref_nx / nx, ref_ny / ny) * h
|
||||
u_highres = None
|
||||
gc.collect()
|
||||
|
||||
#hu_highres = 0.1*np.exp(-(xv**2/size + yv**2/size))
|
||||
v_highres = v_ref + v_amp*0.5*(1.0 + np.cos(np.pi*r/bump_size)) * (r < bump_size)
|
||||
hv = downsample(v_highres, ref_nx/nx, ref_ny/ny)*h
|
||||
|
||||
# hu_highres = 0.1*np.exp(-(xv**2/size + yv**2/size))
|
||||
v_highres = v_ref + v_amp * 0.5 * (1.0 + np.cos(np.pi * r / bump_size)) * (r < bump_size)
|
||||
hv = downsample(v_highres, ref_nx / nx, ref_ny / ny) * h
|
||||
v_highres = None
|
||||
gc.collect()
|
||||
|
||||
dx = width/nx
|
||||
dy = height/ny
|
||||
|
||||
|
||||
dx = width / nx
|
||||
dy = height / ny
|
||||
|
||||
return h, hu, hv, dx, dy
|
||||
|
||||
|
||||
def genShockBubble(nx, ny, gamma, grid=None):
|
||||
def gen_shock_bubble(nx, ny, gamma, grid=None):
|
||||
"""
|
||||
Generate Shock-bubble interaction case for the Euler equations
|
||||
Generate a Shock-bubble interaction case for the Euler equations
|
||||
"""
|
||||
|
||||
|
||||
width = 4.0
|
||||
height = 1.0
|
||||
g = 0.0
|
||||
|
||||
|
||||
rho = np.ones((ny, nx), dtype=np.float32)
|
||||
u = np.zeros((ny, nx), dtype=np.float32)
|
||||
v = np.zeros((ny, nx), dtype=np.float32)
|
||||
E = np.zeros((ny, nx), dtype=np.float32)
|
||||
p = np.ones((ny, nx), dtype=np.float32)
|
||||
|
||||
|
||||
x0, x1, y0, y1, dx, dy = getExtent(width, height, nx, ny, grid)
|
||||
x0, x1, y0, y1, dx, dy = get_extent(width, height, nx, ny, grid)
|
||||
x = np.linspace(x0, x1, nx, dtype=np.float32)
|
||||
y = np.linspace(y0, y1, ny, dtype=np.float32)
|
||||
xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
|
||||
|
||||
#Bubble
|
||||
|
||||
# Bubble
|
||||
radius = 0.25
|
||||
x_center = 0.5
|
||||
y_center = 0.5
|
||||
bubble = np.sqrt((xv-x_center)**2+(yv-y_center)**2) <= radius
|
||||
bubble = np.sqrt((xv - x_center) ** 2 + (yv - y_center) ** 2) <= radius
|
||||
rho = np.where(bubble, 0.1, rho)
|
||||
|
||||
#Left boundary
|
||||
|
||||
# Left boundary
|
||||
left = (xv < 0.1)
|
||||
rho = np.where(left, 3.81250, rho)
|
||||
u = np.where(left, 2.57669, u)
|
||||
|
||||
#Energy
|
||||
|
||||
# Energy
|
||||
p = np.where(left, 10.0, p)
|
||||
E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)
|
||||
|
||||
E = 0.5 * rho * (u ** 2 + v ** 2) + p / (gamma - 1.0)
|
||||
|
||||
bc = BoundaryCondition({
|
||||
'north': BoundaryCondition.Type.Reflective,
|
||||
@ -173,77 +171,74 @@ def genShockBubble(nx, ny, gamma, grid=None):
|
||||
'east': BoundaryCondition.Type.Periodic,
|
||||
'west': BoundaryCondition.Type.Periodic
|
||||
})
|
||||
|
||||
#Construct simulator
|
||||
|
||||
# Construct simulator
|
||||
arguments = {
|
||||
'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E,
|
||||
'rho': rho, 'rho_u': rho * u, 'rho_v': rho * v, 'E': E,
|
||||
'nx': nx, 'ny': ny,
|
||||
'dx': dx, 'dy': dy,
|
||||
'dx': dx, 'dy': dy,
|
||||
'g': g,
|
||||
'gamma': gamma,
|
||||
'boundary_conditions': bc
|
||||
}
|
||||
}
|
||||
return arguments
|
||||
|
||||
|
||||
def genKelvinHelmholtz(nx, ny, gamma, roughness=0.125, grid=None, index=None):
|
||||
def gen_kelvin_helmholtz(nx, ny, gamma, roughness=0.125, grid=None, index=None):
|
||||
"""
|
||||
Roughness parameter in (0, 1.0] determines how "squiggly"
|
||||
the interface betweeen the zones is
|
||||
the interface between the zones is
|
||||
"""
|
||||
|
||||
def genZones(nx, ny, n):
|
||||
|
||||
def gen_zones(nx, ny, n):
|
||||
"""
|
||||
Generates the zones of the two fluids of K-H
|
||||
"""
|
||||
|
||||
zone = np.zeros((ny, nx), dtype=np.int32)
|
||||
|
||||
|
||||
def genSmoothRandom(nx, n):
|
||||
def gen_smooth_random(nx, n):
|
||||
n = max(1, min(n, nx))
|
||||
|
||||
|
||||
if n == nx:
|
||||
return np.random.random(nx)-0.5
|
||||
return np.random.random(nx) - 0.5
|
||||
else:
|
||||
from scipy.interpolate import interp1d
|
||||
|
||||
#Control points and interpolator
|
||||
# Control points and interpolator
|
||||
xp = np.linspace(0.0, 1.0, n)
|
||||
yp = np.random.random(n) - 0.5
|
||||
|
||||
if (n == 1):
|
||||
|
||||
if n == 1:
|
||||
kind = 'nearest'
|
||||
elif (n == 2):
|
||||
elif n == 2:
|
||||
kind = 'linear'
|
||||
elif (n == 3):
|
||||
elif n == 3:
|
||||
kind = 'quadratic'
|
||||
else:
|
||||
kind = 'cubic'
|
||||
|
||||
|
||||
f = interp1d(xp, yp, kind=kind)
|
||||
|
||||
#Interpolation points
|
||||
# Interpolation points
|
||||
x = np.linspace(0.0, 1.0, nx)
|
||||
return f(x)
|
||||
|
||||
|
||||
|
||||
x0, x1, y0, y1, _, dy = getExtent(1.0, 1.0, nx, ny, grid, index)
|
||||
x0, x1, y0, y1, _, dy = get_extent(1.0, 1.0, nx, ny, grid, index)
|
||||
x = np.linspace(x0, x1, nx)
|
||||
y = np.linspace(y0, y1, ny)
|
||||
_, y = np.meshgrid(x, y)
|
||||
|
||||
#print(y+a[0])
|
||||
# print(y+a[0])
|
||||
|
||||
a = genSmoothRandom(nx, n)*dy
|
||||
zone = np.where(y > 0.25+a, zone, 1)
|
||||
a = gen_smooth_random(nx, n) * dy
|
||||
zone = np.where(y > 0.25 + a, zone, 1)
|
||||
|
||||
a = gen_smooth_random(nx, n) * dy
|
||||
zone = np.where(y < 0.75 + a, zone, 1)
|
||||
|
||||
a = genSmoothRandom(nx, n)*dy
|
||||
zone = np.where(y < 0.75+a, zone, 1)
|
||||
|
||||
return zone
|
||||
|
||||
|
||||
width = 2.0
|
||||
height = 1.0
|
||||
g = 0.0
|
||||
@ -252,49 +247,48 @@ def genKelvinHelmholtz(nx, ny, gamma, roughness=0.125, grid=None, index=None):
|
||||
rho = np.empty((ny, nx), dtype=np.float32)
|
||||
u = np.empty((ny, nx), dtype=np.float32)
|
||||
v = np.zeros((ny, nx), dtype=np.float32)
|
||||
p = 2.5*np.ones((ny, nx), dtype=np.float32)
|
||||
p = 2.5 * np.ones((ny, nx), dtype=np.float32)
|
||||
|
||||
#Generate the different zones
|
||||
zones = genZones(nx, ny, max(1, min(nx, int(nx*roughness))))
|
||||
|
||||
#Zone 0
|
||||
# Generate the different zones
|
||||
zones = gen_zones(nx, ny, max(1, min(nx, int(nx * roughness))))
|
||||
|
||||
# Zone 0
|
||||
zone0 = zones == 0
|
||||
rho = np.where(zone0, 1.0, rho)
|
||||
u = np.where(zone0, 0.5, u)
|
||||
|
||||
#Zone 1
|
||||
|
||||
# Zone 1
|
||||
zone1 = zones == 1
|
||||
rho = np.where(zone1, 2.0, rho)
|
||||
u = np.where(zone1, -0.5, u)
|
||||
|
||||
E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)
|
||||
|
||||
_, _, _, _, dx, dy = getExtent(width, height, nx, ny, grid, index)
|
||||
|
||||
|
||||
|
||||
E = 0.5 * rho * (u ** 2 + v ** 2) + p / (gamma - 1.0)
|
||||
|
||||
_, _, _, _, dx, dy = get_extent(width, height, nx, ny, grid, index)
|
||||
|
||||
bc = BoundaryCondition({
|
||||
'north': BoundaryCondition.Type.Periodic,
|
||||
'south': BoundaryCondition.Type.Periodic,
|
||||
'east': BoundaryCondition.Type.Periodic,
|
||||
'west': BoundaryCondition.Type.Periodic
|
||||
})
|
||||
|
||||
#Construct simulator
|
||||
|
||||
# Construct simulator
|
||||
arguments = {
|
||||
'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E,
|
||||
'rho': rho, 'rho_u': rho * u, 'rho_v': rho * v, 'E': E,
|
||||
'nx': nx, 'ny': ny,
|
||||
'dx': dx, 'dy': dy,
|
||||
'dx': dx, 'dy': dy,
|
||||
'g': g,
|
||||
'gamma': gamma,
|
||||
'boundary_conditions': bc
|
||||
}
|
||||
|
||||
return arguments
|
||||
|
||||
}
|
||||
|
||||
def genRayleighTaylor(nx, ny, gamma, version=0, grid=None):
|
||||
return arguments
|
||||
|
||||
|
||||
def gen_rayleigh_taylor(nx, ny, gamma, version=0, grid=None):
|
||||
"""
|
||||
Generates Rayleigh-Taylor instability case
|
||||
Generates a Rayleigh-Taylor instability case
|
||||
"""
|
||||
|
||||
width = 0.5
|
||||
@ -305,43 +299,42 @@ def genRayleighTaylor(nx, ny, gamma, version=0, grid=None):
|
||||
u = np.zeros((ny, nx), dtype=np.float32)
|
||||
v = np.zeros((ny, nx), dtype=np.float32)
|
||||
p = np.zeros((ny, nx), dtype=np.float32)
|
||||
|
||||
|
||||
x0, x1, y0, y1, dx, dy = getExtent(width, height, nx, ny, grid)
|
||||
x = np.linspace(x0, x1, nx, dtype=np.float32)-width*0.5
|
||||
y = np.linspace(y0, y1, ny, dtype=np.float32)-height*0.5
|
||||
|
||||
x0, x1, y0, y1, dx, dy = get_extent(width, height, nx, ny, grid)
|
||||
x = np.linspace(x0, x1, nx, dtype=np.float32) - width * 0.5
|
||||
y = np.linspace(y0, y1, ny, dtype=np.float32) - height * 0.5
|
||||
xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy')
|
||||
|
||||
#This gives a squigly interfact
|
||||
if (version == 0):
|
||||
y_threshold = 0.01*np.cos(2*np.pi*np.abs(x)/0.5)
|
||||
|
||||
# This gives a squiggly interfaction
|
||||
if version == 0:
|
||||
y_threshold = 0.01 * np.cos(2 * np.pi * np.abs(x) / 0.5)
|
||||
rho = np.where(yv <= y_threshold, 1.0, rho)
|
||||
rho = np.where(yv > y_threshold, 2.0, rho)
|
||||
elif (version == 1):
|
||||
elif version == 1:
|
||||
rho = np.where(yv <= 0.0, 1.0, rho)
|
||||
rho = np.where(yv > 0.0, 2.0, rho)
|
||||
v = 0.01*(1.0 + np.cos(2*np.pi*xv/0.5))/4
|
||||
v = 0.01 * (1.0 + np.cos(2 * np.pi * xv / 0.5)) / 4
|
||||
else:
|
||||
assert False, "Invalid version"
|
||||
|
||||
p = 2.5 - rho*g*yv
|
||||
E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0)
|
||||
|
||||
|
||||
p = 2.5 - rho * g * yv
|
||||
E = 0.5 * rho * (u ** 2 + v ** 2) + p / (gamma - 1.0)
|
||||
|
||||
bc = BoundaryCondition({
|
||||
'north': BoundaryCondition.Type.Reflective,
|
||||
'south': BoundaryCondition.Type.Reflective,
|
||||
'east': BoundaryCondition.Type.Reflective,
|
||||
'west': BoundaryCondition.Type.Reflective
|
||||
})
|
||||
|
||||
#Construct simulator
|
||||
|
||||
# Construct simulator
|
||||
arguments = {
|
||||
'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E,
|
||||
'rho': rho, 'rho_u': rho * u, 'rho_v': rho * v, 'E': E,
|
||||
'nx': nx, 'ny': ny,
|
||||
'dx': dx, 'dy': dy,
|
||||
'dx': dx, 'dy': dy,
|
||||
'g': g,
|
||||
'gamma': gamma,
|
||||
'boundary_conditions': bc
|
||||
}
|
||||
}
|
||||
|
||||
return arguments
|
||||
return arguments
|
||||
@ -20,39 +20,38 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
|
||||
from matplotlib.colors import Normalize
|
||||
|
||||
|
||||
def genSchlieren(rho):
|
||||
#Compute length of z-component of normalized gradient vector
|
||||
normal = np.gradient(rho) #[x, y, 1]
|
||||
length = 1.0 / np.sqrt(normal[0]**2 + normal[1]**2 + 1.0)
|
||||
def gen_schlieren(rho):
|
||||
# Compute length of z-component of normalized gradient vector
|
||||
normal = np.gradient(rho) # [x, y, 1]
|
||||
length = 1.0 / np.sqrt(normal[0] ** 2 + normal[1] ** 2 + 1.0)
|
||||
schlieren = np.power(length, 128)
|
||||
return schlieren
|
||||
|
||||
|
||||
def genVorticity(rho, rho_u, rho_v):
|
||||
def gen_vorticity(rho, rho_u, rho_v):
|
||||
u = rho_u / rho
|
||||
v = rho_v / rho
|
||||
u = np.sqrt(u**2 + v**2)
|
||||
u = np.sqrt(u ** 2 + v ** 2)
|
||||
u_max = u.max()
|
||||
|
||||
|
||||
du_dy, _ = np.gradient(u)
|
||||
_, dv_dx = np.gradient(v)
|
||||
|
||||
#Length of curl
|
||||
|
||||
# Length of curl
|
||||
curl = dv_dx - du_dy
|
||||
return curl
|
||||
|
||||
|
||||
def genColors(rho, rho_u, rho_v, cmap, vmax, vmin):
|
||||
schlieren = genSchlieren(rho)
|
||||
curl = genVorticity(rho, rho_u, rho_v)
|
||||
def gen_colors(rho, rho_u, rho_v, cmap, vmax, vmin):
|
||||
schlieren = gen_schlieren(rho)
|
||||
curl = gen_vorticity(rho, rho_u, rho_v)
|
||||
|
||||
colors = Normalize(vmin, vmax, clip=True)(curl)
|
||||
colors = cmap(colors)
|
||||
for k in range(3):
|
||||
colors[:,:,k] = colors[:,:,k]*schlieren
|
||||
colors[:, :, k] = colors[:, :, k] * schlieren
|
||||
|
||||
return colors
|
||||
return colors
|
||||
@ -177,7 +177,7 @@
|
||||
"%%px\n",
|
||||
"\n",
|
||||
"from GPUSimulators.model import EE2DKP07Dimsplit\n",
|
||||
"from GPUSimulators.helpers import InitialConditions\n",
|
||||
"from GPUSimulators.helpers import initial_conditions\n",
|
||||
"\n",
|
||||
"my_context.autotuner = None\n",
|
||||
"\n",
|
||||
@ -192,7 +192,7 @@
|
||||
"print(grid.get_local_rank())\n",
|
||||
"\n",
|
||||
"#arguments = InitialConditions.genShockBubble(nx, ny, gamma, grid=grid)\n",
|
||||
"arguments = InitialConditions.genKelvinHelmholtz(nx, ny, gamma, grid=grid)\n",
|
||||
"arguments = InitialConditions.gen_kelvin_helmholtz(nx, ny, gamma, grid=grid)\n",
|
||||
"#arguments = InitialConditions.genRayleighTaylor(nx, ny, gamma, grid=grid)\n",
|
||||
"\n",
|
||||
"arguments['context'] = my_context\n",
|
||||
@ -260,7 +260,7 @@
|
||||
"%%px\n",
|
||||
"\n",
|
||||
"from GPUSimulators.model import HLL2\n",
|
||||
"from GPUSimulators.helpers import InitialConditions\n",
|
||||
"from GPUSimulators.helpers import initial_conditions\n",
|
||||
"from GPUSimulators.Simulator import BoundaryCondition\n",
|
||||
"\n",
|
||||
"nx = 128\n",
|
||||
|
||||
@ -29,7 +29,7 @@
|
||||
"rc('animation', html='html5')\n",
|
||||
"\n",
|
||||
"from GPUSimulators.common import Timer\n",
|
||||
"from GPUSimulators.helpers import InitialConditions"
|
||||
"from GPUSimulators.helpers import initial_conditions"
|
||||
]
|
||||
},
|
||||
{
|
||||
|
||||
@ -37,7 +37,7 @@ from GPUSimulators import MPISimulator
|
||||
from GPUSimulators.common import run_simulation, get_git_hash, get_git_status
|
||||
from GPUSimulators.gpu import CudaContext
|
||||
from GPUSimulators.model import EE2DKP07Dimsplit
|
||||
from GPUSimulators.helpers import InitialConditions as IC
|
||||
from GPUSimulators.helpers import initial_conditions as IC
|
||||
|
||||
import argparse
|
||||
|
||||
@ -117,7 +117,7 @@ save_times = np.linspace(0, 0.0000999, 2)
|
||||
outfile = "mpi_out_" + str(MPI.COMM_WORLD.rank) + ".nc"
|
||||
save_var_names = ['rho', 'rho_u', 'rho_v', 'E']
|
||||
|
||||
arguments = IC.genKelvinHelmholtz(nx, ny, gamma, grid=grid)
|
||||
arguments = IC.gen_kelvin_helmholtz(nx, ny, gamma, grid=grid)
|
||||
arguments['context'] = cuda_context
|
||||
arguments['theta'] = 1.2
|
||||
arguments['grid'] = grid
|
||||
|
||||
@ -27,7 +27,7 @@ import logging
|
||||
from GPUSimulators import SHMEMSimulatorGroup
|
||||
from GPUSimulators.common import run_simulation
|
||||
from GPUSimulators.model import EE2DKP07Dimsplit
|
||||
from GPUSimulators.helpers import InitialConditions as IC
|
||||
from GPUSimulators.helpers import initial_conditions as IC
|
||||
|
||||
####
|
||||
# Initialize logging
|
||||
@ -79,7 +79,7 @@ logger.info("Running simulation")
|
||||
|
||||
sims = []
|
||||
for i in range(grid.ngpus):
|
||||
arguments = IC.genKelvinHelmholtz(nx, ny, gamma, grid=grid, index=i)
|
||||
arguments = IC.gen_kelvin_helmholtz(nx, ny, gamma, grid=grid, index=i)
|
||||
arguments['context'] = grid.cuda_contexts[i]
|
||||
arguments['theta'] = 1.2
|
||||
|
||||
|
||||
@ -30,7 +30,7 @@ import pycuda.driver as cuda
|
||||
from GPUSimulators.common import run_simulation
|
||||
from GPUSimulators.gpu import CudaContext
|
||||
from GPUSimulators.model import EE2DKP07Dimsplit
|
||||
from GPUSimulators.helpers import InitialConditions as IC
|
||||
from GPUSimulators.helpers import initial_conditions as IC
|
||||
|
||||
import argparse
|
||||
|
||||
@ -82,7 +82,7 @@ save_times = np.linspace(0, 0.5, 10)
|
||||
outfile = "single_gpu_out.nc"
|
||||
save_var_names = ['rho', 'rho_u', 'rho_v', 'E']
|
||||
|
||||
arguments = IC.genKelvinHelmholtz(nx, ny, gamma)
|
||||
arguments = IC.gen_kelvin_helmholtz(nx, ny, gamma)
|
||||
arguments['context'] = cuda_context
|
||||
arguments['theta'] = 1.2
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user