pystencils issueshttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues2021-12-03T16:56:29+01:00https://i10git.cs.fau.de/pycodegen/pystencils/-/issues/13show how to create own opts with sympy.codegen.rewriting.optimize in doc2021-12-03T16:56:29+01:00Stephan Seitzshow how to create own opts with sympy.codegen.rewriting.optimize in docIn `sympy.codegen.rewriting` there is a function `optimize`:
```python
def optimize(expr, optimizations):
""" Apply optimizations to an expression.
Parameters
==========
expr : expression
optimizations : iterable of ``Optimization`` instances
The optimizations will be sorted with respect to ``priority`` (highest first).
Examples
========
>>> from sympy import log, Symbol
>>> from sympy.codegen.rewriting import optims_c99, optimize
>>> x = Symbol('x')
>>> optimize(log(x+3)/log(2) + log(x**2 + 1), optims_c99)
log1p(x**2) + log2(x + 3)
"""
for optim in sorted(optimizations, key=lambda opt: opt.priority, reverse=True):
new_expr = optim(expr)
if optim.cost_function is None:
expr = new_expr
else:
before, after = map(lambda x: optim.cost_function(x), (expr, new_expr))
if before > after:
expr = new_expr
return expr
```
We should use it in `create_kernel` to profit from Sympy's optimizations (currently very few, but some are duplicates from optimizations in pystencils) and to make it easy for users to incorporate their own optimizations into the expressions. `create_kernel` could accept an `Iterable` of `Optimization`s with a default collection of optimizations that pystencils normally uses.
As you can see, it's really easy to implement own optimizations:
```python
def create_expand_pow_optimization(limit):
""" Creates an instance of :class:`ReplaceOptim` for expanding ``Pow``.
The requirements for expansions are that the base needs to be a symbol
and the exponent needs to be an Integer (and be less than or equal to
``limit``).
Parameters
==========
limit : int
The highest power which is expanded into multiplication.
Examples
========
>>> from sympy import Symbol, sin
>>> from sympy.codegen.rewriting import create_expand_pow_optimization
>>> x = Symbol('x')
>>> expand_opt = create_expand_pow_optimization(3)
>>> expand_opt(x**5 + x**3)
x**5 + x*x*x
>>> expand_opt(x**5 + x**3 + sin(x)**3)
x**5 + sin(x)**3 + x*x*x
"""
return ReplaceOptim(
lambda e: e.is_Pow and e.base.is_symbol and e.exp.is_Integer and abs(e.exp) <= limit,
lambda p: (
UnevaluatedExpr(Mul(*([p.base]*+p.exp), evaluate=False)) if p.exp > 0 else
1/UnevaluatedExpr(Mul(*([p.base]*-p.exp), evaluate=False))
))
```In `sympy.codegen.rewriting` there is a function `optimize`:
```python
def optimize(expr, optimizations):
""" Apply optimizations to an expression.
Parameters
==========
expr : expression
optimizations : iterable of ``Optimization`` instances
The optimizations will be sorted with respect to ``priority`` (highest first).
Examples
========
>>> from sympy import log, Symbol
>>> from sympy.codegen.rewriting import optims_c99, optimize
>>> x = Symbol('x')
>>> optimize(log(x+3)/log(2) + log(x**2 + 1), optims_c99)
log1p(x**2) + log2(x + 3)
"""
for optim in sorted(optimizations, key=lambda opt: opt.priority, reverse=True):
new_expr = optim(expr)
if optim.cost_function is None:
expr = new_expr
else:
before, after = map(lambda x: optim.cost_function(x), (expr, new_expr))
if before > after:
expr = new_expr
return expr
```
We should use it in `create_kernel` to profit from Sympy's optimizations (currently very few, but some are duplicates from optimizations in pystencils) and to make it easy for users to incorporate their own optimizations into the expressions. `create_kernel` could accept an `Iterable` of `Optimization`s with a default collection of optimizations that pystencils normally uses.
As you can see, it's really easy to implement own optimizations:
```python
def create_expand_pow_optimization(limit):
""" Creates an instance of :class:`ReplaceOptim` for expanding ``Pow``.
The requirements for expansions are that the base needs to be a symbol
and the exponent needs to be an Integer (and be less than or equal to
``limit``).
Parameters
==========
limit : int
The highest power which is expanded into multiplication.
Examples
========
>>> from sympy import Symbol, sin
>>> from sympy.codegen.rewriting import create_expand_pow_optimization
>>> x = Symbol('x')
>>> expand_opt = create_expand_pow_optimization(3)
>>> expand_opt(x**5 + x**3)
x**5 + x*x*x
>>> expand_opt(x**5 + x**3 + sin(x)**3)
x**5 + sin(x)**3 + x*x*x
"""
return ReplaceOptim(
lambda e: e.is_Pow and e.base.is_symbol and e.exp.is_Integer and abs(e.exp) <= limit,
lambda p: (
UnevaluatedExpr(Mul(*([p.base]*+p.exp), evaluate=False)) if p.exp > 0 else
1/UnevaluatedExpr(Mul(*([p.base]*-p.exp), evaluate=False))
))
```Release 1.0Stephan SeitzStephan Seitzhttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/20Revamp the type system2021-11-22T15:32:24+01:00Michael Kuronmkuron@icp.uni-stuttgart.deRevamp the type system@bauer had planned to do that before he left, but ran out of time. He has some pretty clear ideas of how the type system should be. It should probably be done relatively quickly as it would make quite a few things a lot easier, e.g. the SIMD stuff. It is also necessary before things like #12 can be done properly.
The following test (for test_vectorization.py) does not work currently:
```python
def test_vectorized_loop_counter():
arr = np.zeros((4, 4))
@ps.kernel
def kernel_equal(s):
f = ps.fields(f=arr)
f[0, 0] @= ps.astnodes.LoopOverCoordinate.get_loop_counter_symbol(1)
ast = ps.create_kernel(kernel_equal).compile()(f=arr)
arr_ref = arr.copy()
arr = np.zeros((4, 4))
vectorize(ast, instruction_set=instruction_set)
ps.show_code(ast)
ast.compile()(f=arr)
assert np.allclose(arr_ref, arr)
```
Another problem can be found in the [attached jupyter notebook](/uploads/04f89bc361e00c144f14e5d657d78c4e/Untitled1.ipynb), where in the absence of `type_all_numbers` floating-point comparisons are used for integers in a conditional.
Finally, automatic conversion of arguments to the type expected by the signature of the SIMD intrinsics is not supported, which made https://i10git.cs.fau.de/walberla/walberla/-/merge_requests/414/diffs necessary.@bauer had planned to do that before he left, but ran out of time. He has some pretty clear ideas of how the type system should be. It should probably be done relatively quickly as it would make quite a few things a lot easier, e.g. the SIMD stuff. It is also necessary before things like #12 can be done properly.
The following test (for test_vectorization.py) does not work currently:
```python
def test_vectorized_loop_counter():
arr = np.zeros((4, 4))
@ps.kernel
def kernel_equal(s):
f = ps.fields(f=arr)
f[0, 0] @= ps.astnodes.LoopOverCoordinate.get_loop_counter_symbol(1)
ast = ps.create_kernel(kernel_equal).compile()(f=arr)
arr_ref = arr.copy()
arr = np.zeros((4, 4))
vectorize(ast, instruction_set=instruction_set)
ps.show_code(ast)
ast.compile()(f=arr)
assert np.allclose(arr_ref, arr)
```
Another problem can be found in the [attached jupyter notebook](/uploads/04f89bc361e00c144f14e5d657d78c4e/Untitled1.ipynb), where in the absence of `type_all_numbers` floating-point comparisons are used for integers in a conditional.
Finally, automatic conversion of arguments to the type expected by the signature of the SIMD intrinsics is not supported, which made https://i10git.cs.fau.de/walberla/walberla/-/merge_requests/414/diffs necessary.Release 1.0Markus HolzerMarkus Holzerhttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/29Aligned GPU Array2020-12-10T10:40:02+01:00Markus HolzerAligned GPU ArrayAt the moment only aligned CPU arrays exist. For GPU arrays it should also be possible to pass the alignment and create an array of which.
Possible implementation:
```
base_array = gpuarray.GPUArray((array.size + 8), dtype=array.dtype)
gpu_array = gpuarray.GPUArray(
shape=tuple(array.shape),
dtype=array.dtype,
base=base_array,
gpudata=int(base_array.gpudata) + 24,
strides=tuple(gpuarray._compact_strides(array)),)
```At the moment only aligned CPU arrays exist. For GPU arrays it should also be possible to pass the alignment and create an array of which.
Possible implementation:
```
base_array = gpuarray.GPUArray((array.size + 8), dtype=array.dtype)
gpu_array = gpuarray.GPUArray(
shape=tuple(array.shape),
dtype=array.dtype,
base=base_array,
gpudata=int(base_array.gpudata) + 24,
strides=tuple(gpuarray._compact_strides(array)),)
```https://i10git.cs.fau.de/pycodegen/pystencils/-/issues/33Optimisation of numerical conditioning2021-10-14T22:56:04+02:00Markus HolzerOptimisation of numerical conditioningNumerical conditioning is very important, especially when using single precision. There are several cases that have been shown in the literature where this plays an important role. First storing only the deviation from the equilibrium of the particle distributions of a lattice Boltzmann simulation instead of storing the equilibrium. Second in the advanced lattice Boltzmann methods like the cumulant LBM where it was shown that numerical conditioning is an extremely important task to get a stable simulation.
This task is very tedious and complicated. However, a starting point could be to look at [pyinterval](https://pypi.org/project/pyinterval/) and investigate the error which comes from certain operations. This could lead to knowledge on minimizing errors arising from numerical methods.Numerical conditioning is very important, especially when using single precision. There are several cases that have been shown in the literature where this plays an important role. First storing only the deviation from the equilibrium of the particle distributions of a lattice Boltzmann simulation instead of storing the equilibrium. Second in the advanced lattice Boltzmann methods like the cumulant LBM where it was shown that numerical conditioning is an extremely important task to get a stable simulation.
This task is very tedious and complicated. However, a starting point could be to look at [pyinterval](https://pypi.org/project/pyinterval/) and investigate the error which comes from certain operations. This could lead to knowledge on minimizing errors arising from numerical methods.https://i10git.cs.fau.de/pycodegen/pystencils/-/issues/36Test case for FVM source term and fluctuations2021-11-19T19:11:14+01:00Michael Kuronmkuron@icp.uni-stuttgart.deTest case for FVM source term and fluctuationsPlease create a pystencils test case from your reactive electrokinetics code. We currently don't have test coverage for the discretization of the source term. While you're add it, include the fluctuating EK test as well.Please create a pystencils test case from your reactive electrokinetics code. We currently don't have test coverage for the discretization of the source term. While you're add it, include the fluctuating EK test as well.IngoIngohttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/37C or C++2021-11-19T15:53:42+01:00Jan HönigC or C++I have stumbled into the issue of C vs C++ with `pystencils`. 99% we generate C code. We even would like to use the `restrict` keyword, which is only available in C, not in C++ (although many compilers do support it).
I am not sure if we misuse C/C++ at this point. Generally C only would be beneficial, so we can include our kernels in pretty much any language, since any relevant language has a c-ffi.
What would we miss, if we would restrict ourselves on C only?
So far I have found some C++ allocations (#24), and `std::cout` usage.I have stumbled into the issue of C vs C++ with `pystencils`. 99% we generate C code. We even would like to use the `restrict` keyword, which is only available in C, not in C++ (although many compilers do support it).
I am not sure if we misuse C/C++ at this point. Generally C only would be beneficial, so we can include our kernels in pretty much any language, since any relevant language has a c-ffi.
What would we miss, if we would restrict ourselves on C only?
So far I have found some C++ allocations (#24), and `std::cout` usage.Jan HönigJan Hönighttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/42Logging2021-11-19T15:53:54+01:00Markus HolzerLoggingPython provides very good logging functionality which should be used in pystencils. For example problems in the type-system could be logged on a debug level. Oftentimes it is necessary to add print statements in a debug session to see which expressions are cast to which type etc. This could always be logged on a debug level so we can just always see what happens in the back.
https://docs.python.org/3/howto/logging.htmlPython provides very good logging functionality which should be used in pystencils. For example problems in the type-system could be logged on a debug level. Oftentimes it is necessary to add print statements in a debug session to see which expressions are cast to which type etc. This could always be logged on a debug level so we can just always see what happens in the back.
https://docs.python.org/3/howto/logging.htmlMarkus HolzerMarkus Holzerhttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/45Modularise FD, Runhelper and others2021-11-19T15:54:46+01:00Markus HolzerModularise FD, Runhelper and othersFd and Runhelper should be modules that build upon pystencils. They should not be integrated directly in pystencilsFd and Runhelper should be modules that build upon pystencils. They should not be integrated directly in pystencilsMarkus HolzerMarkus Holzerhttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/46Vectorization revamp2021-11-19T15:54:42+01:00Jan HönigVectorization revamphttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/47Clean extended SymPy functions2021-11-19T16:46:28+01:00Markus HolzerClean extended SymPy functionsThere are some functions which extend a SymPy function. These should all be implemented in a single module/fileThere are some functions which extend a SymPy function. These should all be implemented in a single module/fileRelease 1.0Markus HolzerMarkus Holzerhttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/48Python 3.102021-11-19T16:49:22+01:00Jan HönigPython 3.10We switch to `python3.10` when `conda` and the latest Ubuntu LTS version support that.
Apply pattern matching at important functions.We switch to `python3.10` when `conda` and the latest Ubuntu LTS version support that.
Apply pattern matching at important functions.Markus HolzerMarkus Holzerhttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/50CustomCodeNode should be treated like SympyAssignment in transformations (cur...2021-11-29T19:31:41+01:00itischlerCustomCodeNode should be treated like SympyAssignment in transformations (currently RNG numbers are calculated too often)I don't know if it is a bug or if I am using it wrong, but when using mutliple rng symbols via `next` on e.g. fluxes in different directions, the generated code will calculate each random number when only one is used. so basically for a D3Q19 stencil Philox is called 5 times more then is needed, which slows down the program drastically.
```
import pystencils as ps
import sympy as sp
from pystencils.rng import random_symbol
# Setup
stencil = 19
L=(10,10,10)
dh = ps.create_data_handling(domain_size=L, periodicity=True)
c_field = dh.add_array('c', values_per_cell=1)
j_field = dh.add_array('j', values_per_cell=stencil // 2,
field_type=ps.FieldType.STAGGERED_FLUX)
def grad(f):
return sp.Matrix([ps.fd.diff(f, i) for i in range(dh.dim)])
# genergate flux eq. and set them to 0
fvm_eq = ps.fd.FVM1stOrder(c_field, flux=grad(c_field))
flux=[]
for diff in fvm_eq.discrete_flux(j_field):
flux.append(ps.Assignment(diff.lhs, 0))
flux = ps.AssignmentCollection(flux)
# define rng
rng_symbol_gen = random_symbol(flux.subexpressions, dim=dh.dim)
for i in range(len(flux.main_assignments)):
fluct = next(rng_symbol_gen)
flux.main_assignments[i] = ps.Assignment(flux.main_assignments[i].lhs, fluct)
print(flux)
flux_kernel = ps.create_staggered_kernel(flux).compile()
ps.show_code(flux_kernel)
```I don't know if it is a bug or if I am using it wrong, but when using mutliple rng symbols via `next` on e.g. fluxes in different directions, the generated code will calculate each random number when only one is used. so basically for a D3Q19 stencil Philox is called 5 times more then is needed, which slows down the program drastically.
```
import pystencils as ps
import sympy as sp
from pystencils.rng import random_symbol
# Setup
stencil = 19
L=(10,10,10)
dh = ps.create_data_handling(domain_size=L, periodicity=True)
c_field = dh.add_array('c', values_per_cell=1)
j_field = dh.add_array('j', values_per_cell=stencil // 2,
field_type=ps.FieldType.STAGGERED_FLUX)
def grad(f):
return sp.Matrix([ps.fd.diff(f, i) for i in range(dh.dim)])
# genergate flux eq. and set them to 0
fvm_eq = ps.fd.FVM1stOrder(c_field, flux=grad(c_field))
flux=[]
for diff in fvm_eq.discrete_flux(j_field):
flux.append(ps.Assignment(diff.lhs, 0))
flux = ps.AssignmentCollection(flux)
# define rng
rng_symbol_gen = random_symbol(flux.subexpressions, dim=dh.dim)
for i in range(len(flux.main_assignments)):
fluct = next(rng_symbol_gen)
flux.main_assignments[i] = ps.Assignment(flux.main_assignments[i].lhs, fluct)
print(flux)
flux_kernel = ps.create_staggered_kernel(flux).compile()
ps.show_code(flux_kernel)
```Markus HolzerMarkus Holzerhttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/51Installing Pystencils2021-11-29T10:21:20+01:00Markus HolzerInstalling PystencilsFrom time to time problems arise when installing `pystencils`. The latest problem for example was that `pip` was not installed in a new `conda` environment and then for compiling the boundary C-File was linked with the wrong python version. Problems like these should be collected in an FAQ file or something similar.
Just by comparing with `NumPy` there is quite a long install prerequisite file although the installation of `NumPy` is really straightforward. Something like this I already imagine for `pystencils`:
https://numpy.org/install/From time to time problems arise when installing `pystencils`. The latest problem for example was that `pip` was not installed in a new `conda` environment and then for compiling the boundary C-File was linked with the wrong python version. Problems like these should be collected in an FAQ file or something similar.
Just by comparing with `NumPy` there is quite a long install prerequisite file although the installation of `NumPy` is really straightforward. Something like this I already imagine for `pystencils`:
https://numpy.org/install/Markus HolzerMarkus Holzerhttps://i10git.cs.fau.de/pycodegen/pystencils/-/issues/52Add possiblity to create ".pvd"-files in vtk-writer2021-12-08T11:35:18+01:00Christoph SchwarzmeierAdd possiblity to create ".pvd"-files in vtk-writerCurrently, the files that are written in pystencil's vtk-writer are stored as _vtkImageData_ (`.vti`) only. It would be great if the vtk-writer was able to create a file in _ParaView Data format_ (`.pvd`) which contains the path to each `.vti`-file and the (LBM) time step at which the file was written.
In ParaView, "time" would then actually represent the (LBM) time step instead of being equal to the index, i.e., number of the loaded ".vti"-file.Currently, the files that are written in pystencil's vtk-writer are stored as _vtkImageData_ (`.vti`) only. It would be great if the vtk-writer was able to create a file in _ParaView Data format_ (`.pvd`) which contains the path to each `.vti`-file and the (LBM) time step at which the file was written.
In ParaView, "time" would then actually represent the (LBM) time step instead of being equal to the index, i.e., number of the loaded ".vti"-file.