diff --git a/3D/advection/advection3d.mp4 b/3D/advection/advection3d.mp4 new file mode 100644 index 0000000..6e1d56f Binary files /dev/null and b/3D/advection/advection3d.mp4 differ diff --git a/3D/karman/D3Q19_karman_sphere.py b/3D/karman/D3Q19_karman_sphere.py new file mode 100644 index 0000000..2c58456 --- /dev/null +++ b/3D/karman/D3Q19_karman_sphere.py @@ -0,0 +1,207 @@ +from __future__ import print_function +from __future__ import division +""" +test: True +""" +import pylbm +from six.moves import range +import numpy as np +import sympy as sp +import sys + +X, Y, Z, LA = sp.symbols('X, Y, Z, LA') +rho, qx, qy, qz = sp.symbols('rho, qx, qy, qz', real=True) + +def feq(v, u): + cs2 = sp.Rational(1, 3) + x, y, z = sp.symbols('x, y, z') + vsymb = sp.Matrix([x, y, z]) + w = sp.Matrix([sp.Rational(1,3)] + [sp.Rational(1, 18)]*6 + [sp.Rational(1, 36)]*12) + f = rho + u.dot(vsymb)/cs2 + u.dot(vsymb)**2/(2*cs2**2) - u.norm()**2/(2*cs2) + return sp.Matrix([w[iv]*f.subs([(x, vv[0]), (y, vv[1]), (z, vv[2])]) for iv, vv in enumerate(v)]) + +def bc_rect(f, m, x, y, z, rhoo, uo): + m[rho] = 0. + m[qx] = rhoo*uo + m[qy] = 0. + m[qz] = 0. + +def plot_vorticity(sol, bornes = False): + #### vorticity + ux = sol.m[qx][:,:,3] + uy = sol.m[qy][:,:,3] + vort = np.abs(ux[1:-1, 2:] - ux[1:-1, :-2] + - uy[2:, 1:-1] + uy[:-2, 1:-1]) + if bornes: + return vort.T, 0.0, 0.1, 1 + else: + return vort.T + +def save(sol, im): + x, y, z = sol.domain.x, sol.domain.y, sol.domain.z + h5 = pylbm.H5File(sol.mpi_topo, 'karman', './karman', im) + h5.set_grid(x, y, z) + h5.add_scalar('rho', sol.m[rho]) + qx_n, qy_n, qz_n = sol.m[qx], sol.m[qy], sol.m[qz] + h5.add_vector('velocity', [qx_n, qy_n, qz_n]) + h5.save() + +def printProgress (iteration, total, prefix = '', suffix = '', decimals = 1, barLength = 100): + """ + Call in a loop to create terminal progress bar + @params: + iteration - Required : current iteration (Int) + total - Required : total iterations (Int) + prefix - Optional : prefix string (Str) + suffix - Optional : suffix string (Str) + decimals - Optional : positive number of decimals in percent complete (Int) + barLength - Optional : character length of bar (Int) + """ + formatStr = '{0:.' + str(decimals) + 'f}' + percents = formatStr.format(100 * (iteration / float(total))) + filledLength = int(round(barLength * iteration / float(total))) + bar = '-' * filledLength + ' ' * (barLength - filledLength) + print('\r{0:s} |{1:s}| {2:s}% {3:s}'.format(prefix, bar, percents, suffix), end='', file=sys.stdout, flush=True) + if iteration == total: + print('', end = '\n', file=sys.stdout, flush=True) + +def run(dx, Tf, generator="cython", sorder=None, withPlot=True): + """ + Parameters + ---------- + + dx: double + spatial step + + Tf: double + final time + + generator: pylbm generator + + sorder: list + storage order + + withPlot: boolean + if True plot the solution otherwise just compute the solution + + """ + la = 1 + rhoo = 1. + uo = 0.1 + radius = 0.125 + + Re = 2000 + nu = rhoo*uo*2*radius/Re + + #tau = .5*(6*nu/la/dx + 1) + #print(1./tau) + + s1 = 1.19 + s2 = s10 = 1.4 + s4 = 1.2 + dummy = 3.0/(la*rhoo*dx) + s9 = 1./(nu*dummy +.5) + s13 = 1./(nu*dummy +.5) + s16 = 1.98 + #[0, s1, s2, 0, s4, 0, s4, 0, s4, s9, s10, s9, s10, s13, s13, s13, s16, s16, s16] + s = [0]*4 + [s1, s9, s9, s13, s13, s13, s4, s4, s4, s16, s16, s16, s10, s10, s2] + r = X**2+Y**2+Z**2 + + d_p = { + 'geometry': { + 'xmin': 0, + 'xmax': 2, + 'ymin': 0, + 'ymax': 1, + 'zmin': 0, + 'zmax': 1 + } + } + + dico = { + 'box':{'x':[0., 2.], 'y':[0., 1.], 'z': [0, 1], 'label':[0, 1, 0, 0, 0, 0]}, + 'elements':[pylbm.Sphere((.3,.5+2*dx,.5+2*dx), radius, 2)], + 'space_step':dx, + 'scheme_velocity':la, + 'schemes':[{ + 'velocities': list(range(19)), + 'conserved_moments':[rho, qx, qy, qz], + 'polynomials':[ + 1, + X, Y, Z, + 19*r - 30, + 3*X**2 - r, + Y**2-Z**2, + X*Y, + Y*Z, + Z*X, + X*(5*r - 9), + Y*(5*r - 9), + Z*(5*r - 9), + X*(Y**2 - Z**2), + Y*(Z**2 - X**2), + Z*(X**2 - Y**2), + (-2*X**2 + Y**2 + Z**2)*(3*r - 5), + -5*Y**2 + 5*Z**2 +3*X**2*(Y**2 - Z**2) + 3*Y**4 -3*Z**4, + -53*r + 21*r**2 + 24 + ], + 'relaxation_parameters': s,#[0]*4 + [1./tau]*15, + 'feq':(feq, (sp.Matrix([qx, qy, qz]),)), + 'init':{ + rho:rhoo, + qx: rhoo*uo, + qy: 0., + qz: 0. + }, + }], + 'boundary_conditions':{ + 0:{'method':{0: pylbm.bc.Bouzidi_bounce_back}, 'value':(bc_rect, (rhoo, uo))}, + 1:{'method':{0: pylbm.bc.Neumann_x}}, + 2:{'method':{0: pylbm.bc.Bouzidi_bounce_back}}, + }, + 'parameters': {LA: la}, + 'generator': generator, + } + + sol = pylbm.Simulation(dico, sorder=sorder) + dt = 1./4 + + if withPlot: + #### choice of the plotted field + plot_field = plot_vorticity + + #### init viewer + viewer = pylbm.viewer.matplotlibViewer + fig = viewer.Fig() + ax = fig[0] + ax.xaxis_set_visible(False) + ax.yaxis_set_visible(False) + field, ymin, ymax, decalh = plot_field(sol, bornes = True) + image = ax.image(field, clim=[ymin, ymax], cmap="jet") + + def update(iframe): + while sol.t < iframe * dt: + sol.one_time_step() + image.set_data(plot_field(sol)) + ax.title = "Solution t={0:f}".format(sol.t) + + + #### run the simulation + fig.animate(update, interval=1) + fig.show() + else: + im = 0 + save(sol, im) + while sol.t < Tf: + im += 1 + while sol.t < im * dt: + sol.one_time_step() + #printProgress(im, int(Tf/dt), prefix = 'Progress:', suffix = 'Complete', barLength = 50) + save(sol, im) + + return sol + +if __name__ == '__main__': + dx = 1./256 + Tf = 100. + run(dx, Tf, withPlot=False) diff --git a/3D/karman/karman3d.mp4 b/3D/karman/karman3d.mp4 new file mode 100644 index 0000000..ff6f004 Binary files /dev/null and b/3D/karman/karman3d.mp4 differ diff --git a/3D/rayleigh-benard/rayleigh-benard-3d.mp4 b/3D/rayleigh-benard/rayleigh-benard-3d.mp4 new file mode 100644 index 0000000..9195a5a Binary files /dev/null and b/3D/rayleigh-benard/rayleigh-benard-3d.mp4 differ diff --git a/3D/rayleigh-benard/rayleigh-benard-3d.ogv b/3D/rayleigh-benard/rayleigh-benard-3d.ogv new file mode 100644 index 0000000..b9694dd Binary files /dev/null and b/3D/rayleigh-benard/rayleigh-benard-3d.ogv differ diff --git a/3D/rayleigh-benard/rayleigh-benard.py b/3D/rayleigh-benard/rayleigh-benard.py new file mode 100644 index 0000000..d8718d3 --- /dev/null +++ b/3D/rayleigh-benard/rayleigh-benard.py @@ -0,0 +1,172 @@ +from __future__ import print_function +from __future__ import division +""" +test: True +""" +from six.moves import range +import numpy as np +import sympy as sp +import mpi4py.MPI as mpi +import pyLBM + +X, Y, Z = sp.symbols('X, Y, Z') +rho, qx, qy, qz, T, LA = sp.symbols('rho, qx, qy, qz, T, LA', real=True) + +# parameters +dx = 1./128 +la = 1 +cs = la/np.sqrt(3) + +Tu = -0.5 +Td = 0.5 + +Ra = 1e6 +Pr = 0.71 +g = 9.81 +tau = 1./1.8 +nu = (2*tau-1)/6*la*dx + +diffusivity = nu/Pr +taup = .5+2*diffusivity/la/dx + +DeltaT = Td - Tu +xmin, xmax, ymin, ymax, zmin, zmax = 0., 2., 0., 1., 0., 2. +H = ymax - ymin +beta = Ra*nu*diffusivity/(g*DeltaT*H**3) + +sf = [0]*4 + [1./tau]*15 +sT = [0] + [1./taup]*5 + +def init_T(x, y, z): + #ones = np.ones((x.size, y.size, z.size)) + #T = ( Tu*ones*(y>=.1) + # + Td*ones*(y<.1) + # + 2*Td*ones*(y>=.1)*np.logical_and(y<.25, ((x-1)**2+(z-1)**2)<.1**2) + # ) + #return T + return Td + (Tu-Td)/(ymax-ymin)*(y-ymin) + (Td-Tu) * (0.1*np.random.random_sample((x.shape[0],y.shape[1],z.shape[2]))-0.5) + +def bc_up(f, m, x, y, z): + m[qx] = 0. + m[qy] = 0. + m[qz] = 0. + m[T] = Tu + +def bc_down(f, m, x, y, z): + m[qx] = 0. + m[qy] = 0. + m[qz] = 0. + m[T] = Td + +def save(sol, im): + x, y, z = sol.domain.x, sol.domain.y, sol.domain.z + h5 = pyLBM.H5File(sol.mpi_topo, 'rayleigh-benard', './rayleigh-benard', im) + h5.set_grid(x, y, z) + h5.add_scalar('T', sol.m[T]) + h5.save() + +def feq_NS(v, u): + cs2 = sp.Rational(1, 3) + x, y, z = sp.symbols('x, y, z') + vsymb = sp.Matrix([x, y, z]) + w = sp.Matrix([sp.Rational(1, 3)] + [sp.Rational(1, 18)]*6 + [sp.Rational(1, 36)]*12) + f = rho + u.dot(vsymb)/cs2 + u.dot(vsymb)**2/(2*cs2**2) - u.norm()**2/(2*cs2) + return sp.Matrix([w[iv]*f.subs([(x, vv[0]), (y, vv[1]), (z, vv[2])]) for iv, vv in enumerate(v)]) + +def feq_T(v, u): + c0 = 1#LA + x, y, z = sp.symbols('x, y, z') + vsymb = sp.Matrix([x, y, z]) + f = T/6*(1 + 2*u.dot(vsymb)/c0) + return sp.Matrix([f.subs([(x, vv[0]), (y, vv[1]), (z, vv[2])]) for iv, vv in enumerate(v)]) + +def run(dx, Tf, generator="cython", sorder=None, withPlot=True): + """ + Parameters + ---------- + + dx: double + spatial step + + Tf: double + final time + + generator: pyLBM generator + + sorder: list + storage order + + withPlot: boolean + if True plot the solution otherwise just compute the solution + + """ + r = X**2+Y**2+Z**2 + + dico = { + 'box':{'x':[xmin, xmax], 'y':[ymin, ymax], 'z':[zmin, zmax], 'label':[-1, -1, 0, 1, -1, -1]}, + 'space_step':dx, + 'scheme_velocity':la, + 'schemes':[ + { + 'velocities':list(range(19)), + 'conserved_moments': [rho, qx, qy, qz], + 'polynomials':[ + 1, + X, Y, Z, + 19*r - 30, + 2*X**2 - Y**2 - Z**2, + Y**2-Z**2, + X*Y, + Y*Z, + Z*X, + X*(5*r - 9), + Y*(5*r - 9), + Z*(5*r - 9), + X*(Y**2 - Z**2), + Y*(Z**2 - X**2), + Z*(X**2 - Y**2), + (2*X**2 - Y**2 - Z**2)*(3*r - 5), + (Y**2 - Z**2)*(3*r - 5), + -sp.Rational(53,2)*r + sp.Rational(21,2)*r**2 + 12 + ], + 'relaxation_parameters':sf, + 'feq':(feq_NS, (sp.Matrix([qx, qy, qz]),)), + 'source_terms':{qy: beta*g*T}, + 'init':{rho: 1., qx: 0., qy: 0., qz: 0.}, + }, + { + 'velocities':list(range(1,7)), + 'conserved_moments': [T], + 'polynomials':[1, X, Y, Z, + X**2 - Y**2, + Y**2 - Z**2, + ], + 'feq':(feq_T, (sp.Matrix([qx, qy, qz]),)), + 'relaxation_parameters':sT, + 'init':{T:(init_T,)}, + }, + ], + 'boundary_conditions':{ + 0:{'method':{0: pyLBM.bc.Bouzidi_bounce_back, 1: pyLBM.bc.Bouzidi_anti_bounce_back}, 'value':bc_down}, + 1:{'method':{0: pyLBM.bc.Bouzidi_bounce_back, 1: pyLBM.bc.Bouzidi_anti_bounce_back}, 'value':bc_up}, + }, + 'generator': "cython", + 'parameters': {LA: la}, + } + + sol = pyLBM.Simulation(dico) + + im = 0 + compt = 0 + while sol.t < Tf: + sol.one_time_step() + compt += 1 + if compt == 128: + im += 1 + save(sol, im) + compt = 0 + return sol + +if __name__ == '__main__': + Tf = 400. + run(dx, Tf)