Photoflow - 
A photorealistic
and physically correct
graphical simulator 
of liquid motion

























Designed and implemented by:

Dimitry Rutshtein

Technion - Israeli Institute Of Technology
Computer Science Department

Nov 2003 - Feb 2005
























This project’s goal is to design and implement a program, capable of modeling a photorealistic and physically correct behavior of a viscous incompressible liquid in an arbitrary, user-defined environment, consisting of solid obstacles, valves, sinks and air. The project doesn’t include (at least at this point) a simulation of moving solid objects and air disturbances.

The project was implemented using Visual C++ 6.0 and IRIT solid modeler.


Keywords: Flow simulation, Navier-Stokes equations, fluid particles, numerical solution, explicit scheme, finite difference method, OOP.







Chapter 1: Designer’s Guide (computational model description):

1)      Previous work

2)      Proposed model outline

3)      Computational grid

4)      Calibration phase

5)      Calculation phase

6)      Balancing phase

7)      Distribution phase

8)      Bounding phase

9)      Summary

10)  Design level extensions

11)  References


Chapter 2: Programmer’s Guide

1)      Compilation

2)      Global structure (main function)

3)      Basic types

4)      Error handlers

5)      Computational grid data structure

6)      Field classes

7)      Cell classes

8)      Configuration module

9)      Polygonal to volumetric conversion (Volumetric module)

10)  Computation and output (Environment module)

11)  Implementation level extensions.


Chapter 3: User’s Guide

1)      Program requirements

2)      Configuration file and command line

3)      Geometric input file formats

4)      Running the program

5)      Error messages

6)      Displaying results

7)      Choosing optimal configuration parameters


Chapter 4: Downloads

1)      Testing examples

2)      The program

3)      Source files

4)      Links




Chapter 1: Computational model


(1.1) Previous work:

A photorealistic simulation of a liquid flow has long been an exciting and a challenging goal. Various early simulations included 2-D wave propagation and waterfalls, which made use of some simple formulas of wave mechanics and free fall motion. However, these simulations can only describe some particular, artificially induced effects in a limited number of predefined simple geometric systems. In order to model a general liquid behavior with maximum range of effects in an arbitrary system, much more complicated methods must be used.

A realistic generic model of liquid motion must also be physically correct (at least to certain degree). This way, the entire range of effects is automatically simulated, just by being incorporated in the laws of physics. There are two basic ways of a physical modeling of a liquid: the numerical solution of a system of full Navier-Stokes Equations (including the continuity equation) for a viscous incompressible fluid motion (hereafter “NSE method”) or a numerical solution of Newton’s Equation of motion for massive particles (hereafter, “particle based method”).

Foster and Metaxas [1] propose calculating a 3-D vector velocity field by solving NSE using a finite difference method in a rectangular grid, with boundary conditions specified explicitly after each computational step. The motion of a liquid can then be approximated by massless marker particles, convected in a calculated velocity field, with 2 or 8-point interpolation.

This method, as any other NSE based one, can account for a full range of liquid volume motion phenomena, such as laminar and turbulent flow, waterfalls, eddies and vortices with relatively coarse spatial and temporal resolution (as compared to a visual resolution). However, any surface effects, such as fine ripples and surface tension are completely absent and must be simulated artificially, which is a hard thing to do in a NS grid. Furthermore, liquid surface’s form and position cannot be derived from velocity field alone and require the use of some tracking methods, of which the abovementioned marker particles method is the most convenient and powerful, but space consuming one. And, obviously, the numerical solution is only possible, if boundary conditions are provided, which can be extremely tricky, especially on a liquid-air surface. Foster and Metaxas [1] suggest a complicated scheme with 64 possible configurations in each surface point, but they fail to address the most problematic aspects of that scheme, thus making it unusable.

On the other hand, particle based methods [4, 8 and 9] rely on interactions between particles, which can be described with ODE, which doesn’t require any boundary conditions. This makes it much more convenient for numerical solutions in an arbitrary system. Also, the attractive and repulsive forces between particles (hereafter, “p-p interactions”) are the ones responsible for surface tension and other surface effects, which are incorporated naturally into particle based solutions. However, to achieve accurate solution inside the liquid volume (vortices, eddies, etc.) one requires a much higher spatial (in terms of particle density) and temporal resolution than the one in NSE. Coupled with the fact that p-p interactions’ calculations tend to be quadratic in particles’ number, this renders the particle method impractical for accurate modeling of a liquid volume.



(1.2) Proposed model outline:

One can notice easily, that the abovementioned two methods exactly complement each other. NSE is perfect for internal volume, but can’t be applied to the surface, whereas the particle based method is preferable on the surface and not inside the liquid.

Therefore, it’s seems natural to consider a hybrid method, based on the following rules:

·        Envelope Condition: any volume of liquid is always separated from the air (or, optionally, from solid boundaries as well) by a thin layer of particles, so that liquid and air never come in direct contact.

·        Surface particles’ velocities and locations are calculated using the particle based method.

·        The field inside a liquid is solved using NSE, with boundary conditions determined by the velocities of the surface particles.

As we can see, this hybrid model takes “the best of both worlds”, thus creating a method which is apriori superior to both previous methods in terms of generality and physical correctness. This idea is distantly related to a hybrid surface method, briefly outlined by Foster and Fedkiw [2], but here the hybridization plays much more fundamental role.

Here we can expect a full range of physical phenomena – both internal and on the surface, without any additional modeling on our part. The higher space and time consumption of the particles is compensated by the fact that for a large continuous volume the surface covers only a small portion of the computation grid. The particles can serve not only for determining boundary conditions, but also for tracking the surface, i.e. serve as marker particles as well. This way the NSE method becomes even simpler. Another advantage of the hybrid model is its ability to set boundary conditions at high resolution, thereby enabling a finer shape of the solids. In pure NSE method, the solids must have the same resolution as the computational grid (which is very coarse, as mentioned above), but here the particles can separate liquid volume from the solid boundary as well, which lifts this restriction. However, implementation of this ability is left for the future developers.

The main problem, of course, is to synchronize two very different models. NSE model uses a static rectangular grid with velocity field fixed on it, while particle method requires dynamic moving particles, not bounded to any grid.

We propose the following model:


·        Initialization:

o       At some initial time t0 the entire environment is represented as a 3-D array of cubic cells of edge size dr (and a volume of dr3). This is the computational grid for NSE, where every cell defines velocity and pressure field (see details in 1.3). In addition, each cell is assigned a type according to its contents:

§         S (Solid)     – the entire cell is filled with solid.

§         I (Input)     – the cell is a part of a valve (liquid input)

§         O (Output) – the cell is a part of a sink (liquid output)

§         L (Liquid)   – the entire cell is filled with liquid

§         B (Border) – separates L cells from E cells and contains at least one particle

§         E (Empty)  – contains only air.

o       The particles exist only inside B cells.

o       The state of the environment at any given time consists of:

§         Current type of every cell

§         Current velocity and pressure at every cell

§         Current location, velocity and mass (size) of every particle

o       We wish to determine the state of the environment at time t1 (t1 > t0). To achieve this, we perform a finite number of evolution (computation) steps, each of which evolves the state of a system by some small time difference dt, which is determined according to several considerations (see 1.4).

Each computation/evolution/time step consists of 5 separate phases.

·        Calibration: the new value of dt is calculated.

·        Calculation: the NSE solution is advanced by dt over the entire liquid volume (all the L-cells) and particle solution is advanced by dt for every existing particle.

·        Balancing: the NSE solution is modified in order to accommodate the continuity equation.

·        Distribution: the entire liquid volume is propagated according to the newly calculated velocity field:

For L and I cells, which contain continuous liquid,  the volume of each cell is temporarily filled with newly formed particles, with velocities provided by NSE solution (initial condition for ODE). Then, all the particles (including the ones in B-cells) are propagated according to their own intrinsic velocities.

·        Bounding: the new position of the liquid determines the new geometry: for example, empty B-cells become E-cells, and so on. The new geometry, in turn, determines the boundary conditions (for example on the L-B border) for the next evolution step.


The following sections describe each phase in detail, as well as the general data structure of the environment.



(1.3) Computational grid:

As mentioned above, the entire environment (the Cartesian 3-D space occupied by solid, liquid and air) is partitioned into cubical cells of variable types with edge size dr. The cells are cubical for the sake of simplicity. Since the cubes are relatively small, the number of cells varies in each direction and temporal resolution is determined according to the coarsest spatial resolution, this simplification doesn’t affect the generality of the problem or even the performance.


·        The environment is a box of Dx cells in x direction, Dy – in y direction and Dz – in z direction. Each cell is assigned 3 integer indices (X, Y, Z), which determine its position in the grid.


(1)   0 ≤ X < Dx     0 ≤ Y < Dy     0 ≤ Z < Dz

(2)   A point (x, y, z) belongs to the cell (X, Y, Z) if and only if:

X∙dr ≤ x < (X+1)∙dr    and    Y∙dr ≤ y < (Y+1)∙dr    and    Z∙dr ≤ z < (Z+1)∙dr


·        For each cell, its origin point is defined as the one with the smallest x, y and z values still inside the cell:


(3)  O(X,Y,Z) = (X∙dr, Y∙dr, Z∙dr)


·        Each cell has 6 faces, aligned with the Cartesian axes – 3 backward (closer to the (0,0,0) cell) ones and 3 forward ones. Thus the origin of the cell is located on the intersection of the backward faces. By definition, backward faces belong to the cell and forward faces belong to its forward neighbors.


·        Scalar pressure field P(x, y, z) is defined in the center of each cell. In any arbitrary point it can be determined by interpolation of the closest 8 cell centers or, for a point in a center of a cell’s face, by interpolation of the 2 centers of adjacent cells.


·        Vector velocity field v(x, y, z) is separated into 3 scalar fields: u(x, y, z), v(x, y, z) and w(x, y, z) – velocity components in the x, y and z directions respectively. Each field is defined in the center of a cell’s face, which is perpendicular to the proper axis (see Figure 1). Each component can be determined at an arbitrary point by a linear interpolation of the 8 closest face centers (8-point interpolation) or by interpolation of the 2 nearest face centers - of the cell in which the point is located (2-point interpolation).




(1.4) Calibration phase:

The time step dt must obey certain rules, if a proper and practical solution is to be received:

·        It should be small enough to ensure correctness, because each evolution step creates an error of O(dt3) in fluid’s position and these errors are cumulative.

·        On the other hand, it should be large enough to reduce the number of computation steps, because the time required for each step is virtually independent of dt.

·        It must also obey the convection stability criteria: dt ≤ dr / max(u, v, w). This will ensure the numerical correctness of the NSE solution.

To accommodate the first two rules, there are 2 user-defined parameters: dtmax and dtmin. Time step must always be between these values, but in any case – as large as rule 3 permits:


(4)   dt = min(dtmax ,  dr / max(u, v, w) )


The second rule is accommodated automatically, by truncating the velocity field to dr / dtmin

Although such truncation may create non-physical behavior of the liquid, the effect remains local and visually unnoticeable and occurs only in extreme cases. In any case, it remains in complete control of a user, who determines dtmin.


Complexity: This phase requires exactly one pass over the entire grid, for calculating max(u, v, w).



(1.5) Calculation phase:













Where: u = u(x, y, z, t), v = v(x, y, z, t), w = w(x, y, z, t)       – velocity field components

P = P(x, y, z, t)                                                           – pressure field

g = g(x, y, z, t)                                                            – free force acceleration (usually a constant vector)

r = const                                                                     – density

n = const                                                                     – viscosity


These are the Navier-Stokes equations, which describe the motion of an incompressible viscous fluid in a 3-D Cartesian environment. (5) and (6) are identical, but (6) is written in a more explicit notation.

Despite such impressive exterior, their meaning is actually almost trivial. They simply state that for each reference point (x, y, z) inside the liquid, there are 8 forces which act on a differential volume dV of liquid at that point, causing it to move with acceleration a = F/m, according to the second law of mechanics (Newton’s equation of motion), where F is the sum of all 8 forces. Since each point contains the same mass m = rdV, we can assume m = 1 and use the terms force and acceleration interchangeably for the sake of convenience. Each force can be analyzed from a microscopic (particle based) and macroscopic (differential volume based) point of view. The NSE are based on the macroscopic formulation. These forces (accelerations) are (in the order of appearance):

1) Free force (independent of the liquid external force) a1 = g. Usually it’s just gravity: g = (0, 0, -9.8 m/sec).

2) Pressure induced force. On a microscopic level, pressure is caused by repulsive forces between particles which are too close to each other. On a macroscopic level, the combined forces of many particles can form a pressure difference on each side of the reference point, pushing its volume to the side with less pressure.

3-5) Viscosity or diffusion force. On a microscopic level the attractive forces between the particles cause them to perform inelastic collisions with each other. Particles loose energy in these collisions and tend to move with the average velocity of their neighbors. On the macroscopic level this averaging will cause the liquid’s velocity to change no more than linearly within the volume. The second derivative, which represents the lowest order of deviation from linearity, becomes the measure of the linearization force. n is the viscosity coefficient – the measure of particles’ tendency to stick to each other, i.e. the non-elasticity of collisions. Of course, there are 3 diffusion forces – one for each dimension.

6-8) Convection. As stated above, NS equations are defined relatively to a certain static reference point. However, the liquid performs convection - constantly moves passed this point at alternating velocity. The old liquid leaves the reference point and is replaced by a new liquid, which moves at a different velocity. Since all differential liquid volumes are assumed to be identical, this creates the illusion of acceleration at the reference point, because without tracing each particle (as particles don’t exist in macroscopic description) we can just as well say that the volume never left the reference point, but simply changed its velocity for some reason. Therefore, such formulation generates a fictitious (imaginary) convection force. Again, there are 3 such forces – one for each dimension.


There are several ways of solving these equations numerically and we choose the simplest one – the explicit FTCS (Forward Time, Center Space) scheme of finite difference method. This means that we simply replace each term on the right (spatial) side of the equations with its second order approximation (according to the Taylor series) and on the left (temporal) side – with its first order approximation (there aren’t enough initial conditions for second order approximation anyway).





















Where: dt – as determined in the calibration phase and dx = dy = dz = dr.


These are the finite difference approximations of every term in (6.a). For (6.b) and (6.c) the approximations are basically the same with u replaced with v and w respectively. By combining all these terms we receive 3 algebraic equations with 3 unknowns: u(x, y, z, t+dt), v(x, y, z, t+dt), w(x, y, z, t+dt), which is exactly what we want to find. Their solution is straightforward. Everything else is already contained in the grid, although some values, such as v(x, y, z) must be interpolated from the nearest values (in this case – an average of 4 nearest points).

By applying appropriate equation at the center of each liquid cell’s face, we receive the velocity field for the next iteration up to the numerical error of O(dt2+dt∙dr2). By choosing dt ~ dr2 , the error is reduced to O(dt2).


The velocities of surface particles, however, must be solved in a different way, because NSE cannot be applied to the surface – this would require velocity values from adjacent empty cells, where liquid’s velocity is obviously undefined for the lack of such. However, the underlying principle is the same as in NSE – the second law of Newton [9]. For each particle i of n particles in a liquid:






This is, actually, the same NSE. The first term on the right hand side is gravity, the second – total repulsive force on a particle, which is the microscopic equivalence of the pressure gradient and the last term is the total attraction force, which contributes to the viscous terms of NSE. The convection term is missing, because here the reference point is fixed on the particle and a particle is always motionless in its reference frame, hence no convection. Since the microscopic forces can be presented in vector form, there is no dimensional separation and we get an ODE, instead of a PDE.

Of course, (8) is very difficult and costly to calculate – mainly because it’s linear in particles’ number and also because the repulsion and attraction forces are so complicated, that no exact formula for them exists yet, although there are some very accurate simplified models [9].

We choose an approximate model, which is nevertheless effective, by summing the forces only on the nearest neighbors (in adjacent cells). This is a good approximation, because p-p interactions decay quickly with distance. Also, we assume a = 0 and b = const (user defined). Together with the rejection model, implemented in the distribution phase (see 1.7) this gives a fair first order approximation of the true force generating potential (see figure 2). It’s enough to recreate the basic behavior of particle interactions. Furthermore, we solve (8) in the same resolution as the NSE (both spatial and temporal). This means that we use the same dt as in NSE solution (although physically accurate particle based model generally requires a smaller one) and regard adjacent cells as homogeneous collective particles, rather than a collection of different particles. This reduces the calculations’ complexity from O(n) to O(1) for each particle (linear total complexity) and saves us the trouble of modifying dt in the middle of a time step in order to accommodate the higher resolution of particles.

All these simplifications seem like a big stretch and they are, but amazingly enough, the model is still good enough to produce realistic results. This is because in a turbulent flow, the surface particles exist only for a very short time as separate entities (just a few time steps), before they are reabsorbed by the liquid and the new particles that replace them receive the initial conditions from the accurate NSE solution. Therefore, the errors produced by all the simplifications don’t have a chance to build up and affect the overall solution. As for the laminar (non-turbulent) flow, the errors build up very slowly, because the p-p interactions have negligible effect in this case, relatively to gravity and inertia, which are represented accurately even in this approximate model. At the same time, the surface effects are still represented well, especially the surface tension, which is a direct consequence of the attractive force.

Complexity: As in the previous phase, this one requires exactly one pass over the entire grid, although for B-cells (containing particles) the calculation complexity is linear in particles’ number (in the simplified model).

















(1.6) Balancing phase:

According to problem definition, we wish to model a flow of an incompressible liquid, but in reality no fluid is absolutely incompressible.

The fluid, which obeys (5), flows at different velocities in different locations, so it’s entirely possible that one point will suddenly have higher or lower density than the other. This is where the pressure comes in. Pressure is a function of density: P = P(r), so higher density almost immediately creates higher pressure at the same point. This, in turn creates a strong pressure gradient, which according to (5) generates a negative force in order to stop the incoming liquid and even reverse it. The more incompressible the fluid is, the stronger is the gradient and the resulting reversed acceleration. In nearly incompressible fluids, the incoming volume is reversed almost immediately and returns to the points from which it came, creating increased pressure (and another acceleration) there as well. This process of passing the pressure to neighboring points continues throughout the volume, creating a pressure wave, also called the sound, which expands in all direction at great speed, until reaching the surface, causing it to expand a little, thereby reducing overall density and restoring the pressure balance. This process is called relaxation. By pushing this reasoning to the limit we receive an abstract model of absolutely incompressible fluid, in which the reversal is instantaneous, the accelerations are infinitely large and the pressure waves expand at infinite velocity. This abstract model works very well for many liquids, providing that the typical velocities are much smaller than the speed of sound. For example, in water (which is usually regarded as incompressible) a 1% increase in density creates a pressure difference of 200 atmospheres and the resulting pressure wave expands at over 5000 km/h (3000 mph). In most practical environments this can be regarded as infinitely large. Note however, that incompressibility condition limits the maximum size of the simulation environment by L << vsounddt, otherwise the pressure wave won’t have enough time to propagate throughout the volume within one time step, creating compressed areas in the next time step.

In numerical solution we can’t rely on NSE to uphold absolute incompressibility, for it would require dt = 0 to model infinite accelerations. Therefore, we must balance out the pressure instantly, so that in the next time step the pressure waves will be already long gone and the density will be restored to initial constant value, allowing the use of (5), which only works for constant density. To achieve this, we first use the iterative pressure refinement algorithm, proposed by Foster and Metaxas [1] and then supplement it by explicit particle oriented relaxation.

Each grid cell defines 6 velocity components: u0, v0, w0 in the centers of backward faces and u1, v1, w1 in the centers of forward faces. Assuming uniform velocity throughout the face, the volumes passing through each face are:


(9)   dr2u0dt;   dr2u1dt;   dr2v0dt;   dr2v1dt;   dr2w0dt;   dr2w1dt;


The volume increase inside the cell is therefore:


(10)  DV = dr2dt(u0-u1+v0-v1+w0–w1 )


Without the loss of generality, we assume that DV is positive.

We need to “push” all that excessive volume away to adjacent cells, to relax this one. The “push” must be equal in every direction, because the pressure wave has radial symmetry. Therefore, each velocity component is changed by:


(11) dv = DV/(6dr2dt)


The change is in outward direction (negative for backward faces and positive for forward).

Now, DV= 0 for this cell, but since every component is shared by 2 cells, it may not be true for adjacent cells.

Therefore, we must apply (11) to every cell in the grid. Unfortunately, since the change is symmetric, the already balanced cells may become unbalanced again in the course of iteration. However, the DV will be smaller this time, because only a portion of displaced volume returns back to the cell, so by repeating the iteration several times we can reduce DV in each cell as much as we want (to some acceptable small value).

Foster and Metaxas prescribe very small acceptable values for DV, because large ones can cause significant volume loss after many time steps. However, in our model we perform explicit mass tracking (see 1.7), so we don’t have to perform exact relaxation – only enough to balance the liquid in general, to facilitate further relaxation during the distribution phase.

Complexity: The required number of iterations varies greatly, but since the calculations here are very simple, this phase takes about the same time on average, as the calculation phase.



(1.7) Distribution phase:

By now the velocity field is finally (more or less) set and we can start preparing for the next time-step. But for that we need new boundary conditions and before these conditions can be deduced, the new boundary itself has to be located. This is achieved by propagating the liquid, using the current velocity field.

For each point particle at r(t), moving with some varying velocity v(t) the location after the time step dt can be approximated by expansion in power series:


(12)  r(t+dt) = r(t) + v(r(t),t)∙dt + 0.5∙a(r(t),t)∙dt2 + O(dt3)


Since the acceleration a(r(t),t) has been calculated during calculation phase with O(dt) error (see 1.5), there is no point expanding the series beyond the second order – that will no longer reduce the error.

A common practice is to simplify (12) to: r(t+dt) = r(t) + v(t)∙dt  or  r(t+dt) = r(t) + v(t+dt)∙dt

Both are first order approximations with O(dt2) error, which doesn’t take advantage of the O(dt3) accuracy of the NSE solution and that’s obviously an unnecessary loss, because the acceleration value is readily available.

Note: we don’t even have to store the acceleration, because we can recalculate it from the velocity difference: a(r,t) = ((r,t+dt)–v(r,t))/dt + O(dt), so we can rewrite (12) as:


(13)  r(t+dt) = r(t) + 0.5∙(v(r(t),t)+v(r(t),t+dt))∙dt + O(dt3)


This takes care of the surface particles, but the liquid volume also has to be propagated. For this purpose, we generate a set of particles for every filled cell, uniformly distributed throughout their volumes and propagate them according to (12), with v(r,t) and a(r,t), as provided by interpolation of the existing NSE solution. The particles which end up inside L-cells are reabsorbed by the liquid and others remain separate entities for now.

Note: the particles are generated and propagated for one cell at a time, so the additional particles create very little space complexity overhead, as most of them are reabsorbed (destroyed) immediately.

Mass conservation is very important physical quality in classical mechanics. Foster and Metaxas [1] rely on the pressure refinement phase for maintaining the mass conservation automatically, but we choose to maintain it explicitly.

For each L-cell, we keep the mass of a liquid, currently inside the cell. Each particle also carries some constant mass (can be different for every particle). Each particle leaving an L-cell takes a portion of cell’s mass with it and donates it to the L-cell in which it eventually turns up and absorbed (the absorption can happen immediately or after several time-steps). So, by adding up all the masses inside L-cells and the intrinsic masses of currently existing particles we get the total mass of the liquid, which is guaranteed to remain constant, disregarding valves and sinks. In I and O cells mass obviously must change. Input cells increase their mass by a predefined value before each step (as if the liquid is “created” inside the valve) and Output cells absorb incoming particles, ignoring their mass, as if they are “swallowed” by the sink.

The L-cells can not only absorb particles, but also to reject and return them to their places of origin, if the cell is already full. This technique is another kind of relaxation, only it doesn’t fix the problem of density change – it tries to prevent the problem from occurring in the first place. The rejection technique isn’t perfect in the current implementation and should be extended (see 1.10), but together with balancing phase and proper morphing (see 1.8) it can reduce density variations to an arbitrary small value. In our experiments, even with allowed DV as high as 10%, the average density variations were no more than 3-5% in a turbulent flow (in which it’s very hard to notice density variations anyway) and less than 1% in laminar flow. In standing water, there were no density variations detected. In any case, the mass conservation was always perfect.

Explicit mass tracking also enables us to design a simple and convenient model for pressure. Foster and Metaxas [1] use pressure refinement to update the pressure field. This isn’t the best technique, mainly because it creates a pressure buildup delay – the pressure can still be high even when the extra volume, which created it, is already gone.

Since pressure is a function of density and we can calculate the latter for each cell using mass tracking, calculating pressure is easy. Assuming small enough variations in density, we can approximate P(r) even without knowing what the function is (which is actually extremely complicated):


(14)    P(r) = P0 + a(r /r0 – 1) + O( (r /r0 – 1)2 )

            where:   P0 – atmospheric pressure

                          r0 – normal density

  a = dP(r0)/dr – user defined coefficient, which describes pressure sensitivity to density

     variation (which is, of course, different for every liquid).


This way, the distribution phase automatically sets the new pressure field values, to be used in the next time-step’s calculation phase. With such definition, the pressure always reflects the current liquid distribution – no delay effects. Moreover, it doesn’t even matter what the exact pressure is – only the general trend is important. This is because pressure effects are mostly accounted for during balancing and rejection, and all we need to do now is to account for some small leftover variations, caused by imperfections in balancing and rejection models. This will be done during the next calculation phase and further reduce the errors of balancing and rejection, contributing to even more accurate relaxation.



(1.8) Bounding phase:

With all the liquid in proper place now, we can deduce the contents of each cell. Basically we’re talking about E (Empty), B (Border) and L (Liquid) cells, which can “morph” into each other (all other cell types are static).

Cells with zero mass are labeled E and with non-zero mass – L, unless they are adjacent to E cells, in which case they become B. These are the natural morphs. To this we can add other rules for forced morphs. For example, a cell without adjacent cells can still become B, if it contains too little mass – below some predefined threshold (but still above zero). This is actually an existing physical effect – in a region of small pressure bubbles start to appear, even if it’s deep under water. Another forced morph – B-cell with too much mass can become L-cell and force its E neighbors to become B cells. This is only necessary if balancing and rejection methods are too simple and can produce noticeable deviation from standard density – this condition will prevent this.

The L-cells morph into B-cells by having their volumes converted into a set of particles, just as in distribution phase, but this time the conversion is permanent (at least for this time step). Likewise, B-cells morph into L-cells by absorbing all their particles and becoming uniformly filled. This restores the Envelope Condition (see 1.2), broken by the distribution phase.

With all these conditions, one iteration over every cell may not be enough, because each morph can affect all adjacent cells, even those which we have already passed. Therefore, bounding iteration must be performed several times, until there are no more changes.

Complexity: In all our experiments, the number of iterations in this phase has never exceeded 4.

All that’s left now is to calculate boundary velocities for the new liquid surface. In hybrid model this is very simple and is performed simultaneously with the morphs. There are two possible places for setting the boundary conditions: on the faces between L and other cells or on the opposite sides of these faces (one cell away from the L boundary). The first case is preferable when the method of calculating boundary conditions is accurate and we trust it more than NSE solution. Obviously, there is no need to solve NSE on these faces. The second case is useful when boundary conditions are accurate enough for describing surface evolution, but produce error which is greater than the maximum allowed divergence in relaxation, which may prevent the balancing phase from converging. By relaxing the boundary conditions in this way, we allow relaxation to succeed, without significantly affecting the overall evolution.

For liquid boundary with solids, valves and sinks we employ the first method, using predefined values. For B cells, the boundary conditions are set on the far side (adjacent to air), using the weighted average velocity of all the particles in each cell. The near side retains the values, calculated so far.

As for the E cells, their values are of absolutely no importance, because particles ignore field values and L-cells never come in direct contact with the air. For the sake of elegancy we assign them zero velocity and atmospheric pressure.


This concludes the evolution step and the environment is ready for output in any desired format and/or for the next time-step.



(1.9) Summary:

We’ve presented a hybrid method of photorealistic evolution of viscous incompressible liquid, which is a combination of numerically solving the full Navier-Stokes equations in a uniform Cartesian grid and numerically solving the Newton’s equation of motion for individual particles. The first method is applied to the internal volume and the second one – to the surface. Both methods work in a tight cooperation, with the boundary conditions for the first method set by second one (plus some predefined values for static boundaries) and the initial conditions for the second method are set by the first. Correct mass distribution and conservation is ensured by explicit mass tracking.



(1.10) Design level extensions:

1) The most obvious extension is upgrading the rejection technique. Rejecting a particle for a filled cell isn’t as simple as it sounds, because the current cell’s mass can be used as a measure. The current mass may partially or even fully leave the cell as distribution continues, so the particle will be rejected for nothing and vice versa. Currently we use velocity components of a cell to project the eventual filling (simply by rejecting particles which come in “against the current”). Although fairly successful, this method is obviously prone to errors. We suggest the following upgrade:

- Distribution phase is separated into 2 phases: initial and final.

- Initial distribution is like the current one, except that each particle is always accepted (or allowed to pass through) by any cell (except S, of course).

- Final distribution basically replaces the balancing phase (which may well no longer be required). It performs direct relaxation of the mass distribution: the particles which are too close to each other (inside each other, for example) are merged into one particle with combined mass and weighted average velocity. Too large particles are separated into evenly spaced smaller particles moving at the original velocity. The size of a particle is derived from incompressibility condition (i.e. V(m)=m/r0) and the assumed shape (cubical or spherical, for example). Some particles may leave the cell’s boundaries in the process and will have to be passed along to adjacent cells. The process continues throughout the environment, until all the particles are evenly spaced and the surface particles are of more or less equal sizes. Each L-cell is treated here as a single, large, but undividable (within a single cell) particle which velocity is the average of all its 6 components. Note that components can change now. This redistribution process is reminiscent of garbage collection, algorithms for which may be adaptable to this problem. If so, this may greatly increase efficiency of this phase.


2) Even in initial distribution, L cells can be treated as single particles, thus increasing the efficiency of the initial distribution as well. However, this would reduce velocity accuracy to a 2-point interpolation, which is known to exhibit inaccurate behavior sometimes. But, with the new rejection model this might still work nicely.


3) Particle interactions can be modeled with higher accuracy. For example, making attraction forces proportional to the distance between particles (a second order approximation). Repulsive forces can be added, but with the new rejection model, this would probably be a redundancy. Another idea is to solve particles in a different (smaller) time-step, but this also may be an unnecessary complication.


4) In any pure NSE method, the boundary conditions must be set at relatively coarse computation grid resolution. Therefore, solid objects must have the same resolution, which is not so good. Hybrid model allows much higher resolution of solids. To achieve this, the Envelope Condition must be extended, so that the particles would separate liquid from solids as well. In practical terms, this would add a new heterogeneous cell type – partially filled with solid. Unlike B-cells which are binary heterogeneous (liquid + air), this one is at least trine (liquid + air + solid). This is probably the main complication of the extension – to implement particle acceptance and rejection for this type. In addition, other cell types can be introduced: partially filled with valve and sink and the combination of them all.

With this extension, other surface effects can easily be reproduced, namely friction with solids (which does exist in the current implementation, but it’s not very flexible or accurate).


5) With the flow reaching an optimum performance in terms of realism, the next logical step is to add dynamic solid objects. Current design allows only for static S-cells (including the partially solid ones from (extension 4), but that can be changed, in principle. The most simple ones are small buoyant objects (no bigger than a computational cell) which don’t have noticeable effect on the liquid. They are simply propagated according to local velocity field after each step and included into the environment only at output level. This is enough to simulate various junk, carried by a strong current. Alternatively, the objects can carry their own velocity and update it according to local liquid acceleration (which is derived from the initial and final field values) as well as some other forces, like Archimedes force (but not gravity – that’s already included). This can simulate many interesting situations, like coffee grains descending onto the cup’s bottom, leaves floating and gliding on the surface and so on.

Introducing large objects, which displace large volumes and affect the flow is a tricky business, but with enough time and effort it’s also possible.


6) Another idea is to change the free force to include something besides gravity. Likely candidates are the centrifugal force, Coriolis force, wind blows and walls vibrations. Unfortunately, there is no easy way to define such forces on a user level.


7) The amount of possible extensions is practically unlimited, so there is no point to describe all of them. Here’s just a handful of examples: adding different types of liquids, simulating chemical reactions, heat conduction, boiling and freezing, explosions and shock waves (which require compressible flow) and so on.



(1.11) References:


[1] – “Realistic Animation of Liquids”, by Nick Foster and Dimitri Metaxas, Center for Human Modeling

and Simulation, University of Pennsylvania, Philadelphia, June 3, 1996.


[2] – “Practical Animation of Liquids”, by Nick Foster and Ronald Fedkiw.


[3] – “Controlling Fluid Animation”, by Nick Foster and Dimitris Metaxas, Center for Human Modeling

and Simulation, University of Pennsylvania, Philadelphia.


[4] – “Particle-Based Fluid Simulation for Interactive Applications”, by Matthias Müller, David

Charypar and Markus Gross, Department of Computer Science, Federal Institute of

Technology Zürich (ETHZ), Switzerland.


[5] – “Modeling Realistic Water and Fire”, by Sérgio Leal, Institute of Computer Graphic and Vision,

Technical University of Graz, Austria.


[6] – “IRIT Solid modeler” lecture slides, Center for Graphics and Geometric Computing, Technion


[7] – “Programmers' Manual - IRIT“,


[8] – “Using Charged Particle Systems to Simulate Fluid Flow”, by Jeremy S. Debonet & Chris 

Stauffer, May 2003.


[9] – “Computer Graphics”, by Alla Sheffer and Dani Brunstein, October 2002.




Chapter 2: Programmer’s Guide


(2.1) Compilation:

  • The program can be compiled using Microsoft Visual C++ 6.0 or higher with STL’s ‘list’, ‘iostream’ and ‘fstream’ modules.
  • Compilation also requires a fully installed IRIT solid modeler. Visit to find out how to download and install IRIT (about 100 MB of disk space is required).
  • The VCPP project must have IRIT’s ‘include’ directory among its default include directories and IRIT’s ‘lib’ directory among its default library directories.
  • The following libraries must be added to the linker’s library list (just copy-paste the following list):

bool_lib.lib cagd_lib.lib geom_lib.lib mvar_lib.lib prsr_lib.lib symb_lib.lib trim_lib.lib triv_lib.lib trng_lib.lib user_lib.lib xtra_lib.lib misc_lib.lib ws2_32.lib

  • And finally, the project’s home directory and workspace must contain all the following source files:





Array3D class template – a 3-D array of arbitrary elements


Index and Float type definitions, constants and auxiliary functions

BorderCell.h; BorderCell.cpp

BorderCell and SubGrid classes module

Cell.h; Cell.cpp

Cell class module and CellType definition

Config.h; Config.cpp

Configuration class module

ConfigItems.h; ConfigItems.cpp

GenericItem, ConfigComment and ConfigItem classes module

Curve.h; Curve.cpp

Curve class module

EmptyCell.h; EmptyCell.cpp

EmptyCell class module

Environment.h; Environment.cpp

Environment class module


All error handlers definition

Field.h; Field.cpp

Field and Velocity classes module, EvolutionPhase definition

HetCell.h; HetCell.cpp

HetCell class module

HomCell.h; HomCell.cpp

HomCell class module

InputCell.h; InputCell.cpp

InputCell class module

LiquidCell.h; LiquidCell.cpp

LiquidCell class module

Mstring.h; Mstring.cpp

String class module,  StrList definition

OutputCell.h; OutputCell.cpp

OutputCell class module

Particle.h; Particle.h

Particle and ParticleList classes module

SolidCell.h; SolidCell.cpp

SolidCell class module

StdAfx.h; StdAfx.cpp

Files required by VCPP environment


Vector class module, Array3Vector definition

Volumetric.h; Volumetric.cpp

Volumetric and VoxelStatus classes module


Main function


  • The 41 source files occupy 100 KB and contain 4000 lines of code. The resulting executable is almost 900 KB under Windows XP. Total compilation time is 40 seconds on a Pentium III 800 machine.



(2.2) Global structure (main function):

Main module creates a single global instance of Configuration class and a single local instance of the Environment class. Configuration is accessible from every point in the program, which makes it very convenient to use. Both instances exist during almost the entire program’s run time.

The main function does almost nothing, in accordance with the OOP paradigm.

First it loads a new configuration from file and command line (by calling Redefine member functions, see 2.8) and tests parameters’ validity (Validate member function). This order of loading supports the configuration hierarchy, described in 3.2.

And only then it creates the instance of Environment. Every other functionality is incorporated into Environment’s constructor, so the main function ends here.

The following sections describe all of the program’s modules in a bottom-up order (from basic to general ones).



(2.3) Basic types:

  1. Index = int
  2. Float = double
  3. String – similar to MFC’s CString class, but with a more powerful conversion mechanism – can convert various types to String and be converted to various types (which will be very important).
  4. Vector – can be used as an array or a struct of 3 Floats: x, y and z. Implements many 3-D vector operations, including dot (*) and cross (&) products.
  5. Curve – defines a continuous scalar function f(t) for any positive t. Currently defines only constant functions, but can be easily upgraded to an arbitrary function (see 2.11). Used by InputCell and OutputCell.
  6. Array3D – template for a 3-D array of arbitrary elements of type T. T can be any assignable type for which operators “istream >> T” and “ostream << T” are defined. Operator >> must be able to ignore any white spaces before and after T in the stream. 3-D array is actually implemented as 1-D array, with operator(i, j, k) and STL–compliant iterator to provide fast access to any item. It’s used as a building foundation of SubGrid and Volumetric classes.



(2.4) Error handlers:

Each time an error occurs, the program throws an exception – an object (handler) of type Error or its derived type. There is a derived type for each possible error. Each error object can print an error message, specific for its type, so the main function, which ultimately receives the exception, has only to activate the error handler (via its abstract print function) and the appropriate error message will be printed automatically.



(2.5) Computational grid data structure:


·        Cell – abstract class which defines a cubical region in a 3-D array of such cubes, with specified origin point.

·        FieldElement = (Float∙Float) – defines a scalar field value at certain (unspecified) point in space for 2 different times (initial and final). The value is hidden and can be viewed or modified using virtual GetValue and SetValue functions (see below).

·        Particle – a point object with specified mass, current position and velocity.


Computation grid is implemented as an interlinked 2-way 3-D list-like structure of Cells and Fields. Each cell points to 6 Field Elements, called adjacent fields and each Field Element points back to 2 cells to which it’s adjacent (adjacent cells). Edge Field Elements have only one adjacent cell and their other side points to NULL. The 6 adjacent Fields provide appropriate velocity field values at location corresponding to their position in the structure (see figure 3).



This interlinked structure (with Cells connected via Fields and vise versa) is required for dynamic morphing of cells. There are various cell types (see 2.7), with very different internal structure. Instead of using a union, which is unsafe and inflexible (and violates OOP paradigm), we use abstract classes, in accordance with the OOP paradigm. This means that cells must replace each other in the grid, to allow dynamic type changes at each point (a crucial quality during surface evolution process). If Fields were located inside the Cell (as data members), they’d have to be copied during each replacement. Furthermore, it’s important for a Field not to belong to any specific cell, as it has 2 adjacent cells which it must treat symmetrically (as explained below). For these and other reasons we choose to have a static 3-D net of Field Elements, to which different types of Cells can be easily attached and removed. Since Fields have to “know” their adjacent cells anyway, such approach doesn’t require any additional memory. By traveling along the pointers we can arrive to any Cell or Field Element in the grid.

This is the main data structure, which serves as the foundation for the Environment class.



(2.6) Field classes:

Field Elements define vector velocity field for the entire grid and individual velocities of particles. Both uses require some similar, but also some different operations. Furthermore, we need to make transition from scalar to vector fields. Therefore we define class structure, as described in figure 4.



Field-Element takes care of the speed truncation during SetValue, which upholds stability criteria (see 1.4) and the GetValue function, which returns different values for different evolution phases (see below). This behavior is common to both grid and particles.

Grid-Field-Element is designed to enforce two principles:

  • Locality Principle – a grid velocity component can be modified only from cells adjacent to it.
  • Authority Principle – if two cells are adjacent to a field element, only the one with higher authority can modify the mutual component (if authorities are equal – both can).

These principles are very important for implementation, as will be explained below. These principles should be maintained within the Field Elements, according to OOP paradigm, to prevent programming errors regarding these principles.

  • Grid-Field-Element keeps 2 pointers to each of its adjacent cells (Field Element doesn’t have them) and overrides SetValue function, which now checks whether the modification request is valid (invalid requests are ignored).
  • Particle-Field-Element doesn’t define any new behavior.
  • Particle-Velocity is just a triplet of Particle-Field-Elements, which together define a single vector velocity.
  • Grid-Tracer is the structural basis of the Cell, which connects it to 3 adjacent Grid-Field-Elements. Each Cell has 2 Tracers – one for backward faces and one for forward faces.


The convenience of black-box field values is in their phase-oriented flexibility. Since there are in effect 2 values for every field, the initial and final (for the current time step), every time we use a grid value we have to choose, which value to read or write to. The Field’s GetValue and SetValue functions allow the programmer to ignore this problem entirely – the Field will automatically reference the proper value, according to the current evolution phase (as defined by its static part, see also 2.10). The following table summarizes this behavior.







Swap initial and final values




Calculate new final values (next step’s initial ones)




Fix final values




Final values may change due to current boundary conditions




Calculate boundary conditions for the next step



As will become evident below, this behavior provides correct access to filed values.



(2.7) Cell classes:

Currently, there are 6 specific cell types, separated into Homogeneous and Heterogeneous groups. HomCells don’t have internal structure, i.e. they represent computation level resolution and HetCells are separated to many smaller sub-cells at visualization level resolution. First resolution is used for NSE solution and the second one – for output. This means that a sub-cell should be about the size of a pixel to get adequate visual results.

Obviously, HetCells occupy much more space, so there should be as little of them as possible in the grid.

For each Cell’s subclass we define a numeric type and an authority level for grid field modifications (see 2.6):
































Every sub-class must define it own type and authority, and implement several predefined abstract functions – one for each evolution phase, plus some grid-oriented functions.

To simplify implementation and to improve performance, it is important to make the implementation as universal and natural as possible, i.e. not to rely on the current state of the environment and to avoid type specific conditions as much as possible.

Using OOP paradigm and the Locality Principle (see 2.6), this can be achieved relatively easy. Locality Principle allows each cell type to implement an independent set of evolution functions and OOP based program will automatically call an appropriate set for every cell in the grid, using virtual tables. This means that grid-level algorithms don’t have to care about underlying geometry, which of course makes them simple, homogeneous and fast (see 2.10).

Chapter 1 defines 5 evolution phases: calibration, calculation, balancing, distribution and bounding. Each cell type must implement 5 evolution functions – one for each phase. The tables below present the implementation principles for each type and phase.


Phase 1: Calibration

This phase only sets up field values for the coming time step, so the behavior is the same for every cell type (maximum field value is calculated on the Environment level). Therefore, the ‘Calibrate’ function is implemented once in the root class and not in every sub class. It simply calls Calibrate function for every adjacent field elements (or, as an optimization, just the forward ones). The only exception is the Border cell – it also activates Calibrate function of every particle inside it (recall that particles contain Field elements inside them too).


Phase 2: Calculation

The NSE solution is moved forward in every cell, containing homogeneous liquid and particles’ velocities are updated too, wherever they are (currently, only in B-cells).


Description of ‘Calculate’ function


No change (remains 0).


Set all 6 components to a predefined value, as returned by Curves and reset inflow components (which are negative for positive faces and vice versa) to 0, to prevent suction. Increase mass by intake capacity (calculated from inflow direction’s size).


Set all 6 components to a predefined value, as returned by Curves and reset outflow components (which are positive for positive faces and vice versa) to 0, to prevent ejection.


Evolve using NSE, as specified in 1.5.


Calculate forces which work on every particle and calculate their new velocities, as described in 1.5.


Set all 6 components to 0 (just in case).


Phase 3: Balancing

Returns true iff no balancing was needed.


Description of ‘Balance’ function


Return ‘true’.


Return ‘true’.


Return ‘true’.


Adjust current field values (calculated in Calculation), as specified in 1.6.


Return ‘true’.


Return ‘true’.


Phase 4: Distribution

The “donor” cell sends its particle, which has just left its boundaries (but still formally belongs to this cell) to the adjacent “acceptor” cell – the one which is located on particle’s trajectory.

Grid function ‘Accept’ is an important part of this phase and it is described in the following table.


Description of ‘Distribute’ function




Convert all volume to particles, propagate them and use Accept to locate new hosts for them.




Convert all volume to particles, propagate them and use Accept to locate new hosts for them.


Propagate incorporated particles, use Accept to locate new hosts for them and reset sub-grid.




Accepting a particle:

Reject - return the particle to the donor and prevent it from entering again, by adjusting its velocity accordingly.

Swallow = destroy. Absorb – add its mass to the cell’s and only then destroy.


Description of ‘Accept’ function


Reject always


Absorb particles coming from donors of higher authority, reject otherwise.


Swallow always


If donor’s authority is at least L or field values are compatible with particle’s direction – absorb, otherwise reject.


Keep the particle or allow it to pass through, to its final destination.


Allow to pass to final destination, or if that’s the destination – morph to B and keep.


Phase 5: Bounding

If any morphing occurs, the current cell deletes itself (performs suicide) and Bound returns a pointer to the new cell, which replaced the current one.


Description of ‘Bound’ function




Incorporate all the mass, accepted during Distribution, into the cells internal mass




Incorporate all the mass, accepted during Distribution, into the cells internal mass. If still no mass inside the cell, perform LàE morph. If the mass is too low (as indicated by vapor_factor) or there is at least one E neighbor, perform LàB morph.


Incorporate all the accepted particles (by moving them to internal list). If the list is empty, perform BàE morph. Otherwise, set all adjacent field elements to the average of all the particles’ velocities. If the mass is high enough, as indicated by cond_lo_factor, and there are no E neighbors perform BàL morph. If the mass is too high, as indicated by cond_hi_factor, first morph all E neighbors (and transfer one particle to each neighbor, to make the morph stable) to B and only then perform BàL morph on itself.




We leave it to the reader, to integrate those tables over the entire environment and see that the global outcome of each function is as it should be, according to Chapter 1.

Other Cell functions include some simple revealers and the Output function (to save the cell’s information to PHF file). They are simple enough to speak for themselves – see source files for details.



(2.8) Dynamic Configuration Module:


This project requires a very powerful and flexible configuration module. The range of simulated phenomena and the requirement for a maximum possible generality creates a great deal of variable parameters, which must remain user-defined. A full list of such parameters can’t be determined exactly beforehand and can change constantly throughout the program’s development and even afterwards. Therefore, the configuration module must provide the programmer with the possibility of introducing new configuration parameters and removing obsolete ones with a minimum amount of effort, including that of loading the new parameter from configuration file and command line. Needless to say, the module must also be convenient to use by both the programmer and the user.



The module is implemented as a Configuration class – a simple structure, which specifies all of the desired parameters (defined publicly). It also contains a linked list of pointers to GenericItems. A GenericItem is an abstract class, which can provide all the information about a certain parameter – how to initialize, load and maintain it. Basically, there is exactly one item in the list for every parameter, but that’s not an obligation.

In principle, there has to be a separate subclass of GenericItem for every possible parameter type, but the module also contains a subclass template, which can support any type, as long as it can be converted to String and back. This includes integer, floating point, character, string, boolean and a list of strings. See 2.11 on how to add new parameters (and their ConfigItems) and new types of parameters (e.g. Vector and Curve).

The template implements a number of abstract functions, declared in GenericItem and a constructor:




Locates and loads a new value from a given list of strings (see below)


Outputs the current value of the parameter to the specified output stream


Prints a predefined help string for that parameter


Resets parameter’s value and internal variables


To implement member functions, a ConfigItem must receive a reference to the corresponding parameter, a default value, com-string and help string (in this order) during construction. PrintValue and PrintHelp simply output the reference and help string respectively.

‘Redefine’ is the core of the loading procedure. It receives a reference to a list of strings and searches for a string among those which equals com-string. The next string in the list (after com-string) should contain a new value for the corresponding parameter. That’s why it’s important that String can convert itself to any relevant type. After that, Redefine removes both strings from the list, so they wouldn’t be traversed again.

All defined ConfigItems instances are inserted into the list.

The Configuration’s Redefine function simply traverses this list, calling Redefine of every item with a given list of strings, loaded from a file or from a command line (see 3.2 for details). In the end, this list must be empty, otherwise there is a format error (some strings weren’t identified by any item).

Likewise, Print and PrintHelp traverse items list and call PrintValue and PrintHelp of every item.

There is also a special kind of item – ConfigComment. It’s not connected to any variable and simply removes all the comments (inside /* */ pair) from the list of strings.



With this module, any parameter ‘x’ can be accessed from any point in the program using the syntax ‘config.x’ – doesn’t get easier than that. And adding/removing a parameter requires adding/removing exactly 2 simple lines of code. See 2.11 for details.

‘Validate’ function defines some desirable relations between configuration parameters. Definition is via VALIDATE macro, which works exactly like _ASSERT macro, except that it exists in Release mode as well.



(2.9) Polygonal to volumetric conversion (Volumetric module)


Before the program can start working, it must create a geometric environment, based on the geometric input files. The PHF files (as described in 3.3) are very simple to load, because they are the exact textual representation of the computational grid and their parsing is trivial.

However, building a PHF file for a complex environment can be tedious and this will have to be done all over again if Dx, Dy, Dz or K configuration parameters are changed. Therefore, a more convenient and flexible IRIT files should also be used. But before they can be used, they must be converted to a volumetric form (which is highly unusual for IRIT) to accommodate the computational grid. For example, if an IRIT file defines a solid sphere as a polygonal mesh, it must be approximated to its volumetric (voxelized) form and incorporated into the grid. The Volumetric module does exactly that.



Volumetric is derived from Array3D of type VoxelStatus, which is a set of scanning flags and a pointer to a surface normal (see below). The array has exactly the same sizes as the computational grid (although it occupies much less memory) and each array item corresponds to a computational cell.

The member function LoadPolygonalShape receives an IRIT object (of type IPObjectStruct), which must represent a closed (or nearly closed) surface in a form of a polygonal mesh. In the current implementation, the surface should be properly scaled and translated, to be inside the environment (for example, 0 < x < dr∙Dx). The parts of the surface beyond the scope of the environment will be cut off. The function also receives scan flags, which extend the conversion’s flexibility (see below).

When the function returns, the 3-D array contains the voxelized form of the surface and/or volume inside it, with each array item representing a single voxel.

The conversion algorithm is as follows:

1)      Convert all the polygons to triangles (this function is readily available in IRIT) and calculate a normal for each one (unless already provided).

2)      Create a voxelized approximation of each triangle, by superimposing a dense enough net of points on the triangle (as shown in figure 6) and “filling” all the voxels (by raising a ‘border’ flag) which contain at least 1 such point (determining to which voxel belongs a point is very easy). This will create a volumetric representation of the surface. Each voxel keeps an average normal of all the triangles that intersect with that voxel.

3)      If the surface is all we need, then the algorithm is complete. However, most of the time we’d like to fill the entire volume inside the surface. This is done by a 3-D variation of the scan-line algorithm, as shown in figure 7. The scan fills all the voxels it encounters after a ‘border’ voxel with a normal pointing at the scan direction (meaning, this is the “incoming” border), until the opposite border is reached. For a completely closed surface this will be enough, but small imperfections in the original mesh may cause the scan to “leak” through the surface (see figure 7c). This can be fixed in most cases, by performing another scan from a different direction and using both scans to determine the most probable shape of the volume. If that still doesn’t give the right results, another scan can be performed and so on. In a 3-D Cartesian system, there are 6 possible scan directions. The user can decide which scans to perform, by setting a proper flag in function’s parameters and viewing the results to see what scan combination provides the best results (64 possible combinations altogether). Of course, there is no need to try all 64 combinations – the user can guess the proper one relatively easy or he can activate all flags from the start, for the maximum power. Another idea is to add automatic combination finder. A good combination is the weakest one, which doesn’t create leaks (too many scans require more time and may cut more than necessary). Leaks can be detected even on the user level, so this extension can be done without altering the underlying algorithm.

The scan flags are as follows:

XN – scan along X-axis, in the negative direction.

YN – scan along Y-axis, in the negative direction.

ZN – scan along Z-axis, in the negative direction.

XP – scan along X-axis, in the positive direction.

YP – scan along Y-axis, in the positive direction.

ZP – scan along Z-axis, in the positive direction.

Each scan has its own flag (with the same name), which is set for the voxels, filled during that scan, so each voxel contains 8 flags (including the border and full). The ‘full’ flag represents the final outcome of all scans. For example, if all 6 scans have been performed, then:

full = border or [(XP or XN) and (YP or YN) and (ZP or ZN)]

This way, scans from opposite directions complement each other (in case one of them misses something) and scans from perpendicular directions leave only the intersection, thereby “closing the leak”. If only some of the scans have been performed, then only the relevant flags would be considered.

Notes: Voxels contain only the pointers to normals to conserve space, as most of the time there won’t be any normal defined for a voxel.






(2.10) Computation and output (Environment module)

The Environment class provides the main functionality for the entire project: input, computation and output. The following sections describe each function.


a) Constructor (used only once):

  • Creates default geometry – a solid box (Dx by Dy by Dz cells) filled only with air. The velocity field is set to 0 everywhere and pressure field equals air pressure in every cell.
  • Loads geometry input files (if any), using the Volumetric module to convert polygonal meshes into voxelized volumes (whereas PHF files define voxels directly). Each grid cell, which corresponds to a filled voxel morphs into a cell of proper type (e.g. solid, liquid, input, output). PHF defines cell type for every voxel, while IRIT defines a single type per object. All the files are loaded one by one, creating a union of all the files’ geometries.
  • Replaces all L-cells, which are adjacent to E-cells, with B-cells, resetting the Envelope Condition.
  • Calls Output function to store the initial state of the environment (for t = 0).
  • Runs a PHF editor or a loop of Evolve function, whichever is requested by the user. The evolution loops run until the environmental time reaches some user-defined limit.
  • When the constructor returns, the program ends.


b) Evolve:

The evolution step closely follows the design specifications (see Chapter 1), where each loop on the grid is performed first along the x, then along y and then along z axes. With evolution phases functions of the Cell, as defined in 2.7, the process is very simple:

  • Set phase to Calibration and calculate dt (time step), as specified in 1.4. This requires a single loop on every cell, during which the Calibrate function is also used to swap new field values with old ones.
  • Set phase to Calculation and call Calculate function of every cell.
  • Set phase to Balancing and call Balance function of every cell. Repeat the process until all calls return true (which means that the entire grid is balanced).
  • Set phase to Distribution and call Distribute function of every cell.
  • Set phase to Bounding and call Bound function of every cell. Repeat the process until the are no morphs detected (a morph is detected when Bound returns a different pointer than its host cell)
  • Advance environmental timer by dt and call Output if necessary (if the time has come).

Notes: by setting a proper phase (a static variable of FieldElement), the programmer can use GetValue and SetValue functions to access/modify the velocity field, without any second thought. The Field classes will make sure that the field behaves properly.


c) Output:

            The output function can save the current geometry in 2 forms (user can request each or even both):

  • The entire current geometry in a PHF file. Its format is described in 3.3 and the output process is very simple – to call Output function once for every cell.
  • The current liquid surface as an IRIT polygonal mesh. This one relies on the IsLiquidAt function, which returns true for any point in the environment containing liquid (false otherwise). The process is also pretty simple and it’s called the silhouette method: traverse all the faces between cells (including B’s sub-cells) and if one of them is full of liquid and the other isn’t, add a square representing this face to the mesh. It’s easy to see that in the end, the resulting mesh will contain a polygonal surface (of big and small squares), perfectly enclosing any volume of liquid (creating its 3-D silhouette). By setting the representation resolution fine enough, the surface will appear smooth. Of course one can apply a corner cutting or some another scheme to make it even more smooth.


d) EditPHF – a simple editor command line program, which creates a new PHF file for later input.


e) Destructor – is supposed to destroy the environment and de-allocate all the memory, occupied by the grid, but

since the destructor is called only once, at the very end, there is no point to it. The operating system will clean up all of the program’s context much faster then destroying cells one by one.



(2.11) Implementation level extensions.

1) Adding new configuration parameters:

  • To add a parameter ‘x’ of any type ‘T’, add “T x;” line to the Configuration class declaration (Config.h) and “RegisterItem(x, def_val, com_str, help_str);” to the constructor (Config.cpp). Optionally, define valid values for ‘x’ in the Validate function (Config.cpp).

Example:          Index Dx;                                                                     //  In Config.h

                              RegisterItem(Dx,10,"-Dx","x-size of the mesh");            //  In Config.cpp (c-tor)

                                    VALIDATE(Dx >= 3);                                                 //  In Config.cpp (‘Validate’ function)

  • To remove an item, reverse the preceding process.
  • The item can be of any type T, as long as String::String(T) and String::ConvertTo(T) are defined for it. If not, then these definitions must be added to String class (mstring.h)


2) To upgrade Curve into a true function in time, one has to re-implement all 5 functions in Curve.cpp. Of course, one may have to add more private members to Curve’s declaration (Curve.h), but the interface must remain the same. All other files remain unchanged.


3) The sub-cell array inside B-cells is somewhat redundant – it’s used only by IsLiquidAt and only during initial environment construction or during output. Creating these arrays only at these moments and destroying right after may reduce computation time, during which the heavy B-cells are constantly created and destroyed, hampering performance. Of course this may complicate the otherwise beautifully simplistic Evolve function.


4) Apply Matrix to IRIT objects before volumetric conversion to scale and translate it properly, so it would be entirely within the environment.




Chapter 3: User Guide


(3.1) Program requirements:

  • The program exerts a heavy load on the hardware, which is proportional to the user-defined spatial and temporal resolution. A practical environment can easily contain tens, even hundreds of thousands of cells. A grid of N cells usually occupies no less than 120N bytes (with about the same additional size for the particles) and must not exceed the size of an operating memory for the program to work. Output usually requires at least several megabytes per frame in IRIT format (sometimes up to 100MB per frame) and at least 2N bytes per frame in PHF format. This requires a lot of available disk space. Displaying the results efficiently, normally requires a graphics card with OpenGL.
  • The program has only been tested in Windows NT and XP environments.



(3.2) Configuration file and command line:

There are 2 ways to change the default configuration parameters: via configuration file or command line. The command line has higher precedence – it can override the configuration file and default values.

The format for the file and command line is exactly the same – a set of string pairs: <switch> <value>, separated by white spaces (blanks, ends-of-line, etc.).

  • <Switch> defines the parameter to which the new <value> will be assigned.
  • Boolean values are: true, false, yes, no, y, n, on, off. Can be either in upper or lower case.
  • Lists of string values are: {S1,S2,S3}, where Si is the i-th string in the list. No white spaces are allowed inside the string list.

The list of strings is used to specify a list of geometric input files. The files are loaded one by one in the specified order and the resulting geometry is the union of all the files, with the intersection’s type determined by the last relevant file.

  • Numeric types’ values are defined in the normal fashion.
  • There must be at least one white space between switch and value strings.
  • The amount of white spaces (including ends-of-lines) between parameters and between switch and value is unlimited.
  • Each switch string must appear no more than once in every configuration override (i.e. it can appear no more than twice – once in the file and once in the command line). This is done for the security reasons only – not to allow accidental multiple redefinitions.
  • Configuration file must be named “config.txt”.


The full list of switches is presented in the following list (you can get it by running the program with  /HELP ON):


/* <text> */                                         -- comments in the configuration file

-Dx <value>                                       -- x-size of the mesh (in unit cells)

-Dy <value>                                       -- y-size of the mesh (in unit cells)

-Dz <value>                                       -- z-size of the mesh (in unit cells)

-dr <value>                                         -- meters per cell

-K <value>                                         -- number of subcells per cell

-Nx <value>                                       -- number of generated particles in x dir.

-Ny <value>                                       -- number of generated particles in y dir.

-Nz <value>                                        -- number of generated particles in z dir.

-max_dt <value>                                -- optimal evolution step time [sec]

-min_dt <value>                                 -- minimal evolution step time [sec]

-finish_time <value>                          -- total simulation length [sec]

-frame_time <value>                        -- time between frame outputs [sec]

-min_dr <value>                                 -- minimum finite distance [m]

-interp8 <value>                                 -- perform 8-point field interpolation

-Gx <value>                                       -- external force x-acceleration component [m/sec^2]

-Gy <value>                                       -- external force y-acceleration component [m/sec^2]

-Gz <value>                                        -- external force z-acceleration component [m/sec^2]

-p_factor <value>                               -- dependence of pressure on density [(m/sec)^2]

-div_factor <value>                            -- determines allowed cell divergence during balancing (0..1)

-vapor_factor <value>                       -- min. filling of enclosed L-cell (forced vaporization) [0..1)

-cond_lo_factor <value>                    -- max. filling of enclosed B-cell (normal condensation) []

-cond_hi_factor <value>                    -- max. filling of B adjacent to E (forced condensation) []

-balancing_rounds <value>               -- max. number of balancing rounds (even if not converged)

-P0 <value>                                        -- air pressure [Pa]

-ro0 <value>                                       -- liquid normal density [kg/m^3]

-viscosity <value>                              -- liquid viscosity [m^2/sec]

-friction_factor <value>                     -- no slip coefficient of solids [0..1]

-cohesion_acc <value>                      -- L-p attractive force acceleration [m/sec^2]

-adhesion_acc <value>                      -- S-p attractive force acceleration [m/sec^2]

-progress_report <value>                  -- basic progress report

-extra_progress_report <value>       -- extended progress report

-statistics <value>                              -- some statistical analysis after each evolution step

-input_file_list <value>                      -- input files: *.itd or *.phf

-output_name_mask <value>            -- common part of output file names

-binary_output <value>                     -- output to a binary (as opposed to text) IRIT frame files

-PHF_output <value>                        -- output frames to PHF files

-IRIT_output <value>                       -- output frames to IRIT files

-binary_output <value>                     -- output to a binary (as opposed to text) IRIT frame files

-draw_sinks <value>                          -- draw Output cells as filled with liquid

-draw_ground <value>                       -- draw ground rectangle (for visual reference)

-XNScan <value>                               -- perform x-wise linescan with neg. step

-YNScan <value>                               -- perform y-wise linescan with neg. step

-ZNScan <value>                               -- perform z-wise linescan with neg. step

-XPScan <value>                               -- perform x-wise linescan with pos. step

-YPScan <value>                                -- perform y-wise linescan with pos. step

-ZPScan <value>                                -- perform z-wise linescan with pos. step

/HELP <value>                                  -- print this help screen and exit

/EDITOR <value>                             -- activate PHF-file editor module



(3.3) Geometric input files formats:

Being designed specifically for this program, the PHF (PHotoFlow) files are more convenient and compact than IRIT files. However, they must also be created specifically for every simulation environment.

The recommended way is via integrated PHF editor. To run the editor, type “photoflow /EDITOR ON” at the command line and press Enter. Alternatively, you can add “/EDITOR ON” to the configuration file.

This will launch a simple built-in editor, with which you can edit the grid, cell by cell. Grid’s sizes are defined by the current configuration as specified above. After you launch the editor, just follow the on-screen instructions. When you quit the editor, the program terminates, creating <name>0.phf file, containing the resulting grid (<name> is defined by the “output_name_mask” parameter). That file can be renamed (keeping the extension of course) and used as the input file for any environment with the same Dx, Dy, Dz and K. The “PHF_output” option must be turned on to run the editor.


Another way is to edit the PHF file directly, using any text file editor (e.g. Notepad). Although not recommended for beginners, this is sometimes a much faster way to insert small modifications into the input geometry. The PHF format is as follows:


<Dx> <Dy> <Dz> <K> (<type> [<parameters>])*


where <Di> and <K> are numeric values of the proper configuration parameters, followed by cell types and their parameters at every grid point. The order of the cells corresponds to the tri-linear traversal of the grid, first along x-axis, then along y-axis and finally, along z-axis. That is, the first cell in the list is the (0,0,0) cell, than (0,0,1), (0,0,2), … , (0,0,Dx-1), (0,1,0), … … … , and finally: (Dx-1,Dy-1,Dz-1).

The following table describes the format for each cell type:








No parameters



Vx Vy Vz

Input direction – a vector of 3 floating point numbers



Vx Vy Vz

Output direction – a vector of 3 floating point numbers




No parameters




K3 zeros or ones, describing array of sub-cells (empty | filled)




No parameters


All parameters appear in textual format and can be separated from each other by an arbitrary amount of white spaces. It is convenient to write cells in rectangular blocks, representing horizontal grid levels.

Example (7x5x3 grid, surrounded by solid walls, with single input cell in the center, ejecting liquid in (1,0,0) direction, with K = 10):


7 5 3 10


6 6 6 6 6 6 6

6 6 6 6 6 6 6

6 6 6 6 6 6 6

6 6 6 6 6 6 6

6 6 6 6 6 6 6


6 6 6 6          6 6 6

6 1 1 1          1 1 6

6 1 1 5 1 0 0 1 1 6

6 1 1 1          1 1 6

6 6 6 6 6 6 6


6 6 6 6 6 6 6

6 6 6 6 6 6 6

6 6 6 6 6 6 6

6 6 6 6 6 6 6

6 6 6 6 6 6 6



As for the IRIT files, they have to comply to certain conditions:

Only IPObjects of polygonal mesh type are accepted, everything else is ignored.

IPObjects can define different types of geometry, using “CellType” integer object parameter. The numeric values for each type are provided in the cell types table (see 2.7). If CellType parameter isn’t provided, 6 (Solid) is assumed by default.

Objects representing valves (I) and sinks (O) must provide 3 Curve type parameters “Vx”, “Vy” and “Vz”.



(3.4) Running the program:

To run the program, simply activate the executable. For example, type “photoflow” at the command line and press Enter. No other files are necessary, except input and configuration files, if requested by the user.

The valid input files are of *.itd, *.ibd and *.phf types and a single config.txt configuration file.

If the “-progress_report” is ON, then the program will print a notification before each evolution step and before each phase of that step. Some messages are also printed during initial environment construction.

If the “-extra_progress_report” is ON, the program will print some additional notifications during each phase.

If the “-statistics” is ON, then the program will print some grid statistics (current number of cell of each type, total volume of liquid, etc.) after each evolution step.



(3.5) Displaying results:

Currently, there is no graphical viewer for PHF files, but when viewed with a text editor, the geometry is often evident. Also, the program can serve as a converter from PHF to IRIT files. For example, the following command line will convert in.phf into out0.itd:

Photoflow   –input_file_list {in.phf}    –output_name_mask  out   –IRIT_output  YES


IRIT files can be viewed with one of the IRIT viewers: wntdrvs.exe or wntgdrvs.exe.



(3.6) Error messages:




What to do


Out of memory

Close other programs, try to decrease grid size and/or resolution (larger dr, smaller K).


Unknown or redundant parameter set: <S>

Switch is either undefined or appears more than once in a configuration file or command line. Fix or remove the switch and its value.


<S> has missing or improper argument after it

Each switch string must be followed with a value of a proper type. A comment must be enclosed with /* */ which are separated from the comment string(s) by at least one white space.


Unable to open <S>

An input file must be located in the same directory as the program, or a full path to the file must be provided with its name.


Error loading file <S>

Usually occurs while trying to load a PHF file with unmatched Dx, Dy, Dz or K values. Adjust these values or use a different PHF. Also occurs at any anomaly in the input file format. Also occurs at any disk I/O error.


IPObject attribute wasn't provided

IRIT files representing Input and Output geometry must specify their parameters, as described in 3.3.


Invalid parameter values. Must have: <S>

Validation rule is violated. See 3.7.



(3.7) Choosing optimal configuration parameters:

First of all, the configuration parameters must obey certain validation rules, as follows:

  • Dx >= 3
  • Dy >= 3
  • Dz >= 3
  • min_dr > 0.0
  • dr > min_dr*2
  • K >= 3
  • min(Nx,Ny,Nz) >= K
  • min_dt > 0.0
  • max_dt > min_dt
  • frame_time >= max_dt
  • finish_time > frame_time
  • vapor_factor >= 0.0 && vapor_factor < 1.0
  • cond_lo_factor > vapor_factor
  • cond_hi_factor > cond_lo_factor
  • P0 > 0.0
  • ro0 > 0.0
  • viscosity >= 0.0
  • viscosity <= 0.5*sqr(dr)/max_dt
  • friction_factor >= 0.0 && friction_factor <= 1.0
  • adhesion_acc >= 0.0
  • cohesion_acc >= 0.0
  • div_factor > 0.0 && div_factor < 1.0
  • progress_report || !extra_progress_report
  • PHF_output || IRIT_output
  • balancing_rounds >= 0


These rules are vital for the computation results to make sense, but there are other rules, which cannot be so rigidly defined:

  • Note that each computation step adds a cumulative error, which can be as high as O(max_dt3 + max_dt2∙dr2). For optimal performance, it’s recommended to choose such max_dt and dr, so that dr2 ~ max_dt. In this case, the error will be up to O(max_dt3) per step. After n time steps, the error may increase to O(sqrt(n)∙max_dt3). Therefore, max_dt should take the finish_time into account. The larger is the latter, the smaller should be the former to minimize the cumulative error.
  • On the other hand, max_dt can’t be too small, because it would increase the computation time, but that’s up to the user. Also, the valves have a certain threshold of max_dt, below which they can no longer work. This threshold depends on particle density, so increasing Nx, Ny or Nz can solve the problem even for small max_dt, although it also increases memory consumption. It’s up to the user to find a proper compromise between all those considerations.
  • The pure physical parameters (P0, ro, viscosity, Gx, Gy, Gz) are the easiest to choose, because they are known for every realistic environment. Density and viscosity values can be found in various scientific tables. The default values are given for the water.
  • The recommended value for div_factor is 1/max(Nx,Ny,Nz) but that can be changed, if a more (or less) strict incompressibility is required.
  • Adhesion_acc and cohesion_acc should be equal, unless you desire to simulate effects, where surface tension is most evident, like a small liquid volume assembling itself into a spherical droplet. In this case, cohesion_acc (which causes the surface tension) should be larger than adhesion_acc. Normally, both should be the same order of magnitude as the gravity acceleration Gz.
  • vapor_factor, cond_lo_factor and cond_hi_factor can be used to calibrate the incompressibility condition. They should be very close to 1 (while still obeying the validation rules) for best results, but sometimes the values far from 1 can also be useful.



Chapter 4 – Downloads

(4.1) Testing examples:

All the following tests were conducted on a common variety Pentium IV, 2.8 GHz with 512 MB of RAM.

The resulting images are unprocessed (simple processing would have made the images even more realistic) and only basically rendered with Wntgdrvs viewer. The frames were created using AT Screen Thief 3.6 and the movies were assembled on Windows Movie Maker 5.1. The movies run at a slower rate than the actual one, to allow the viewer to see every aspect of the flow more clearly.

Note that only the liquid is shown (including the one inside valves and sinks) – the solid obstacles and the valves themselves are invisible.


a) Tub:

Description: The initial environment of 1,331 cells (11x11x11) is a cubical box (tub), partially filled with standing water up to 3 cells deep. A valve near the top of the box (above the water) starts ejecting water in horizontal direction, at 0.5 m/sec. The system is evolved from t = 0 to t = 1 second. The visual resolution is 66x66x66. The movie runs 6 times slower than the real flow rate.

tub.wmv (450 KB)



Comments: Notice the shape of the water jet (until it hits the front wall of the tub). Despite the horrible computational resolution, it clearly gets the proper parabolic shape. If you look at the front water line, you can see a wave spreading from the center, reflecting from the walls and performing a wall to wall oscillation as a standing wave – a perfect physical correctness (even if hardly visible at this resolution). If you look closely, you can see that the initial depth of the tub decreases a little at the beginning. This is an expected numerical error, caused by the extremely coarse resolution, but this error is contained within a depth of one cell and, obviously, no numerical solution can guarantee correctness within computational cell. Actually, with the proper initial geometry this effect can be fixed, but in this case it’s beneficial, because it compresses the outer layer of the water, which facilitates wave generation. Without this compression, the water would have been less turbulent.

Performance: The simulation ran for 6 minutes, performing 150 steps and generating 26 frames. Maximum memory usage – 22 MB, required disk space – 78 MB.

Input files:


tub.phf (3 KB)

config.txt (1 KB)




b) Drop:

Description: The initial environment of 15,625 cells (25x25x25) is an empty cubical box, with a water blob hanging in mid air, near the top. The system is evolved from t = 0 to t = 0.4 seconds. The visual resolution is 125x125x125. The movie runs 12.5 times slower than the real rate.

drop.wmv (200 KB)



Comments: The cohesion forces (attraction between particles) are evident in this example. The resulting surface tension changes the initial unnaturally cubical shape of the blob into a more spherical one (although there isn’t enough time before the blob hits the floor). The process of shape adjustment creates expected, yet funny looking oscillations of the blob, which are also responsible for the asymmetry of the final splat. The 8 cubes in the corners of the box are added for proper scaling (these are just inert valves) – without them the scaling would change, creating an interesting, but distracting illusion of camera motion, as the blob falls down. The physical correctness is most evident in the fact that it takes between 0.18 to 0.19 seconds for the blob to reach the bottom, which is consistent with the theoretical value. This example shows how small the overhead of the non-liquid environment is. Even though the overall resolution is about 40 times larger, than in the previous example, the simulation runs almost 5 times faster!

Performance: The simulation ran for 32 seconds, performing 71 steps and generating 41 frames. Maximum memory usage – 10 MB, required disk space – 23 MB.

Input files:


drop.phf (32 KB)

config.txt (1 KB)




c) Waterfall:

Description: The environment of 23,750 cells (50x19x25) contains 13 elevated valves, generating a river, which flows at 1 m/sec off a cliff into a riverbed below, with raised sides (like a large tube). On the other side, the water disappears into a large (visible) sinkhole (the sinkhole has the same shape as the riverbed’s profile). The system is evolved from t = 0 to t = 2 seconds. The visual resolution is 250x95x125. The movie runs 3 times slower than the real flow rate.

waterfall.wmv (810 KB)



Comments: The surface has a high degree of friction, which causes the front portion of the flow to slow down and creates a rising wave. The cohesion forces are deliberately turned off, to increase the splashing effects. This was the heaviest test ever performed on the program to date.

Performance: The simulation ran for 28 minutes, performing 560 steps and generating 51 frames. Maximum memory usage – 130 MB, required disk space – 756 MB.

Input files:


waterfall.phf (48 KB)

config.txt (1 KB)





(4.2) The program:

This zip archive contains the main executable – photoflow.exe and the default configuration file (optional) – config.txt: (371 KB)





(4.3) Source files:

This zip archive contains 41 C++ source files (*.h and *.cpp) of the main executable:

project. zip (40 KB)





(4.4) Links:

·        IRIT Manual –

·        Another project in flow simulation –

·        Ron Fedkiw homepage (with flow simulation video clips) –

·        Water simulation and rendering –

·        An attempt to break the myth surrounding the Coriolis Force:





Created: Jan. 26, 2005                         Last modified: Feb. 4, 2005                      Technion – Computer Science Department