[aspect-devel] Progress in writing the mantle convection code Aspect

[Wolfgang Bangerth]

On 10/03/2013 04:58 PM, Eric Heien wrote:
> I'll take an initial stab at a file format for reading particle data (as
> opposed to an initializing function).  This will also be useful for
> particle checkpointing, otherwise we would have to store lots of particle
> data in the checkpoint file.

I think it would actually be simplest if we just put it into the checkpoint 
file. It needs to be stored somewhere anyway, and the checkpoint file already 
has a well-defined file format that is easy to use. Why not use it?


>  My plan is to support ASCII and HDF5 formats
> in the same way as is currently done - columns of ID, X, Y, Z, scalar0,
> scalar1, etc in a single file.  For now I'll just allow up to 3 scalars,
> more than that will generate an error.  Can you think of any potential
> problems with this?

You could define the number of scalars associated with each particle within 
the file somewhere at the top (or with each particle, if that were easier). I 
have no problem with it if we just error out if there are more than 3 
properties, but I don't think this should be hard-coded in the file format.


>> Maybe some pseudocode will help (think of a single processor, it's not
>> much more complicated on multiple processors):
>>
>> std::vector<double> relative_cell_probability (n_active_cells); int index
>> = 0; for (cell=....; ++cell, ++index) relative_cell_probability(index) =
>> cell->volume() * probability_density(cell->center());
>>
>> // normalize all probabilities so that the sum == 1
>> relative_cell_probability /= relative_cell_probability.l1_norm();
>>
>>
>> // compute a cumulative probability. the last value of // cumulative_c_p
>> will be one after this std::vector<double> cumulative_cell_probability
>> (n_active_cells); cumulative_cell_probability[0] =
>> relative_cell_probability[0]; for (i=1...) cumulative_cell_probability[i]
>> = cumulative_cell_probability[i-1] + relative_cell_probability[i];
>>
>> // draw P particles for (p=0...P) { // find which cell to put a particle
>> in double r = rand();   // between 0 and 1 int index = first element of
>> cumulative_c_p so that cumulative_c_p[i] > r;
>>
>> // place a particle in cell 'index' }
>>
>> Makes sense?
>
> Yes, this is similar to what I do in my current setup.  I was thinking the
> original proposal was to assign each cell a fixed number of particles based
> on the probability density.  Do we want to ensure the probability density
> function covers the domain exactly?  It shouldn't be a big problem because
> of the renormalization but the results can be unexpected if the
> distribution doesn't match the domain.

I was thinking the probability density is simply a function f(\vec x) defined 
on all of \R^dim. If it wants to return zero in some parts of space, then 
that's no problem.


> Also, this scheme (and the one I already implemented) is still not
> completely random for multiple processors.  This is because each processor
> will need to generate an exact number of particles within the local
> subdomain, but this isn't a huge problem.

How do you do it right now?


 > Another thing to be careful of
> is to ensure each processor has a different random seed, otherwise you can
> get identical particle distributions  on different processors if they have
> congruent geometry (as can commonly happen with structured meshes and 2^N
> core counts).

That's an interesting problem. Good point.


> One thing I'm confused about is the "place a particle in cell" line - if
> the distribution is relatively smooth over a cell then approximating it
> with a uniform distribution is fine.  But if it changes significantly  I
> don't see how to place particles within a cell without requiring the
> cumulative distribution function.  Any thoughts about this?

It's a valid point, in particular if you want to restrict particles to a 
subdomain (i.e., your probability function is zero/one). But I think that if 
your mesh is not fine enough to adequately represent this, then you're going 
to be out of luck on other levels as well.

So my idea was to just evaluate the probability density at the cell center and 
assume it to be constant within each cell. We could, if we wanted to, improve 
this doing a quadrature over the cell, but let's not worry too much about this 
right now.

Best
  W.

-- 
------------------------------------------------------------------------
Wolfgang Bangerth               email:            bangerth at math.tamu.edu
                                 www: http://www.math.tamu.edu/~bangerth/