From 77dc93fd3c1ba81a469f208d08b4ddd5366125e7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9=20R=2E=20Brodtkorb?= Date: Thu, 8 Nov 2018 22:13:37 +0100 Subject: [PATCH] Added helper files --- GPUSimulators/helpers/InitialConditions.py | 261 +++++++++++++++++++++ GPUSimulators/helpers/Visualization.py | 61 +++++ 2 files changed, 322 insertions(+) create mode 100644 GPUSimulators/helpers/InitialConditions.py create mode 100644 GPUSimulators/helpers/Visualization.py diff --git a/GPUSimulators/helpers/InitialConditions.py b/GPUSimulators/helpers/InitialConditions.py new file mode 100644 index 0000000..724823c --- /dev/null +++ b/GPUSimulators/helpers/InitialConditions.py @@ -0,0 +1,261 @@ +# -*- coding: utf-8 -*- + +""" +This python module implements Cuda context handling + +Copyright (C) 2018 SINTEF ICT + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . +""" + + +from GPUSimulators.Simulator import BoundaryCondition +import numpy as np + +def genShockBubble(nx, ny, gamma): + """ + Generate Shock-bubble interaction case for the Euler equations + """ + + width = 4.0 + height = 1.0 + dx = width / float(nx) + dy = height / float(ny) + 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) + + x_center = 0.5 + y_center = 0.5 + x = np.linspace(0.5*dx, nx*dx-0.5*dx, nx, dtype=np.float32) - x_center + y = np.linspace(0.5*dy, ny*dy-0.5*dy, ny, dtype=np.float32) - y_center + xv, yv = np.meshgrid(x, y, sparse=False, indexing='xy') + + #Bubble + radius = 0.25 + bubble = np.sqrt(xv**2+yv**2) <= radius + rho = np.where(bubble, 0.1, rho) + + #Left boundary + left = (xv - xv.min() < 0.1) + rho = np.where(left, 3.81250, rho) + u = np.where(left, 2.57669, u) + + #Energy + p = np.where(left, 10.0, p) + E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0) + + #Estimate dt + scale = 0.45 + max_rho_estimate = 5.0 + max_u_estimate = 5.0 + dx = width/nx + dy = height/ny + dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate)) + + bc = BoundaryCondition({ + 'north': BoundaryCondition.Type.Reflective, + 'south': BoundaryCondition.Type.Reflective, + 'east': BoundaryCondition.Type.Periodic, + 'west': BoundaryCondition.Type.Periodic + }) + + #Construct simulator + arguments = { + 'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E, + 'nx': nx, 'ny': ny, + 'dx': dx, 'dy': dy, 'dt': dt, + 'g': g, + 'gamma': gamma, + 'boundary_conditions': bc + } + return arguments + + + + + + + +def genKelvinHelmholtz(nx, ny, gamma, roughness=0.125): + """ + Roughness parameter in (0, 1.0] determines how "squiggly" + the interface betweeen the zones is + """ + + def genZones(nx, ny, n): + """ + Generates the zones of the two fluids of K-H + """ + zone = np.zeros((ny, nx), dtype=np.int32) + + dx = 1.0 / nx + dy = 1.0 / ny + + def genSmoothRandom(nx, n): + assert (n <= nx), "Number of generated points nx must be larger than n" + + if n == nx: + return np.random.random(nx)-0.5 + else: + from scipy.interpolate import interp1d + + #Control points and interpolator + xp = np.linspace(0.0, 1.0, n) + yp = np.random.random(n) - 0.5 + f = interp1d(xp, yp, kind='cubic') + + #Interpolation points + x = np.linspace(0.0, 1.0, nx) + return f(x) + + + + x = np.linspace(0, 1, nx) + y = np.linspace(0, 1, ny) + + _, y = np.meshgrid(x, y) + + #print(y+a[0]) + + a = genSmoothRandom(nx, n)*dy + zone = np.where(y > 0.25+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 + dx = width / float(nx) + dy = height / float(ny) + g = 0.0 + gamma = 1.4 + + 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) + + #Generate the different zones + zones = genZones(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 + 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) + + #Estimate dt + scale = 0.9 + max_rho_estimate = 3.0 + max_u_estimate = 2.0 + dx = width/nx + dy = height/ny + dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate)) + + + bc = BoundaryCondition({ + 'north': BoundaryCondition.Type.Periodic, + 'south': BoundaryCondition.Type.Periodic, + 'east': BoundaryCondition.Type.Periodic, + 'west': BoundaryCondition.Type.Periodic + }) + + #Construct simulator + arguments = { + 'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E, + 'nx': nx, 'ny': ny, + 'dx': dx, 'dy': dy, 'dt': dt, + 'g': g, + 'gamma': gamma, + 'boundary_conditions': bc + } + + return arguments + + + +def genRayleighTaylor(nx, ny, gamma, version=0): + """ + Generates Rayleigh-Taylor instability case + """ + width = 0.5 + height = 1.5 + dx = width / float(nx) + dy = height / float(ny) + g = 0.1 + + rho = np.zeros((ny, nx), dtype=np.float32) + u = np.zeros((ny, nx), dtype=np.float32) + v = np.zeros((ny, nx), dtype=np.float32) + p = np.zeros((ny, nx), dtype=np.float32) + + x = np.linspace(0.5*dx, nx*dx-0.5*dx, nx, dtype=np.float32)-width*0.5 + y = np.linspace(0.5*dy, ny*dy-0.5*dy, 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) + rho = np.where(yv <= y_threshold, 1.0, rho) + rho = np.where(yv > y_threshold, 2.0, rho) + 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 + else: + assert False, "Invalid version" + + p = 2.5 - rho*g*yv + E = 0.5*rho*(u**2+v**2) + p/(gamma-1.0) + + #Estimate dt + scale = 0.9 + max_rho_estimate = 3.0 + max_u_estimate = 1.0 + dx = width/nx + dy = height/ny + dt = scale * min(dx, dy) / (max_u_estimate + np.sqrt(gamma*max_rho_estimate)) + + bc = BoundaryCondition({ + 'north': BoundaryCondition.Type.Reflective, + 'south': BoundaryCondition.Type.Reflective, + 'east': BoundaryCondition.Type.Reflective, + 'west': BoundaryCondition.Type.Reflective + }) + + #Construct simulator + arguments = { + 'rho': rho, 'rho_u': rho*u, 'rho_v': rho*v, 'E': E, + 'nx': nx, 'ny': ny, + 'dx': dx, 'dy': dy, 'dt': dt, + 'g': g, + 'gamma': gamma, + 'boundary_conditions': bc + } + + return arguments \ No newline at end of file diff --git a/GPUSimulators/helpers/Visualization.py b/GPUSimulators/helpers/Visualization.py new file mode 100644 index 0000000..8314952 --- /dev/null +++ b/GPUSimulators/helpers/Visualization.py @@ -0,0 +1,61 @@ +# -*- coding: utf-8 -*- + +""" +This python module implements Cuda context handling + +Copyright (C) 2018 SINTEF ICT + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . +""" + + + +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) + schlieren = np.power(length, 128) + return schlieren + + +def genVorticity(rho, rho_u, rho_v): + u = rho_u / rho + v = rho_v / rho + u = np.sqrt(u**2 + v**2) + u_max = u.max() + + du_dy, _ = np.gradient(u) + _, dv_dx = np.gradient(v) + + #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) + + colors = Normalize(vmin, vmax, clip=True)(curl) + colors = cmap(colors) + for k in range(3): + colors[:,:,k] = colors[:,:,k]*schlieren + + return colors \ No newline at end of file