From 24847cbb2567f508a7c30b39c6fb7ba6379d1adc Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Sat, 2 Nov 2019 17:18:32 +0100 Subject: Restructure LBM model / lattice distinction --- boltzgen.py | 11 ++- boltzgen/__init__.py | 2 +- .../kernel/template/bounce_back_boundary.cl.mako | 2 +- .../kernel/template/bounce_back_boundary.cpp.mako | 2 +- .../kernel/template/collide_and_stream.cl.mako | 2 +- .../kernel/template/collide_and_stream.cpp.mako | 2 +- boltzgen/kernel/template/momenta_boundary.cl.mako | 4 +- boltzgen/kernel/template/momenta_boundary.cpp.mako | 2 +- boltzgen/lbm/__init__.py | 80 +--------------------- boltzgen/lbm/lattice/D2Q9.py | 7 ++ boltzgen/lbm/lattice/D3Q19.py | 18 +++++ boltzgen/lbm/lattice/D3Q27.py | 7 ++ boltzgen/lbm/lattice/D3Q7.py | 18 +++++ boltzgen/lbm/lattice/__init__.py | 6 ++ boltzgen/lbm/lattice/characteristics.py | 27 ++++++++ boltzgen/lbm/model/BGK.py | 76 ++++++++++++++++++++ boltzgen/lbm/model/D2Q9.py | 7 -- boltzgen/lbm/model/D3Q19.py | 18 ----- boltzgen/lbm/model/D3Q27.py | 7 -- boltzgen/lbm/model/D3Q7.py | 18 ----- boltzgen/lbm/model/__init__.py | 7 +- boltzgen/lbm/model/characteristics.py | 27 -------- boltzgen/utility/__init__.py | 5 ++ 23 files changed, 183 insertions(+), 172 deletions(-) create mode 100644 boltzgen/lbm/lattice/D2Q9.py create mode 100644 boltzgen/lbm/lattice/D3Q19.py create mode 100644 boltzgen/lbm/lattice/D3Q27.py create mode 100644 boltzgen/lbm/lattice/D3Q7.py create mode 100644 boltzgen/lbm/lattice/__init__.py create mode 100644 boltzgen/lbm/lattice/characteristics.py create mode 100644 boltzgen/lbm/model/BGK.py delete mode 100644 boltzgen/lbm/model/D2Q9.py delete mode 100644 boltzgen/lbm/model/D3Q19.py delete mode 100644 boltzgen/lbm/model/D3Q27.py delete mode 100644 boltzgen/lbm/model/D3Q7.py delete mode 100644 boltzgen/lbm/model/characteristics.py diff --git a/boltzgen.py b/boltzgen.py index 82adafe..4317c68 100755 --- a/boltzgen.py +++ b/boltzgen.py @@ -8,7 +8,8 @@ argparser = argparse.ArgumentParser( argparser.add_argument('target', help = 'Target language (currently either "cl" or "cpp")') -argparser.add_argument('--lattice', required = True, help = 'Lattice type (D2Q9, D3Q7, D3Q19, D3Q27)') +argparser.add_argument('--lattice', required = True, help = 'Lattice type ("D2Q9", "D3Q7", "D3Q19", "D3Q27")') +argparser.add_argument('--model', required = False, help = 'LBM model (currently only "BGK")') argparser.add_argument('--layout', required = True, help = 'Memory layout ("AOS" or "SOA")') argparser.add_argument('--index', required = False, help = 'Cell indexing ("XYZ" or "ZYX")') argparser.add_argument('--precision', required = True, help = 'Floating precision ("single" or "double")') @@ -21,13 +22,17 @@ argparser.add_argument('--extras', action = 'append', nargs = '+', default args = argparser.parse_args() -lattice = eval("lbm.model.%s" % args.lattice) +if args.model is None: + args.model = "BGK" + +lattice = eval("lbm.lattice.%s" % args.lattice) +model = eval("lbm.model.%s" % args.model) if args.index is None: args.index = 'XYZ' generator = Generator( - model = LBM(lattice, tau = float(args.tau), optimize = not args.disable_cse), + model = model(lattice, tau = float(args.tau), optimize = not args.disable_cse), target = args.target, precision = args.precision, index = args.index, diff --git a/boltzgen/__init__.py b/boltzgen/__init__.py index 01fc32c..d319fdc 100644 --- a/boltzgen/__init__.py +++ b/boltzgen/__init__.py @@ -1,5 +1,5 @@ +import boltzgen.lbm.lattice import boltzgen.lbm.model -from boltzgen.lbm import LBM from boltzgen.geometry import Geometry from boltzgen.kernel.generator import Generator diff --git a/boltzgen/kernel/template/bounce_back_boundary.cl.mako b/boltzgen/kernel/template/bounce_back_boundary.cl.mako index 8bb4f6c..4766bba 100644 --- a/boltzgen/kernel/template/bounce_back_boundary.cl.mako +++ b/boltzgen/kernel/template/bounce_back_boundary.cl.mako @@ -10,7 +10,7 @@ __kernel void bounce_back_boundary_gid(__global ${float_type}* f_next, % endfor <% - subexpr, assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = True)) + subexpr, assignment = model.collision(f_eq = model.equilibrium(resolve_moments = True)) %> % for i, expr in enumerate(subexpr): diff --git a/boltzgen/kernel/template/bounce_back_boundary.cpp.mako b/boltzgen/kernel/template/bounce_back_boundary.cpp.mako index 91dcfa0..c37233e 100644 --- a/boltzgen/kernel/template/bounce_back_boundary.cpp.mako +++ b/boltzgen/kernel/template/bounce_back_boundary.cpp.mako @@ -10,7 +10,7 @@ void bounce_back_boundary( ${float_type}* f_next, % endfor <% - subexpr, assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = True)) + subexpr, assignment = model.collision(f_eq = model.equilibrium(resolve_moments = True)) %> % for i, expr in enumerate(subexpr): diff --git a/boltzgen/kernel/template/collide_and_stream.cl.mako b/boltzgen/kernel/template/collide_and_stream.cl.mako index 303b21b..6b55cd2 100644 --- a/boltzgen/kernel/template/collide_and_stream.cl.mako +++ b/boltzgen/kernel/template/collide_and_stream.cl.mako @@ -10,7 +10,7 @@ __kernel void collide_and_stream_gid(__global ${float_type}* f_next, % endfor <% - subexpr, assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = True)) + subexpr, assignment = model.collision(f_eq = model.equilibrium(resolve_moments = True)) %> % for i, expr in enumerate(subexpr): diff --git a/boltzgen/kernel/template/collide_and_stream.cpp.mako b/boltzgen/kernel/template/collide_and_stream.cpp.mako index 47dd6dd..2c52b54 100644 --- a/boltzgen/kernel/template/collide_and_stream.cpp.mako +++ b/boltzgen/kernel/template/collide_and_stream.cpp.mako @@ -10,7 +10,7 @@ void collide_and_stream( ${float_type}* f_next, % endfor <% - subexpr, assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = True)) + subexpr, assignment = model.collision(f_eq = model.equilibrium(resolve_moments = True)) %> % for i, expr in enumerate(subexpr): diff --git a/boltzgen/kernel/template/momenta_boundary.cl.mako b/boltzgen/kernel/template/momenta_boundary.cl.mako index 4e29601..7c1e3df 100644 --- a/boltzgen/kernel/template/momenta_boundary.cl.mako +++ b/boltzgen/kernel/template/momenta_boundary.cl.mako @@ -1,6 +1,6 @@ <% moments_subexpr, moments_assignment = model.moments() -collision_subexpr, collision_assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = False)) +collision_subexpr, collision_assignment = model.collision(f_eq = model.equilibrium(resolve_moments = False)) %> <%def name="momenta_boundary(name, param)"> @@ -30,7 +30,7 @@ __kernel void ${name}_momenta_boundary_gid( const ${float_type} ${ccode(expr)} % endfor -% for i in range(0,descriptor.q): +% for i, expr in enumerate(collision_assignment): preshifted_f_next[${layout.pop_offset(i)}] = f_next_${i}; % endfor } diff --git a/boltzgen/kernel/template/momenta_boundary.cpp.mako b/boltzgen/kernel/template/momenta_boundary.cpp.mako index 6d05288..c41d07e 100644 --- a/boltzgen/kernel/template/momenta_boundary.cpp.mako +++ b/boltzgen/kernel/template/momenta_boundary.cpp.mako @@ -1,6 +1,6 @@ <% moments_subexpr, moments_assignment = model.moments() -collision_subexpr, collision_assignment = model.bgk(f_eq = model.equilibrium(resolve_moments = False)) +collision_subexpr, collision_assignment = model.collision(f_eq = model.equilibrium(resolve_moments = False)) %> <%def name="momenta_boundary(name, param)"> diff --git a/boltzgen/lbm/__init__.py b/boltzgen/lbm/__init__.py index 3f608cb..276c0bb 100644 --- a/boltzgen/lbm/__init__.py +++ b/boltzgen/lbm/__init__.py @@ -1,78 +1,2 @@ -from sympy import * -from sympy.codegen.ast import Assignment - -import boltzgen.utility.optimizations as optimizations -from boltzgen.lbm.model.characteristics import weights, c_s - - -def assign(names, definitions): - return list(map(lambda x: Assignment(*x), zip(names, definitions))) - -class LBM: - def __init__(self, descriptor, tau, optimize = True): - self.descriptor = descriptor - self.tau = tau - self.optimize = optimize - - if self.tau <= 0.5: - raise Exception('Relaxation time must be larger than 0.5') - - self.f_next = symarray('f_next', descriptor.q) - self.f_curr = symarray('f_curr', descriptor.q) - - if not hasattr(descriptor, 'w'): - self.descriptor.w = weights(descriptor.d, descriptor.c) - - if not hasattr(descriptor, 'c_s'): - self.descriptor.c_s = c_s(descriptor.d, descriptor.c, self.descriptor.w) - - def moments(self, optimize = None): - if optimize is None: - optimize = self.optimize - - rho = symbols('rho') - u = Matrix(symarray('u', self.descriptor.d)) - - exprs = [ Assignment(rho, sum(self.f_curr)) ] - - for i, u_i in enumerate(u): - exprs.append( - Assignment(u_i, sum([ (c_j*self.f_curr[j])[i] for j, c_j in enumerate(self.descriptor.c) ]) / sum(self.f_curr))) - - if optimize: - return cse(exprs, optimizations=optimizations.custom, symbols=numbered_symbols(prefix='m')) - else: - return ([], exprs) - - def equilibrium(self, resolve_moments = False): - rho = symbols('rho') - u = Matrix(symarray('u', self.descriptor.d)) - - if resolve_moments: - moments = self.moments(optimize = False)[1] - rho = moments[0].rhs - for i, m in enumerate(moments[1:]): - u[i] = m.rhs - - f_eq = [] - - for i, c_i in enumerate(self.descriptor.c): - f_eq_i = self.descriptor.w[i] * rho * ( 1 - + c_i.dot(u) / self.descriptor.c_s**2 - + c_i.dot(u)**2 / (2*self.descriptor.c_s**4) - - u.dot(u) / (2*self.descriptor.c_s**2) ) - f_eq.append(f_eq_i) - - return f_eq - - def bgk(self, f_eq, optimize = None): - if optimize is None: - optimize = self.optimize - - exprs = [ self.f_curr[i] + 1/self.tau * (f_eq_i - self.f_curr[i]) for i, f_eq_i in enumerate(f_eq) ] - - if optimize: - subexprs, f = cse(exprs, optimizations=optimizations.custom) - return (subexprs, assign(self.f_next, f)) - else: - return ([], assign(self.f_next, exprs)) +from . import model +from . import lattice diff --git a/boltzgen/lbm/lattice/D2Q9.py b/boltzgen/lbm/lattice/D2Q9.py new file mode 100644 index 0000000..e3ac9de --- /dev/null +++ b/boltzgen/lbm/lattice/D2Q9.py @@ -0,0 +1,7 @@ +from sympy import Matrix +from itertools import product + +d = 2 +q = 9 + +c = [ Matrix(x) for x in product([-1,0,1], repeat=d) ] diff --git a/boltzgen/lbm/lattice/D3Q19.py b/boltzgen/lbm/lattice/D3Q19.py new file mode 100644 index 0000000..e9e6eec --- /dev/null +++ b/boltzgen/lbm/lattice/D3Q19.py @@ -0,0 +1,18 @@ +from sympy import Matrix, Rational, sqrt + +d = 3 +q = 19 + +c = [ Matrix(x) for x in [ + ( 0, 1, 1), (-1, 0, 1), ( 0, 0, 1), ( 1, 0, 1), ( 0, -1, 1), + (-1, 1, 0), ( 0, 1, 0), ( 1, 1, 0), (-1, 0, 0), ( 0, 0, 0), ( 1, 0, 0), (-1,-1, 0), ( 0, -1, 0), ( 1, -1, 0), + ( 0, 1,-1), (-1, 0,-1), ( 0, 0,-1), ( 1, 0,-1), ( 0, -1,-1) +]] + +w = [Rational(*x) for x in [ + (1,36), (1,36), (1,18), (1,36), (1,36), + (1,36), (1,18), (1,36), (1,18), (1,3), (1,18), (1,36), (1,18), (1,36), + (1,36), (1,36), (1,18), (1,36), (1,36) +]] + +c_s = sqrt(Rational(1,3)) diff --git a/boltzgen/lbm/lattice/D3Q27.py b/boltzgen/lbm/lattice/D3Q27.py new file mode 100644 index 0000000..6fb1f80 --- /dev/null +++ b/boltzgen/lbm/lattice/D3Q27.py @@ -0,0 +1,7 @@ +from sympy import Matrix +from itertools import product + +d = 3 +q = 27 + +c = [ Matrix(x) for x in product([-1,0,1], repeat=d) ] diff --git a/boltzgen/lbm/lattice/D3Q7.py b/boltzgen/lbm/lattice/D3Q7.py new file mode 100644 index 0000000..04e16a3 --- /dev/null +++ b/boltzgen/lbm/lattice/D3Q7.py @@ -0,0 +1,18 @@ +from sympy import * + +q = 7 +d = 3 + +c = [ Matrix(x) for x in [ + ( 0, 0, 1), + ( 0, 1, 0), (-1, 0, 0), ( 0, 0, 0), ( 1, 0, 0), ( 0, -1, 0), + ( 0, 0,-1) +]] + +w = [Rational(*x) for x in [ + (1,8), + (1,8), (1,8), (1,4), (1,8), (1,8), + (1,8) +]] + +c_s = sqrt(Rational(1,4)) diff --git a/boltzgen/lbm/lattice/__init__.py b/boltzgen/lbm/lattice/__init__.py new file mode 100644 index 0000000..c956759 --- /dev/null +++ b/boltzgen/lbm/lattice/__init__.py @@ -0,0 +1,6 @@ +from . import D2Q9 +from . import D3Q7 +from . import D3Q19 +from . import D3Q27 + +from . import characteristics diff --git a/boltzgen/lbm/lattice/characteristics.py b/boltzgen/lbm/lattice/characteristics.py new file mode 100644 index 0000000..b68afeb --- /dev/null +++ b/boltzgen/lbm/lattice/characteristics.py @@ -0,0 +1,27 @@ +from sympy import * + +# copy of `sympy.integrals.quadrature.gauss_hermite` sans evaluation +def gauss_hermite(n): + x = Dummy("x") + p = hermite_poly(n, x, polys=True) + p1 = hermite_poly(n-1, x, polys=True) + xi = [] + w = [] + for r in p.real_roots(): + xi.append(r) + w.append(((2**(n-1) * factorial(n) * sqrt(pi))/(n**2 * p1.subs(x, r)**2))) + return xi, w + +# determine weights of a d-dimensional LBM model on velocity set c +# (only works for velocity sets that result into NSE-recovering LB models when +# plugged into Gauss-Hermite quadrature without any additional arguments +# i.e. D2Q9 and D3Q27 but not D3Q19) +def weights(d, c): + _, omegas = gauss_hermite(3) + return list(map(lambda c_i: Mul(*[ omegas[1+c_i[iDim]] for iDim in range(0,d) ]) / pi**(d/2), c)) + +# determine lattice speed of sound using directions and their weights +def c_s(d, c, w): + speeds = set([ sqrt(sum([ w[i] * c_i[j]**2 for i, c_i in enumerate(c) ])) for j in range(0,d) ]) + assert len(speeds) == 1 # verify isotropy + return speeds.pop() diff --git a/boltzgen/lbm/model/BGK.py b/boltzgen/lbm/model/BGK.py new file mode 100644 index 0000000..f5a34e2 --- /dev/null +++ b/boltzgen/lbm/model/BGK.py @@ -0,0 +1,76 @@ +from sympy import * +from sympy.codegen.ast import Assignment + +from boltzgen.utility import assign +import boltzgen.utility.optimizations as optimizations + +from boltzgen.lbm.lattice.characteristics import weights, c_s + +class BGK: + def __init__(self, descriptor, tau, optimize = True): + self.descriptor = descriptor + self.tau = tau + self.optimize = optimize + + if self.tau <= 0.5: + raise Exception('Relaxation time must be larger than 0.5') + + self.f_next = symarray('f_next', descriptor.q) + self.f_curr = symarray('f_curr', descriptor.q) + + if not hasattr(descriptor, 'w'): + self.descriptor.w = weights(descriptor.d, descriptor.c) + + if not hasattr(descriptor, 'c_s'): + self.descriptor.c_s = c_s(descriptor.d, descriptor.c, self.descriptor.w) + + def moments(self, optimize = None): + if optimize is None: + optimize = self.optimize + + rho = symbols('rho') + u = Matrix(symarray('u', self.descriptor.d)) + + exprs = [ Assignment(rho, sum(self.f_curr)) ] + + for i, u_i in enumerate(u): + exprs.append( + Assignment(u_i, sum([ (c_j*self.f_curr[j])[i] for j, c_j in enumerate(self.descriptor.c) ]) / sum(self.f_curr))) + + if optimize: + return cse(exprs, optimizations=optimizations.custom, symbols=numbered_symbols(prefix='m')) + else: + return ([], exprs) + + def equilibrium(self, resolve_moments = False): + rho = symbols('rho') + u = Matrix(symarray('u', self.descriptor.d)) + + if resolve_moments: + moments = self.moments(optimize = False)[1] + rho = moments[0].rhs + for i, m in enumerate(moments[1:]): + u[i] = m.rhs + + f_eq = [] + + for i, c_i in enumerate(self.descriptor.c): + f_eq_i = self.descriptor.w[i] * rho * ( 1 + + c_i.dot(u) / self.descriptor.c_s**2 + + c_i.dot(u)**2 / (2*self.descriptor.c_s**4) + - u.dot(u) / (2*self.descriptor.c_s**2) ) + f_eq.append(f_eq_i) + + return f_eq + + def collision(self, f_eq, optimize = None): + if optimize is None: + optimize = self.optimize + + exprs = [ self.f_curr[i] + 1/self.tau * (f_eq_i - self.f_curr[i]) for i, f_eq_i in enumerate(f_eq) ] + + if optimize: + subexprs, f = cse(exprs, optimizations=optimizations.custom) + return (subexprs, assign(self.f_next, f)) + else: + return ([], assign(self.f_next, exprs)) diff --git a/boltzgen/lbm/model/D2Q9.py b/boltzgen/lbm/model/D2Q9.py deleted file mode 100644 index e3ac9de..0000000 --- a/boltzgen/lbm/model/D2Q9.py +++ /dev/null @@ -1,7 +0,0 @@ -from sympy import Matrix -from itertools import product - -d = 2 -q = 9 - -c = [ Matrix(x) for x in product([-1,0,1], repeat=d) ] diff --git a/boltzgen/lbm/model/D3Q19.py b/boltzgen/lbm/model/D3Q19.py deleted file mode 100644 index e9e6eec..0000000 --- a/boltzgen/lbm/model/D3Q19.py +++ /dev/null @@ -1,18 +0,0 @@ -from sympy import Matrix, Rational, sqrt - -d = 3 -q = 19 - -c = [ Matrix(x) for x in [ - ( 0, 1, 1), (-1, 0, 1), ( 0, 0, 1), ( 1, 0, 1), ( 0, -1, 1), - (-1, 1, 0), ( 0, 1, 0), ( 1, 1, 0), (-1, 0, 0), ( 0, 0, 0), ( 1, 0, 0), (-1,-1, 0), ( 0, -1, 0), ( 1, -1, 0), - ( 0, 1,-1), (-1, 0,-1), ( 0, 0,-1), ( 1, 0,-1), ( 0, -1,-1) -]] - -w = [Rational(*x) for x in [ - (1,36), (1,36), (1,18), (1,36), (1,36), - (1,36), (1,18), (1,36), (1,18), (1,3), (1,18), (1,36), (1,18), (1,36), - (1,36), (1,36), (1,18), (1,36), (1,36) -]] - -c_s = sqrt(Rational(1,3)) diff --git a/boltzgen/lbm/model/D3Q27.py b/boltzgen/lbm/model/D3Q27.py deleted file mode 100644 index 6fb1f80..0000000 --- a/boltzgen/lbm/model/D3Q27.py +++ /dev/null @@ -1,7 +0,0 @@ -from sympy import Matrix -from itertools import product - -d = 3 -q = 27 - -c = [ Matrix(x) for x in product([-1,0,1], repeat=d) ] diff --git a/boltzgen/lbm/model/D3Q7.py b/boltzgen/lbm/model/D3Q7.py deleted file mode 100644 index 04e16a3..0000000 --- a/boltzgen/lbm/model/D3Q7.py +++ /dev/null @@ -1,18 +0,0 @@ -from sympy import * - -q = 7 -d = 3 - -c = [ Matrix(x) for x in [ - ( 0, 0, 1), - ( 0, 1, 0), (-1, 0, 0), ( 0, 0, 0), ( 1, 0, 0), ( 0, -1, 0), - ( 0, 0,-1) -]] - -w = [Rational(*x) for x in [ - (1,8), - (1,8), (1,8), (1,4), (1,8), (1,8), - (1,8) -]] - -c_s = sqrt(Rational(1,4)) diff --git a/boltzgen/lbm/model/__init__.py b/boltzgen/lbm/model/__init__.py index c956759..fb8f12b 100644 --- a/boltzgen/lbm/model/__init__.py +++ b/boltzgen/lbm/model/__init__.py @@ -1,6 +1 @@ -from . import D2Q9 -from . import D3Q7 -from . import D3Q19 -from . import D3Q27 - -from . import characteristics +from .BGK import BGK diff --git a/boltzgen/lbm/model/characteristics.py b/boltzgen/lbm/model/characteristics.py deleted file mode 100644 index b68afeb..0000000 --- a/boltzgen/lbm/model/characteristics.py +++ /dev/null @@ -1,27 +0,0 @@ -from sympy import * - -# copy of `sympy.integrals.quadrature.gauss_hermite` sans evaluation -def gauss_hermite(n): - x = Dummy("x") - p = hermite_poly(n, x, polys=True) - p1 = hermite_poly(n-1, x, polys=True) - xi = [] - w = [] - for r in p.real_roots(): - xi.append(r) - w.append(((2**(n-1) * factorial(n) * sqrt(pi))/(n**2 * p1.subs(x, r)**2))) - return xi, w - -# determine weights of a d-dimensional LBM model on velocity set c -# (only works for velocity sets that result into NSE-recovering LB models when -# plugged into Gauss-Hermite quadrature without any additional arguments -# i.e. D2Q9 and D3Q27 but not D3Q19) -def weights(d, c): - _, omegas = gauss_hermite(3) - return list(map(lambda c_i: Mul(*[ omegas[1+c_i[iDim]] for iDim in range(0,d) ]) / pi**(d/2), c)) - -# determine lattice speed of sound using directions and their weights -def c_s(d, c, w): - speeds = set([ sqrt(sum([ w[i] * c_i[j]**2 for i, c_i in enumerate(c) ])) for j in range(0,d) ]) - assert len(speeds) == 1 # verify isotropy - return speeds.pop() diff --git a/boltzgen/utility/__init__.py b/boltzgen/utility/__init__.py index 8927868..fa9c760 100644 --- a/boltzgen/utility/__init__.py +++ b/boltzgen/utility/__init__.py @@ -1,2 +1,7 @@ from . import optimizations from . import ndindex + +from sympy.codegen.ast import Assignment + +def assign(names, definitions): + return list(map(lambda x: Assignment(*x), zip(names, definitions))) -- cgit v1.2.3