Archive for December, 2007
Entity Crisis has some advice about efficiency.
Snipplr has some code that I might borrow from. I’m just learning python, so I’m eating up all the sample code I can at the moment.
I need to mention a set of “Slides from INF3330 lectures” I came across, by Hans Petter Langtangen, Dept. of Informatics, Univ. of Oslo & Simula Research Laboratory, August 2007. I’ll dig up the URL. There are more than a few code snippets that will come in quite handy – in particular I love the file-to-code idea – defining the actual kernel function in the input file, anyone?
Pylab’s Matplotlib is easy to use, and does a good enough job. However, pyglet sounds amazing. It’s designed as a game engine, but IMO that makes it well suited to an SPH exploratory program.
Now, there is no need to have only one viewer, but in my experience it is a headache to try to integrate more than one windowing/gui system with a model. I’m open to the idea if it turns out no one system will do everything. So, I give you the requirements for the application’s windowing/display system:
- Easy way to save animation as a movie
- Fast enough in real time that the overhead is negligible compared to number crunching of particle interactions
- Free, open source
- Supports 3D, and advanced 2D
- Scientific plotting capabilities
- Minimal need for OS/gui framework specific code
And the candidates:
- Matplotlib (pylab). Requires choosing a backend: GTK, Tk, Wx?
- Some combination of the above
In my last post I touched on the question of data structures for particle systems. It boils down to the question of whether the basic data structure is an individual particle, or a system of particles.
The code you find scattered about the web from the graphics community prefers for the basic data structure to be a particle. I think I’ll find differently when I look at molecular dynamics codes. I’m compiling a list here of the design approach taken by several codes. They are all either production quality publicly available codes, or archetypal design patterns detailed in published work.
First up, NAMD. There is an atom data structure. However neither position or velocity appear to be contained in it. I have to say the comments are pretty bloody sparse!
Next, GADGET uses structures for particles, such that they are referenced in the form particle.position[dimension]. They copy the data (or pointers, not sure) before doing computations though – the force calculation for particle i goes something like
pos = particle(i).position
vel = particle(i).velocity
force = do_something_with(pos,vel)
Not sure why, there is probably a reason.
For now, the group of particles will be one object, with the arrays of various properties as vectors. This leaves scope to introduce multiple groups (e.g. for two mixing fluids). However some of the C functions I need to wrap take arrays of positions, masses etc as arguments, so some scattering will be required. If it turns out the scattering costs too much, I’ll move to a simple r,v type system.
The paper that most influential for me when I built my first toy particle systems is this paper by Andrew Witkin. Now I need to decide whether I am going to follow this pattern, or use another.
The defining features of the Witkin model as I see are:
- Particles are objects(structures), aggregated into a particle system object
- The system state, and its derivatives, are scattered to a ‘state vector’ for some operations
- The solver methods know about the Particle data structure
Now, I’ve dreamed about writing a particle system that uses a completely generic ODE solver, one that only knows about the state vector (just a list of variables) and its derivative, and nothing about a “particle” object. The problem is that any multi-step integration method needs to call the subroutine that calculates the time derivatives of each variable. Unless you want to write this subroutine in a kind of cypher. Imagine s is our state vector – s is position in the x direction, s position in the y direction. It’s much easier to write a position dependent force in terms like particle[i].position rather than having to remember which state vector index corresponds to which dynamical variable.
The point I always come to when considering this (it’s come up a few times now) is this: is it really worth worrying about having a generic ODE solver, given that the solver code is not *that* long, and that the chances of ever wanting to use it for a system with oodles of degrees of freedom that is not a particle system is incredibly slim?
So, it’s OK for the solver to know about particles. Why then do we need to scatter our particle data to a vector. Because it’s faster to do things like increment the whole vector, rather than deal with multiple levels of indirection.
But is this still the case if we have our variables not linked to a particle structure, but instead a particle *system* structure. That is
I still don’t know the answer, still researching. I believe I have C code somewhere that implements both, so a few tests might be in order.
This is my latest google time project. The vision is this – a gui application, with real-time interaction, that displays sph particles interacting in two or three dimensions, rendered in a variety of ways (surface, particles, interpolated..). Splashing, condensation, boiling, sloshing, exploding are all possible. Using a defined netcdf interface, a system state can be loaded from a file produced by a different application (e.g. my number-crunching fortran SPH code), and examined, allowed to evolve while being observed, etc etc.
More mundand objectives are:
1. To build a basic but extensible particle dynamics system with python
2. To convert my old particle system work in C to a library that can be python-wrapped and integrated
3. Use a well separated MVC pattern to enforce the independence of the model and the view.
4. To keep it general enough to be able to use the same framework for basic MD style interactions, and more complex SPH style interactions.
My work on the C particle system hit a bit of a dead end, partly due to the combination of time constraints, and the inability to do anything both quickly and well. Nevertheless the ODE solvers are probably quite suitable to be wrapped and called by python. I’m enforcing a journal for this project, and a slow, steady-as-she-goes development philisophy.
I’m going to break it down into a series of exercises, and post code as I write it. I’ll keep a reasonably up-to-date version on the google code pyticles project.