diff options
| -rw-r--r-- | standalone_cpp_codegen.py | 32 | ||||
| -rw-r--r-- | template/standalone.mako | 109 |
2 files changed, 141 insertions, 0 deletions
diff --git a/standalone_cpp_codegen.py b/standalone_cpp_codegen.py new file mode 100644 index 0000000..b38898a --- /dev/null +++ b/standalone_cpp_codegen.py @@ -0,0 +1,32 @@ +import sympy +from mako.template import Template +from pathlib import Path + +from simulation import Geometry +from symbolic.generator import LBM + +import symbolic.D2Q9 as D2Q9 + +lbm = LBM(D2Q9) + +moments = lbm.moments(optimize = False) +collide = lbm.bgk(f_eq = lbm.equilibrium(), tau = 0.6) + +geometry = Geometry(512, 512) + +program_src = Template(filename = str(Path(__file__).parent/'template/standalone.mako')).render( + descriptor = lbm.descriptor, + geometry = geometry, + + steps = 100, + + moments_subexpr = moments[0], + moments_assignment = moments[1], + collide_subexpr = collide[0], + collide_assignment = collide[1], + + float_type = 'double', + ccode = sympy.ccode +) + +print(program_src) diff --git a/template/standalone.mako b/template/standalone.mako new file mode 100644 index 0000000..6750de9 --- /dev/null +++ b/template/standalone.mako @@ -0,0 +1,109 @@ +#include <array> +#include <cstdint> +#include <memory> + +<% +def pop_offset(i): + return i * geometry.volume +%> + +void equilibrilize(${float_type}* f_next, + ${float_type}* f_prev, + const std::size_t gid) +{ + ${float_type}* preshifted_f_next = f_next + gid; + ${float_type}* preshifted_f_prev = f_prev + gid; + +% for i, w_i in enumerate(descriptor.w): + preshifted_f_next[${pop_offset(i)}] = ${w_i}.f; + preshifted_f_prev[${pop_offset(i)}] = ${w_i}.f; +% endfor +} + +<% +def neighbor_offset(c_i): + return { + 2: lambda: c_i[1]*geometry.size_x + c_i[0], + 3: lambda: c_i[2]*geometry.size_x*geometry.size_y + c_i[1]*geometry.size_x + c_i[0] + }.get(descriptor.d)() + +%> + +void collide_and_stream( ${float_type}* f_next, + const ${float_type}* f_prev, + const int* material, + const std::size_t gid) +{ + const int m = material[gid]; + + if ( m == 0 ) { + return; + } + + ${float_type}* preshifted_f_next = f_next + gid; + const ${float_type}* preshifted_f_prev = f_prev + gid; + +% for i, c_i in enumerate(descriptor.c): + const ${float_type} f_curr_${i} = preshifted_f_prev[${pop_offset(i) + neighbor_offset(-c_i)}]; +% endfor + +% for i, expr in enumerate(moments_subexpr): + const ${float_type} ${expr[0]} = ${ccode(expr[1])}; +% endfor + +% for i, expr in enumerate(moments_assignment): + ${float_type} ${ccode(expr)} +% endfor + +% for i, expr in enumerate(collide_subexpr): + const ${float_type} ${expr[0]} = ${ccode(expr[1])}; +% endfor + +% for i, expr in enumerate(collide_assignment): + const ${float_type} ${ccode(expr)} +% endfor + +% for i in range(0,descriptor.q): + preshifted_f_next[${pop_offset(i)}] = f_next_${i}; +% endfor +} + +int main() +{ + auto f_a = std::make_unique<${float_type}[]>(${geometry.volume*descriptor.q}); + auto f_b = std::make_unique<${float_type}[]>(${geometry.volume*descriptor.q}); + auto material = std::make_unique<int[]>(${geometry.volume}); + + ${float_type}* f_prev = f_a.get(); + ${float_type}* f_next = f_b.get(); + + for (int iX = 0; iX < ${geometry.size_x}; ++iX) { + for (int iY = 0; iY < ${geometry.size_x}; ++iY) { + if (iX == 0 || iY == 0 || iX == ${geometry.size_x-1} || iY == ${geometry.size_y-1}) { + material[iY*${geometry.size_x} + iX] = 0; + } else { + material[iY*${geometry.size_x} + iX] = 1; + } + } + } + + for (std::size_t iCell = 0; iCell < ${geometry.volume}; ++iCell) { + equilibrilize(f_prev, f_next, iCell); + } + + for (std::size_t iStep = 0; iStep < ${steps}; ++iStep) { + if (iStep % 2 == 0) { + f_next = f_a.get(); + f_prev = f_b.get(); + } else { + f_next = f_b.get(); + f_prev = f_a.get(); + } + + for (std::size_t iCell = 0; iCell < ${geometry.volume}; ++iCell) { + collide_and_stream(f_next, f_prev, material.get(), iCell); + } + } + + return 0; +} |