diff options
18 files changed, 106 insertions, 95 deletions
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/model/D2Q9.py b/boltzgen/lbm/lattice/D2Q9.py index e3ac9de..e3ac9de 100644 --- a/boltzgen/lbm/model/D2Q9.py +++ b/boltzgen/lbm/lattice/D2Q9.py diff --git a/boltzgen/lbm/model/D3Q19.py b/boltzgen/lbm/lattice/D3Q19.py index e9e6eec..e9e6eec 100644 --- a/boltzgen/lbm/model/D3Q19.py +++ b/boltzgen/lbm/lattice/D3Q19.py diff --git a/boltzgen/lbm/model/D3Q27.py b/boltzgen/lbm/lattice/D3Q27.py index 6fb1f80..6fb1f80 100644 --- a/boltzgen/lbm/model/D3Q27.py +++ b/boltzgen/lbm/lattice/D3Q27.py diff --git a/boltzgen/lbm/model/D3Q7.py b/boltzgen/lbm/lattice/D3Q7.py index 04e16a3..04e16a3 100644 --- a/boltzgen/lbm/model/D3Q7.py +++ b/boltzgen/lbm/lattice/D3Q7.py 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/model/characteristics.py b/boltzgen/lbm/lattice/characteristics.py index b68afeb..b68afeb 100644 --- a/boltzgen/lbm/model/characteristics.py +++ b/boltzgen/lbm/lattice/characteristics.py 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/__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/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))) |