dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Jan 13, 2009 at 03:31:33PM -0700, Bartosz Sawicki wrote:
> On 12/01/09 05:41 PM, Johan Jansson wrote:
> > Bartosz Sawicki wrote:
> > 
> >> On 09/01/09 09:38 AM, Johan Jansson wrote:
> >>> Hi!
> >>>
> >>> This looks interesting, could you please post references to the 
> >>> papers you mentioned about the algorithm?
> >>
> >> Basic concept is from Rivara, but I'd like to avoid recursive 
> >> algorithms, which has usually efficiency problems.
> >>
> > 
> > I don't think recursive algorithms have efficiency problems, what do you 
> > mean by that? Quicksort is a recursive algorithm for example, and it's 
> > one of the fastest sorting algorithms.
> 
> You are right, recursive can be very fast. I've found idea of iterative 
> Rivara, which seems easy to implement and robust - so I implemented it.
> I'm not convinced that iterative is better than recursive, or opposite.
> 
> >> I've found in some papers [1], where they use similar algorithm but 
> >> with  iteration procedure. Propably something like this is also 
> >> implemented in FEtk.org (GAMer package). I write "probably" because 
> >> the code is still unpublished, and I can't find any documentation 
> >> about it.
> >>
> >> [1] N. A. Baker, D. Sept, M. J. Holst, J. A. McCammon: The adaptive 
> >> multilevel finite element solution of the Poisson–Boltzmann equation 
> >> on massively parallel computers,IBM J. RES. & DEV. VOL. 45 NO. 3/4 
> >> MAY/JULY 2001
> > 
> > Thanks for the reference, I'll look this paper up.
> 
> On this paper they only mention that they use iterative bisection 
> algorithm. I could find anything more about iterative version of Rivara.
> 
> >> Algorithm with I've implemented. It wasn't directly taken from any paper:
> >>
> >>   while( there are any marked cells ):
> >>      forech cell which is marked:
> >>         find the longest edge
> >>         if the edge hasn't been bisected:
> >>            bisect longest edge e
> >>            remember bisected edge and new vertex
> >>         else :
> >>            use vertex which was remembered
> >>         bisect cell
> >>      foreach bisected edge:
> >>         mark all cells which incidence the edge
> >>
> >> The weak point of the algorithm is the storing of bisected edges, and 
> >> serching in them. I use double index STL mapping for that purpose - 
> >> which works pretty well, but that's the slowest part of the algorithm.
> >>
> >> I have some ideas, which can make it faster, but it hasn't been tested 
> >> yet. Stay tuned ...
> >>
> >> The mesh quality is fine, because _only the longest edge_ is bisected.
> > 
> > 
> > Ok, good! I see that the algorithm you've implemented is indeed an 
> > iterative variant of the Rivara algorithm. I think this could be very 
> > useful to know. But as you mention your implementation has some 
> > performance issues, could it also be because of the mesh 
> > re-initialization (connectivity) in each iteration?
> 
> You are deeply right. mesh.init(1) takes a lot of time - I think that it 
> is possible to completely eliminate replace edge MeshFunction with map. 
> It should definitely speed up the algorithm.
> 
> 
> > Here is a bundle of the Rivara refinement:
> > 
> > http://www.csc.kth.se/~jjan/dolfin-transfer/rivara.bundle
> > 
> > Which exists as:
> > 
> > dolfin/mesh/RivaraRefinement.cpp
> > 
> > together with a performance test in:
> > 
> > sandbox/refinement-perf
> > 
> > You can run the test as:
> > 
> > ./demo rivara
> > 
> > for the recursive Rivara refinement and:
> > 
> > ./demo iterative-Rivara
> > 
> > for the iterative Rivara refinement.
> > 
> > I realize the iterative version is an early implementation, but these 
> > are the times I get on my computer:
> > 
> > recursive Rivara: 2.3s
> > iterative Rivara: 416s
> > 
> > So the recursive implementation is over 100 times faster...
> 
> That's really impressive. But I would rather conclude that my 
> implementation is 100 times slower.
> As you've mentioned, I haven't used tricky, dynamic mesh structure, so a 
> lot of time is lost in every iteration for new mesh object creation. 
> Computation of mesh connectivity for edges is also time consuming.
> 
> I'll try to improve efficiency, but I'm afraid that your time is 
> unbeatable.
> 
> 
> > Is it really only the map which makes the iterative version so slow, or 
> > could it be something simple? Perhaps a re-implementation with a dynamic 
> > mesh as in my implementation would speed it up?
> 
> I'm sure you are right, but do you think that it is worth to have two 
> implementation of nearly the same algorithm?
> 
> cheers,
> BArtek
> 
> 
> BTW. What is dolfin standard for real numbers, "double" or "real"? The 
> most of you use "double" but sometimes "real" is used.

It should be "double" anywhere except for in the ODE solvers which
work with both regular double precision or with higher precision using
GMP, as defined in

  dolfin/common/real.h:

We used to have "real" everywhere but then realized that it doesn't
really work with most our external dependencies anyway.

-- 
Anders

Attachment: signature.asc
Description: Digital signature