Fortran - Representing structures for simulations

Author

Representing structures for simulations

justabeginner

2005-04-24, 8:57 am

Hello. I am planning on programming my own small smooth particle
hydrodynamics code. I have no prior experience with programming such
code, or even other meshed methods. This SPH code will be interfaced
with another program I have written to study hypervelocity impact.

I guess the main question I have at the moment is how I should
represent structures in my code. How is it normally approached, in both
meshless and mesh-free methods?

For example, say you have a 2D plate of material. I could discretise
space and represent the plate in a 2D memory array. The size of my
array would be proportional to the size of the plate, or the level of
discretisation. This method would be a memory hog for large structures
or fine resolutions, but it would be faster (assuming the structure can
be stored in RAM) because you immediately know which elements/particles
are adjacent to one another. You don't need to run adjacency checks to
determine which sections affect which - you can tell just by the array
co-ordinates.

Alternatively, I could store element/particle information in memory in
a manner unrelated to their spatial positions (like representing a 3D
structure in a very long 1D array). Spatial co-ordinates are stored as
entries in the array. This means that I will need to do adjacency
checks, which is CPU intensive.

Suppose you had to represent a cube of material (so no voids). The
amount of information you must store is the same, regardless of which
storage scheme you choose. Therefore the best scheme is to store it in
a 3D array, since you'd be using the same memory (theoretically) in
either case, but losing on CPU costs for the 1D array scheme. However,
if the material has lots of voids, then the best option may be the long
1D array, because it doesn't allocate memory to wastefully defining
voids.

Are there alternative means of storing structures in memory I have not
thought about? What thoughts do you guys have?

Thanks.

Jon Harrop

2005-04-24, 8:57 am

justabeginner wrote:
> Are there alternative means of storing structures in memory I have not
> thought about?

Yes, trees, graphs, queues, stacks...

--
Dr Jon D Harrop, Flying Frog Consultancy
http://www.ffconsultancy.com

justabeginner

2005-04-24, 3:57 pm

not[color=darkred]
>
> Yes, trees, graphs, queues, stacks...
>

Thanks, from what I've read - correct me if my understanding is wrong:

Lists contain data entries which points to where in memory the
subsequent entry is. IMHO these are OK for sequential data (like diary
entries, for example) but not OK for >1D data (like spatial
distribution of materials in a composite).

Stacks are data structures where the last entry to go in is the first
out.

Queues are (like in real life) data structures where the first entry is
the first out.

Trees order data in branching clusters, and the program follows one
cluster to another without looping back to previous clusters (still not
100% clear on this, but I think it would be clearer if I had to program
one myself).

Lattices are like trees, but allow looping back to prior data sets.

Graphs...I'm quoting: "If your data structure has loops so that you can
return to some data item by following a sequence of pointers from it,
then technically you have a directed graph, often just called a graph."

It still seems that none of them are particularly suited to
representing physical structures. For my problem, I will need to:
1) Record the states (velocity, temperature, phase, etc.) of different
particles;
2) Be able to quickly run proximity searches to locate neighbouring
particles.

It still seems like the best option is to record a N-dimensional
structure in an N-dimensional array, with the N+1th dimension storing
the state variables.

Gordon Sande

2005-04-24, 3:57 pm

justabeginner wrote:
>
> not
>
>
>
<snip>
>
> It still seems that none of them are particularly suited to
> representing physical structures. For my problem, I will need to:
> 1) Record the states (velocity, temperature, phase, etc.) of different
> particles;
> 2) Be able to quickly run proximity searches to locate neighbouring
> particles.
>
> It still seems like the best option is to record a N-dimensional
> structure in an N-dimensional array, with the N+1th dimension storing
> the state variables.
>

A fairly natural setup would be to have a table of attributes against
particle ID. That is 2d with one for ID and one for attribute number.
The practicalities of mixing data type confuse this considerably so
there may be several such tables.

To find the particle by position one would use a kd-tree, or one of
the related data structures, to allow searches for neighbouring
particles which are known by ID. The kd-tree is effectively an
index structure for a 2d table of ID vrs positions.

Notice that you stated what operations were of interest rather than
just asking for some data structure that could represent some data.
Some of your comments could be taken as an indication of a bit of
confusion between the F90 notions of rank and dimension. A 3d position
of ( x, y, z ) is a rank 1 array of dimension 3. The Fortran subscript
goes with the physical dimension. Same word in different context has
different meaning. (Take two Aspirins and get a good nights sleep!)

NuclearWizard

2005-04-24, 8:57 pm

I would check out some of the papers on SPH simulations to decide how
you should proceed. The Los Alamos National Lab web site has some info
including sweet simulation animations:
http://qso.lanl.gov/~clf/main.html
Although I am not familiar with SPH, but it seems like you have a
material with various space/time dependent properties as well as
particles which you are tracking through this mesh.

I would have real/complex arrays for the unknowns, a structure for the
geometry, structure for the fundamental data, and inquiry procedures so
that given particle with a certain direction, at a certain point in
space, with a certain energy you may determine the properties of the
media it is in and will be in until the next interaction is sampled.
Some setup procedure would form the system of equations (matrices,
tensors, whatever) using inquiries to the geometry and fundamental data
structures. Then a solver operates on the system of equations and
returns one or more of the unknowns. Then the next setup procedure and
solver combination is run to get some more unknowns and the simulation
continues until some convergence criteria for all unknowns is reached.

Jon Harrop

2005-04-24, 8:57 pm

justabeginner wrote:
> Lists contain data entries which points to where in memory the
> subsequent entry is. IMHO these are OK for sequential data (like diary
> entries, for example) but not OK for >1D data (like spatial
> distribution of materials in a composite).

Yes. Lists are ideal when you want to accumulate a sequence of elements
which is of unknown length.

> Stacks are data structures where the last entry to go in is the first
> out.
>
> Queues are (like in real life) data structures where the first entry is
> the first out.

Yes and yes.

> Trees order data in branching clusters, and the program follows one
> cluster to another without looping back to previous clusters

That is a fairly good description. Trees are recursive data structures which
contain 0 or more child trees.

In the interests of correctness and performance, trees are often restricted
to having either zero or two children, known as binary trees.

Balanced binary trees are a very important class of data structures which
are much faster than arrays for many useful operations (such as inserting
an element into an ordered container).

There is a huge volume of literature on trees as they are one of the most
basic and fundamental data structures.

> (still not
> 100% clear on this, but I think it would be clearer if I had to program
> one myself).

If this is for a big project then you should definitely study trees before
starting to code.

> Lattices are like trees, but allow looping back to prior data sets.

I would call that a graph.

> Graphs...I'm quoting: "If your data structure has loops so that you can
> return to some data item by following a sequence of pointers from it,
> then technically you have a directed graph, often just called a graph."

Not quite, a graph with "loops" is a cyclic graph. There are many other
types of graph.

A tree is a type of directed graph. Sometimes you want several either
distinct or partially overlapping trees (a kind of tree with 2 or more
roots), this is a directed acyclic graph (DAG). This data structure pops up
all over the place, e.g. inside compilers.

In general, a graph is a set of vertices and a set of edges joining pairs of
vertices. If the edges have a direction associated with them (i.e. edge "x"
goes from A to B but not from B to A) then it is a directed graph.

There is also a large volume of literature on graphs.

> It still seems that none of them are particularly suited to
> representing physical structures. For my problem, I will need to:
> 1) Record the states (velocity, temperature, phase, etc.) of different
> particles;

An array would be ideal for this.

> 2) Be able to quickly run proximity searches to locate neighbouring
> particles.

An array of particles is not likely to be ideal for this. You may wish to
use an array of lists representing the lists of particles in each region of
space (known as voxels in 3D, good if "neighbouring" means within a
constant cutoff radius), or a tree representing a recursive subdivision of
space into progressively smaller parts (good if you have long-range
interactions, e.g. Coulombic, where you want to group the effects of large
numbers of distant particles).

> It still seems like the best option is to record a N-dimensional
> structure in an N-dimensional array, with the N+1th dimension storing
> the state variables.

This is very unlikely to be correct. One could say that you are not seeing
the forest for the trees: you correctly identified several low-level
optimisations which could speed up your program by 10-20% but have failed
to consider any of the possible optimisations which could speed up your
program up by several orders of magnitude.

I strongly suggest you dig out some suitable books. Firstly, study
asymptotic algorithmic complexity, then study common data structures and
algorithms.

I recommend "Introduction to algorithms" by Cormen and Leiserson.

If you are willing to consider other languages, then I'd also recommend my
own book "Objective CAML for Scientists":

http://www.ffconsultancy.com/produc..._for_scientists

This may be a good idea as the most relevant data structures are very
difficult to implement in Fortran but trivial in many other languages.

--
Dr Jon D Harrop, Flying Frog Consultancy
http://www.ffconsultancy.com

Brooks Moses

2005-04-25, 8:57 pm

justabeginner wrote:
> Hello. I am planning on programming my own small smooth particle
> hydrodynamics code. I have no prior experience with programming such
> code, or even other meshed methods. This SPH code will be interfaced
> with another program I have written to study hypervelocity impact.
>
> I guess the main question I have at the moment is how I should
> represent structures in my code. How is it normally approached, in both
> meshless and mesh-free methods?

I presume you mean "meshed" and meshless methods -- meshless and
mesh-free are the same thing.

In meshed methods, there are two types of meshes: structured and
unstructured. (Note that this distinction is approximate; there are
meshes that are somewhat vaguely both.)

A structured mesh is, roughly speaking, one that can be stored in a
multidimensional array. This does not imply fixed resolution, or even a
cartesian grid; the grid can be distorted in various ways (e.g., a polar
grid) to adjust the resolution in various regions.

An unstructured mesh is one that's made up of individual grid elements
that can be randomly aligned; these are usually triangles or
quadrilaterals in 2D, or tetrahedrons or six-faced blocks in 3D. These
would typically be stored as a 1D array of grid cells, with the
adjacency information also stored. (Depending on arrangement, the
storage may be a 1D array of vertices with pointers to the cells that
contain them, and a 1D array of grid cells which contain pointers to the
relevant vertices). Because the adjacency information is precomputed
and stored, it is not especially CPU-intensive to handle.

For large complicated geometries where a pure structured mesh may not be
appropriate, one common approach is a "block" method, where the geometry
is broken up into subdomains, and each subdomain is discretized in a
structured mesh. Again, the adjacency information between subdomains is
stored.

Now, in mesh-free methods, things become trickier -- and simpler. In
the purest sense, because there is no mesh, there is no topological
structure, and therefore the particles _must_ be stored in a long 1D
array.

As can be learned from the unstructured meshed methods, however, long 1D
arrays are not CPU intensive if you have a good way of handling the
adjacency data. The problem is complicated, however, by the fact that
the particles are generally not fixed, and so the adjacency information
changes. Thus, some cleverness is called for. I can't claim to be
anywhere near an expert on the matter, but here are some things that
occur to me from having read a few papers on SPH and related subjects:

One simple method is storing, with each particle, a list of all of the
particles within X distance from it. At each timestep (or after each
suitable number of timesteps), you iterate through the particles, and
check all its neighbors and see if they're still within that distance.
Additionally, you check all the neighbors of its neighbors, and see if
they've become neighbors of it. And then you update the lists, and
continue. This, of course, requires that particles not be changing
positions too quickly, or you might not notice some new neighbors.
Also, it starts becoming very CPU intensive unless the neighbor lists
are small, or if particles interact with other particles that aren't on
their neighbor lists.

Another method is to overlay the method with a grid of some sort
(generally a structured mesh), and with each grid cell store a list of
particles that are in that cell. Then, you need only do adjacency
checks between particles in the grid cell and particles in that cell or
adjacent cells.

There are undoubtably many more methods, and means of refining these
methods so that they work orders of magnitude better than naive
approaches might; I second the recommendation of another poster to look
into papers in the field. In this day and age, with access to decent
online databases (if you're at a university, they should have access to
something like the "ISI Web of Science"; failing that, Google Scholar is
probably sufficient though much poorer), it shouldn't take more than a
day or two to find much of the relevant literature in the field.

Beyond all the program design questions, there's the programming-level
issue of storage. Unless you're going with a structured mesh (where the
adjacency calculations are trivial), you need to store two kinds of
information: particle adjacency information, and particle state
information. Those need different kinds of storage.

Particle state information will virtually always be stored in a 1D array
or set of 1D arrays, indexed by particle number. (There may be a
strategy for re-using particle numbers if a particle goes out of the
domain; I won't get into that. There may also be a strategy for how you
choose your particle numbers to make memory accesses work better, but
you absolutely should not be worrying about that at this stage.) I
would suggest defining a type called ParticleState or something like
that, containing named elements for each component of the state;
alternately, you could use a 2D array with the second dimension being an
index for the state components, but that's rather less clear.

Particle adjacency information is a different matter. This is where the
lists, trees, graphs, lattices, and/or other more-complicated things
that Jon Harrop mentioned are useful. Functionally, they probably get
implemented in Fortran as a lot of arrays containing integers that
"point" to other elements in the array or to elements in the particle
state array, though they may contain actual Fortran pointers to other
arrays, depending on how you implement your lists.

You may also find it useful to read (or at least skim) a book about data
structures, so that you get an idea of what sorts of things are
possible.

- Brooks

--
The "bmoses-nospam" address is valid; no unmunging needed.

Brooks Moses

2005-04-25, 8:57 pm

justabeginner wrote:
> It still seems that none of them are particularly suited to
> representing physical structures. For my problem, I will need to:
> 1) Record the states (velocity, temperature, phase, etc.) of different
> particles;
> 2) Be able to quickly run proximity searches to locate neighbouring
> particles.

The way to quickly run proximity searches is to have most of the work
already done, by storing a set of information such that you don't need
to search through all the particles. This set of data is something that
you'll want to put a fair bit of work into figuring out how to store,
because it doesn't look anything like a fixed-dimensional array.

> It still seems like the best option is to record a N-dimensional
> structure in an N-dimensional array, with the N+1th dimension storing
> the state variables.

This is called a structured mesh. I was under the impression you were
wanting to implement a meshless method. :)

Implementing things in an N-dimensional array works well, _only_ if the
points at which the states are stored are lined up in rows and columns
in a very structured manner. If the data points are not structured in
that way, then there is no obvious way to map them into the cells of the
N-dimensional array.

Also, as mentioned earlier, don't use the N+1th dimension for the state
variables. Define a data type that contains named components for each
state variable, and make an array of that -- it's much clearer to read.

- Brooks

--
The "bmoses-nospam" address is valid; no unmunging needed.

Greg Lindahl

2005-04-27, 3:58 am

In article <1114328588.607230.82890@g14g2000cwa.googlegroups.com>,
justabeginner <nom_de_plume79@yahoo.co.uk> wrote:

>I guess the main question I have at the moment is how I should
>represent structures in my code. How is it normally approached, in both
>meshless and mesh-free methods?

I would go look at othe SPH codes for examples. I think there are
several freely-available versions out there. I suspect that most
people here don't know what SPH is and won't give useful answers after
reading your terse description of how it works.

-- greg

Sponsored Links