lettucecfd / lettuce Goto Github PK

In the .gitignore "docs/make.bat" and "docs/Makefile" are specifically ignored.
Nevertheless, both files exist in the repo.

The .gitignore should either not ignore these files or the files should be removed from the repo.

The way we initialize the grid may lead to confusion and inconsistencies. Since we rely on the standard indexing 'xy' of numpy.meshgrid, the x-axis lies along index 1, while the y-axis lies along index 0. (The z-axis lies always along index 2)

For example in the shear-layer test case, the final velocity field had to be set to:

# switching to ij/matrix-indexing -> 1st entry: i = -y (i/line index is going down instead of up like y), 2nd entry: x = j (column index)

u = np.array([-uy, ux])

We could circumvent this situation by setting numpy.meshgrid(indexing='ij'). As a result, the x-axis would be index 0 and so on. I think this would be more intuitive.

Although, the indexing could be set individually for each test flow, we should change it consistently in Lettuce.
What do you think? Any concerns?

Hi
When I ran "python setup.py test" on Ubuntu. The prcoess was interuptted by this error
[subprocess.CalledProcessError: Command '['.\CONDA\envs\lettuce\python.exe', 'setup.py', 'install']' returned non-zero exit status 1.]
How can I deal with it?
Thank you for your help.

Rho is required for the calculation of equilibrium and u (and u_eq).
Therefor it should be pre calculated for both.

Atm:

class lattices
    def u(self, f):
        return self.j(f) / self.rho(f)

class BGKCollision:
    def __call__(self, f):
        rho = self.lattice.rho(f)
        u_eq = 0 if self.force is None else self.force.u_eq(f)
        u = self.lattice.u(f) + u_eq
        feq: torch.Tensor = self.lattice.equilibrium(rho, u)
        Si = 0 if self.force is None else self.force.source_term(u)
        return f - 1.0 / self.tau * (f - feq) + Si

My suggestion:

class lattices
    def u(self, f, rho_):
        return self.j(f) / rho_

class BGKCollision:
    def __call__(self, f):
        rho_ = self.lattice.rho(f)
        u_eq = 0 if self.force is None else self.force.u_eq(f)
        u = self.lattice.u(f, rho_) + u_eq
        feq: torch.Tensor = self.lattice.equilibrium(rho_, u)
        Si = 0 if self.force is None else self.force.source_term(u)
        return f - 1.0 / self.tau * (f - feq) + Si

Hi,
I executed the A-first-example.ipynp, varying num_steps as 10000 and 20000. Surprisingly, the output remained consistent across both values. Strikingly, even in the absence of initial simulation steps, the figure replicated the outcomes from num_steps=10000 and 20000. Curiously, despite reaching an energy of zero at 20000 steps, the output figure didn't reflect this and remained vibrant instead of being a plain dark blue plane, as expected.

Furthermore, I experimented with different reynolds_number values 4 and 4000000. However, the resulting figures were consistently identical, showing no variation based on this parameter.

Could you kindly provide an explanation for these observations?

Can we do Multi-GPU over MPI?

Here is a summary of distributed computing capabilities
https://pytorch.org/docs/stable/distributed.html

AFAIK there are no "distributed tensors" or anything like that. But there are things like sending and receiving tensors via MPI. This means that we have to set up the ghost layers in a distributed simulation by hand (only in the streaming step and grid generation).

This requires a version of PyTorch that is built with MPI (I don't think they have these on conda, so we need to build it by hand).

Here is a list of things that would be required on the way to a usable multinode implementation:

• A grid class. Instead of being defined as a np.array inside the flow, the grid needs to be an instance of a class RegularGrid.
• When an MPI process accesses the grid, it should only see the rectangular part of the grid that is associated with its domain plus the ghost layers. In this way, all the distribution functions etc. have the dimension of the local domain.
• The grid class also stores information about the ghost layers and the processes that are associated with the neighbor domain.
• The streaming operator has to do the roll operation plus the communication of information in the ghost layer via torch.distributed.send() and torch.distributed.recv().
• When computing global quantities (energy, ...) we need to do a torch.distributed.allreduce().
• Distributed I/O of vtk files is possible (see NATriuM) , but I don't know if it is supported by pyevtk, yet. If not, it should be quite easy to add.

Here is an example of distributed computing in pytorch https://www.glue.umd.edu/hpcc/help/software/pytorch.html

At this point, I am not sure how efficient and feasible such an MPI implementation would be. We should try this outside lettuce first for a mini example before proceeding.

First, thanks for this awesome Pytorch-based LBM library.

I am new to LBM and have a question about grid_fine_to_coarse function. I do not get the idea why we set the relaxation parameter as 2 * tau_coarse / tau_fine as follow:
f_coarse = f_eq + 2 * tau_coarse / tau_fine * f_neq

Thanks

TODO: Add examples as jupyter notebooks and convert the existing one from .py to .ipynb as well.
See #22

My simulation in porous_medium.ipynb returned:

/home/philipp/lettuce/lettuce/simulation.py:122: UserWarning: where received a uint8 condition tensor. This behavior is deprecated and will be removed in a future version of PyTorch. Use a boolean condition instead. (Triggered internally at /opt/conda/conda-bld/pytorch_1678411187366/work/aten/src/ATen/native/TensorCompare.cpp:493.)
  self.f = torch.where(self.collision.no_collision_mask, self.f, self.collision(self.f))

This originates here:

I don't know if you're already on it, just wanted to let you know.

Atm not all files are pep conform.
It would be best practise if we would reformat the code.

I can't seem to find reference to any neural network model being used as the relaxation parameters (unless the output was hardcoded, but even then, it seems generalised for different flow types as the stencils are interchangeable). An indication in the code of where this happens would be good. Additionally, I also can't seem to find GPU acceleration through torch.autograd, was this feature cut?

Hi,

I want to initialize proper Dirichlet Boundary Conditions, e.g. of the Form

\u_D = \u_{avg} * \frac{6y(H-y)} {H^2} * w(t),

with w(t) given as
w(t) = 0.5 - 0.5cos(2pi*t) if t < 0.5 and 1 if t > 0.5.

In the Krüger et. al. book, there is a short chapter on proper initialization. They mention the method by Mei et. al. However, this approach is only for one single time step (most likely t=0). How to handle the above case. Use the iterative approach of Mei, repeat it for subsequent time steps and set the initial condition for every time step separately? Seems a little bit odd to me.

Furthermore, there is a method in the Simulation class, that initializes f such that it matches the given moments. It looks pretty similar to the approach of Mei et. al., with the only difference that the termination criterion is pressure-dependent, Mei et. al. propose a density-dependent criterion.

Any suggestions on that?

Best,
Robert

remove deprecated reporters (MaxUReporter,...) that have been replaced by ObservableReporter

• add a couette flow example
• add an analytic_solution
• update boundary condition for top and bottom plates (currently too low at < 1e-6)

The main parameter beside resolution is the steps parameter.
This indicates (for me at least) that as you make more steps the benchmark is more meaningful.

That is somehow true. But this also highly depends on the resolution.
E.g.:

1. You benchmark with a resolution of 4048 and 1000 steps. The result of the benchmark is 400 Mlups.

Basically you told the benchmark to calculate 4,048 Nodes * 4,048 Nodes * 1,000 Steps = 16,386,304,000 Nodes
The benchmark took 16,386,304,000 Nodes / 400,000,000 lups = ~41 Seconds

2. You benchmark with a resolution of 16 and 1000 steps. The result of the benchmark is 350 Mlups.

Basically you told the benchmark to calculate 16 Nodes * 16 Nodes * 1,000 Steps = 256,000 Nodes
The benchmark took 256,000 Nodes / 350,000,000 lups = ~0.0007314 Seconds

There one big issue in comparing these two benchmarks:

• The time the second benchmark took is so short that you would say it is not accurate at all. There is just too much noise in the performance of a GPU for such a short benchmark to be accurate
• On the other hand the first benchmark seems much more accurate.

What I described is not a Bug but it leads the user to misuse the api.
I would consider to either change the benchmark parameter to Seconds or to TotalNodes (to process).

I'm facing a challenge: I'd like to simulate the TaylorGreenVortex2D with an added force. I'm uncertain about how to modify the code to achieve this.
Any assistance would be greatly appreciated.

where

All the best,

TGV2D with initialize_pressure seems to produce larger errors than without.

ctx = {"device": torch.device("cuda:0"), "dtype": torch.float64, "use_native": False}

lattice = lt.Lattice(lt.D2Q9, **ctx)
flow = lt.TaylorGreenVortex2D(
    resolution=32,
    reynolds_number=1000,
    mach_number=0.05,
    lattice=lattice
)

In file lettuce/boundary.py:BounceBackBoundary.__init__

"I actually don't think these lines will ever throw a LettuceException. Rather, they will raise AssertionError which does not seem to be the desired behavior."
Originally posted by @Olllom in #80 (comment)

Hi, i am a new lettuce users. When i run the example of porous_medium, there is an error:

 class PorousMedium2D(lt.Obstacle2D):

TypeError: function() argument 1 must be code, not str

I check the code, and notice that the lt.Obstacle2D is a function defined in obstacle.py:

def Obstacle2D(resolution_x, resolution_y, reynolds_number, mach_number, lattice, char_length_lu):
warnings.warn("Obstacle2D is deprecated. Use Obstacle instead", DeprecationWarning)
shape = (resolution_x, resolution_y)
domain_length_x = resolution_x / char_length_lu
return Obstacle(shape, reynolds_number, mach_number, lattice, domain_length_x=domain_length_x)

Therefore, in my python, class PorousMedium2D can not inherit from the lt.Obstacle2D as the type of lt.Obstacle2D is not a class.

I wanna know what i can do to solve the problem.

Thanks in advanced for your answer!

Hi,

the current implementation of the HDF5 writer saves the flow object as metadata. However, this limits the grid size since the metadata should not exceed 64kb (if I remember correctly). E.g. for large resolutions on the obstacle flow, the mask is limiting the usage. (if my assumption is correct). A quick fix would be, to save the mask into a separate dataset and set it to None before pickling the flow. Therefore, the assertion test has to be removed from the mask setter. :)

Best
Robert

Dear lettuce friends,

I have been using your code for a few weeks. I have been very successful with the 2D single phase model.

Now, I would like to extend the code to perform flow sims in 3D. In my first attempt I notice that the iteration speed did not scale very well, I wanted to see if you guys have some ideas on how to implement the model better.

o. In the example below, I have a lot of solid. Is it possible to mask the operations there to speed-up the iterations? Or can you think about something else? The current problem 256^3 takes about 1 hr on a A100 (for reference, a parallel code in ~4 cores would finish in 1 minute)

o. I'm also getting weird pressure values. Can you spot anything dumb in my implementation?

https://github.com/je-santos/lettuce/blob/bug/examples/3d_frac.ipynb

Many thanks for all your help, your code is great.

-Javier

Hi,

I am quite unsure about the implementation of the obstacle class in 2d and 3d. The grid is initialized with stop_x and stop_y by dividing the resolution with the characteristic length which is given in lattice units. However, this changes the geometry I am using when I increase the resolution. It may be, that I got something wrong since I am from a FEM background, but this seems odd to me. I would be happy about an explanation. Shall I subclass the obstacle class if I want to change that (e.g. for using a fixed geometry)?

def grid(self):
    stop_x = self.resolution_x / self.units.characteristic_length_lu
    stop_y = self.resolution_y / self.units.characteristic_length_lu

    x = np.linspace(0, stop_x, num=self.resolution_x, endpoint=False)
    y = np.linspace(0, stop_y, num=self.resolution_y, endpoint=False)

    return np.meshgrid(x, y, indexing='ij')

Best.
Robert

• Fix comments and typos in example notebooks.
• Convert existing .py example into notebook (or introduce the vtk writer in a new example notebook).

@McBs

The computational bottleneck in LBM simulations is memory access. The current code framework handles that very naively. To boost the performance, the collision and streaming routines can be combined, like for example in the esoteric twist method.
This would require to step beyond the standard routines offered by PyTorch and write custom C++/CUDA PyTorch extensions:
https://pytorch.org/tutorials/advanced/cpp_extension.html

I tested the porous_medium.ipynb and it was pretty fast on the GPU. Can it be used to simulate a two-phase flow? How to set the parameters? Thank you!

The way I understand this is that the network minimizes the difference between fine and coarse grid results. So once the network is trained you can apply to a field computed on a corase grid to recover the accurate results. Is this right?

Dear,
I followed the installation steps on my Linux workstation.

When I run the command python setup.py install

This error showed up:
setuptools.extern.packaging.version.InvalidVersion: Invalid version: 'master'

Thanks in advance

As recent tests showed, it would be good to have a test case that runs the convergence tool to identify the root of other tests failing.

Flow5 is a great tool for foils. It is the successor of XFoil and Xflsr5.
https://flow5.tech/

It can use with opencascade cad with stL, iges and more.

Do you see here a fast solution for coupling your 2d-Tool to improve IBL-solver
to a better solutions?
https://flow5.tech/index.php/flow5/new-features/2d-ibl-solver/

Some examples don't work with the new setup. Common errors are

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[9], line 1
----> 1 u = flow.units.convert_velocity_to_pu(lattice.u(simulation.f)).numpy()
      2 u_norm = np.linalg.norm(u,axis=0)
      3 plt.imshow(u_norm)

TypeError: can't convert cuda:0 device type tensor to numpy. Use Tensor.cpu() to copy the tensor to host memory first.

and some issues caused by use_native=True

Now that we have PEP 8-compatible formatting (#80), we should try to keep it. It would be good to add a PEP8 hook to the PR workflow.

Dear lettucers,

Would it be feasible to run a single-phase flow through a porous medium using the current implementation?

Something like this:

https://github.com/je-santos/MPLBM-UT/tree/master/examples/singlephase_and_tortuosity_DRPcarbonate

Thanks a lot for this project, it's pretty sweet

Hey there, I was curious if there was documentation/recommendations on best practice for preparing an stl to work in lettuce for LBM aerodynamic testing. I'm curious if there's any recommended tools for creating a mesh from an stl and if there is any guidance on a workflow in lettuce.

thanks!

`import lettuce as lt

class PorousMedium2D(lt.Obstacle2D):
...
`
lt.Obstacle2D is recognized by Python as a function， PorousMedium2D class cannot be compiled.The error prompt is ‘TypeError: function() argument 'code' must be code, not str‘.
what can I do ?

As discussed before, we might want to do ourselves a pleasure and get rid of the LatticeAoS class (abstracting all the indexing has taken too much effort up to this point and we want to keep the code as simple as possible).

Before making the final decision, Dominik is going to run some benchmarks on the Volta nodes. If the difference in performance turns out to be marginal (or in favor of the structure-of-array Lattice), let's get rid of the thing asap.

TODO (@dominikwilde): Please post the benchmark results here.

Tolerances are too loose for CUDA and float32.

One solution would be that the dtype_device fixture also returns a "base tolerance" that is suitable for the given compute context.

如果出现类似“The detected CUDA version (12.3) mismatches the version that was used to compile
PyTorch (11.8). Please make sure to use the same CUDA versions.”，请找到检测版本的代码注释掉，就ok

According to the diagram in Krüger page 67 (Fig. 3.2) "An overview of one cycle of the LB algorithm" the proper order of algorithm steps is:

increment time > output > collision > propagation > boundaries

however in lettuce we implemented:

self.i += 1
self.f = self.streaming(self.f)
self.f = torch.where(self.no_collision_mask, self.f, self.collision(self.f))
for boundary in self._boundaries:
  self.f = boundary(self.f)
self._report()

so increment time > propagation > collision > boundaries > output

which doesn't agree, so I think it means that the result of boundaries is streamed away instead of replacing the streaming step for the respective fs (because of the no stream masks they wont be overwritten but will instead be on both n and n-1 nodes)#
also I think we are reporting a state where the time step isn't fully over

Please have a look and see if we should fix this @Olllom @McBs @dominikwilde

Hi,

I was trying to benchmark the lettuce single-phase porous media model with my existing codes, and I found something curious, which might be caused by some dumb mistake on my end.

In the simulation results in this notebook, in the last plot, you can see that the fluid velocity increases with distance from the inlet, theoretically which should be constant through the domain.

Let me know if you can spot the mistake.

Many thanks!!!!

examples
in obstacle.ipynb, How to add a fluid source of air pollen particles to add new fluids to the space?

I would like to understand how the main code is able to access the ml part, which AFAIU does the collision model per slide # 7 of this deck

Hi, here I will be listing things from my fork that I think might be beneficial for lettuce as a whole, as discussed previously.

(Still WIP on this list)

• Option to cancel simulation once any entry of f is NaN: https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/simulation.py#L91 with a delay so the output of the VTK-reporter that shows where NaN is as mask, has time to work: https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/reporters.py#L60

• If more inspiration or something for the parallelisation is needed: https://github.com/MartinKliemank/lettuce/blob/lettuceMPI/lettuce/distributed.py; https://github.com/MartinKliemank/lettuce/blob/lettuceMPI/lettuce/grid.py

• Option to have domain in CPU-RAM and collide two halfs separately in GPU (a bit faster than CPU only, larger domain than GPU only), idea by McBs: https://github.com/MartinKliemank/lettuce/blob/69df1fadc837a00286cc1a101c4bc53b21ea540f/lettuce/simulation.py#L261

Reporter / Observables

• VTK-Reporter expanded to be able to output bounce-back mask as VTK-file: https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/reporters.py#L82

• VTK-Reporter that only outputs pressure to save disk space (maybe integrate this with the normal VTK reporter and offer a choice of parameters to write?): https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/reporters.py#L98

• Drag Coefficient Observable (probably needs proper validation but was somewhat tested): https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/observables.py#L119

• Reporter for pressure at given coordinates: https://github.com/MartinKliemank/lettuce/blob/69df1fadc837a00286cc1a101c4bc53b21ea540f/lettuce/observables.py#L367, old version: https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/observables.py#L153

• Reporter for velocity at given coordintes: https://github.com/MartinKliemank/lettuce/blob/69df1fadc837a00286cc1a101c4bc53b21ea540f/lettuce/observables.py#L382

• Reporter for real time duration of a time step (maybe not the greatest idea): https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/observables.py#L139

• Reporter for current time step number (i): https://github.com/MartinKliemank/lettuce/blob/69df1fadc837a00286cc1a101c4bc53b21ea540f/lettuce/reporters.py#L309

Boundary conditions:

Sadly I am still not sure if the implementation of many these is proper, so they might need proper validation (there are even more boundaries there, but I'm not sure how far along they were :x). Most of these have atleast been running and produced probable results^^

• Synthetic Eddy Inlet (Buffa 2021): https://github.com/MartinKliemank/lettuce/blob/69df1fadc837a00286cc1a101c4bc53b21ea540f/lettuce/boundary.py#L283 works quite well I guess^^ needs a lot of inputs that will be in the upload to the Haus-Paper

• Equilibrium Boundary that extrapolates u from neighbour (pressure outlet): https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/boundary.py#L272

• Equilibrium Boundary that extrapolates rho from neighbouring nodes (velocity inlet): https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/boundary.py#L300

• Nonequilibrium extrapolation outlet (Krüger): https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/boundary.py#L466

• Nonequilibrium extrapolation Inlet: https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/boundary.py#L518

• Parent class for boundaries using direction attribute: https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/boundary.py#L30

• Equilibrium Extrapolation boundary (Krüger?): https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/boundary.py#L344

• ZeroGradientOutlet (apparently not a great idea to use): https://github.com/MartinKliemank/lettuce/blob/cbeb758247046d6901fa2a6478b7cc0892c5c31e/lettuce/boundary.py#L370

As far as I understand it, lattice.einsum is basically a wrapper that allows to omit the "..." in the dimension-equation (feel free to correct me if it has a greater purpose than that).

• In my opinion leaving out the "..." makes the code less readable, since it doesn't show that there are more dimensions than listed.
• Lattice.einsum also doesn't work if you only want to use it on parts of the domain, which makes it unusable in boundary conditions and the likes.

So I think it would be better to phase out lattice.einsum and return to the bare torch.einsum.