bugchannelflow: a spectral fluid-flow simulator - Bugs: bug #13389, gmm++ instead of std::vector,...

 
 

bug #13389: gmm++ instead of std::vector, lightweight fft, different class hierarchy

Submitted by:  None
Submitted on:  Mon 13 Jun 2005 03:18:12 PM UTC  
 
Category: NoneSeverity: 3 - Normal
Item Group: NoneStatus: None
Privacy: PublicAssigned to: None
Open/Closed: Open

Add a New Comment(Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

Tue 28 Jun 2005 03:49:01 PM UTC, comment #3:

Hi, Paolo. Thanks for your comments and suggestions. Balancing generality of computation, straighforwardness of user interface, and simplicity in implementation has probably been the trickiest issue in developing Channelflow. I'm not sure how well Channelflow will generalize beyond spectral approximation of Fourier x Chebyshev x Fourier expansions, or degenerate cases of those. The Flowfield class, which is pretty central to everything Channelflow does, is hardwired for Fourier x Chebyshev x Fourier expansions. One could write a similar class for another form of expansion without too much trouble, or one could write a general framework for building such expansions out of elements. My guess if that the latter approach would require years of deep thought and hard work, if done in C++.

Before Channelflow, I worked on a general-geometry, general-representation PDE approximation C++ code, which collapsed under the weight of its own generality. So I decided to err in the direction of specifity towards my problem domain for Channelflow.

I'm also not sure about integration of Channelflow with other packages, particularly ones like ATLAS that make heavy use of sophisticated template constructions. I'm looking into interfacing Channelflow with Python, in order to make Channelflow scriptable and to provide graphics and GUI capabilities. I don't yet have a complete grasp of the issues involved, but it looks like this will be simpler if I keep the template usage low and self-contained, and the object reference/copy model simlpe as well.

On the other hand, I appreciate your interest in working with/on Channelflow, and I don't want to turn you away by rulign anything out. You might have to bear with me while I come to understand your issues.

Other people I'm talking to about Channelflow are mostly interested in parallelization and multithreading.

So with that in mind, responses to your specific comments:

(1) I'd like to hold off on migrating from my home-grown Array to the STL class until I understand the Python integration issues more completely. The Array class is written and it works, so there's no pressure to get rid of it besides cleanliness (which I agree with but it's an aesthetic consideration!).

(2) The linear algebra in core Channelflow code is very simple (banded tridiags) and complete, so I don't see a pressing need for deep integration between Channelflow and GMM++. However, it could be very useful for users to be able to combine the two packages. I wonder if there are any name-space collisions. (Making Channelflow name-space friendly is something I haven't gotten around to yet.)

I've been combining Channelflow and Octave linear algebra libraries in my research codes, but I've found that binary Octave packages (RPMs) always have missing function definitions and I've been forced to tinker with the source packages to make it work. It's tricky, and so I haven't supported this combination to other users. I wonder if GMM++ integration would be easier.

(5) I want to keep the object memory-management model as simple as possible, for the sake of future Python integration and because my experience with reference-counting smart-pointers code was painful and bitter.

(6) I'm planning to rename ChebyCoeff to ChebyExpan or ChebyExpansion. I plan to implement a number of your suggestions for it, like building in a quick-and-dirty FFT, and moving the gridpts function to a member function. But I'm not inclined to generalize it beyond finite-domain Chebyshev expansions, because it won't pay off for my problem domain, and it seesm it would it would complicate the code significantly.

Can you implement infinte-domain expansions as a separate class? Would this satisfy your needs? If not, how would you integrate them, and at what cost?

(7) Doxygen would be great. I've concentrated on thea Latex manual since it's good for mathemtical expressions. But Doxygen would be good for a programmer's reference manual.

ok, thanks again for your thoughts!

John

John F. Gibson <gibson17>
Project Administrator
Mon 20 Jun 2005 07:45:54 AM UTC, comment #2:

(1) std::vector<double> access method; I used at() because I am lazy; operator[] give the fastest access from within the class; the best thing probably is to write a compact inline access method as you say, maybe:
inline double Distribution1::AT(int n) {
#ifdef VERIFY
return at(n)
#else
return operator[n];
#endif;
}

(2) GMM++ has precisely the vector class with mathematical vector methods you mention. Registering the Distribution1/ChebyCoeff as a vector for GMM++ allow one to access norms, addistion, scaling ...
Provided the data layout is FORTRAN-like fast BLAS-like algorithms (ATLAS...) can be turned on if required, achieving better performance than any self-written code.
It is interesting to see that GMM++ has in some places bound checking enabled/disabled through the macro definition GMM_VERIFY

(4) 2d FFT; to be honest at the moment I'm more in need of 1d transforms so I have not looked into the 2d code ...

(5) Domain pointers; you are perfectly right in pointing out the risk of segmentation faults; this same problem arises in Windows proramming if you free a Handle and then use it again; personally I a not very scared, but for the general public you could use the facilities modern C++ offers to avoid those problems, i.e. shared_ptr.

Additional tought: one desireable feature is resampling, like if I resize the Domain, all Distribution1 based on that resize their "data_" (and interpolate if Physical or truncate the series expansion / leave the new terms 0 if Spectral).
This could be implemented having a std::list<Distribution1 *> in the Domain, so that it knows how to deal with that.
What do you think ?

Final comment: I plan to use portions of your code (at the moment from the file "chebyshev.cpp") modified along the lines I have shown you within a commercial software development project. I will in that case redistribute the modified source code portion and notify you of any issue I may find with the original code, according to the GPL.
I would be happy to contribute my bit as a developer for channelflow, but it depends how much of my ideas get into it ...
The point where I most agree with you is where you write "Channelflow falls short of a providing a language for spectral simulation, due to the scope of the problem domain and to the difficulty of presenting a clean syntax through C++ class
libraries." I have looked at a few C++ projects for FEM or pseudospectral, and invariably they were trying to misuse the language following the attractive dream of copying the equations from a book into C++; better clean and fast and let the user write her proper for-loops.

Having said that, all my comments go in the direction of widening the scope a bit.
In fact the classes you develop for "incompressible fluid flow on a periodic, rectangular, wallbounded domain" with a little adptation can be used for any problem where pseuso-spectral methods are applicable.
Calling ChebyCoeff "Distribution1" paves the way to implementing different transforms or applying it to infinite domains or whatever.
And of course for documentation with little effort and Doxygen we can make a much better job (and achieving the dream of literate programming ...).

Anonymous
Wed 15 Jun 2005 05:38:50 PM UTC, comment #1:

Thanks! I have long meant to include a quick-and-dirty FFT so that Chebyshev transforms can be done without constructing an FFTW plan. I also agree that the cheybpoints function should be eliminated, in favor of Vector ChebyCoeff::gridpoints() and Real y ChebyCoeff::gridpoint(int i).

I'm wary of your plan to have object X contain an internal pointer to object Y. If Y goes out of scope before X, this produces a hanging pointer and potential seg fault.

The ChebyCoeff/ChebyTransform (and most everything else) do have a lot of assertions to check things like objects' N, a, and b match. The checks are turned on in the debugging libraries.

The main reason for distinct ChebyCoeff and ChebyTransform classes is the relatively high cost of FFTW plan generation. In many places in the code, one ChebyTransform is built, optimized, and used repeatedly on many ChebyCoeff objects. Getting rid of chebypoints will eliminate the possibility of errors due to changes in a,b. But I think coordinating N between ChebyCoeff and ChebyTransform will remain necessary, for efficiency's sake.

I'll try to incorporate your FFT in the next release.

John

John F. Gibson <gibson17>
Project Administrator
Mon 13 Jun 2005 03:18:12 PM UTC, original submission:

In channelflow for 1-D distributions there are currently three distinct objects:

ChebyTransform
chebypoints
ChebyCoeff

When needed, these objects are connected together with ad-hoc instructions:

ChebyTransform t(N);
ChebyCoeff f(N, a, b, Physical);
t.chebyfft(f);
Vector y = chebypoints(N, a, b);
for (int n=0; n<N; ++n)
f[n] = my_func(y[n]);

But what happens if a and b between the call to the ChebyCoeff constructor and the call to chebypoints change ? Or if N changes ?
Everything screws up or you should add a lot of checking code !

The alternative is to assemble higher-level objects and rename as follows:
ChebyTransform + chebypoints = Domain
ChebyCoeff = Distribution

The code becomes:

Domain t(N, a, b); // initialize fftw etc.
Distribution f(t, Physical);
f.chebyfft(); // has a pointer to the domain so calls fftw_execute but using plan
for (int n=0; n<t.N(); ++n)
f[n] = my_func(t.y[n]);

In the attached ZIP I have tested this approach.
Note that it is still a prototype, pretty sure not all the costructors work propertly (what is not tested does not work...).

There is also the following additional stuff:
- fast evaluation of chebyshev Tn @x and its derivative
- MMT now part of domain_MMT object and can be used as (painful) alternative to FFT (to use only for the sake of cross-checking)
- integration with GMM++ (http://www-gmm.insa-toulouse.fr/getfem/gmm_intro) for matrix manipulation rather than use std::vector<double>
- -DUSE_FFTW compile switch required to turn on FFTW3 code, otherwise it will default to a very quick and dirty R2CFFT derived from Public domain code for computing the FFT contributed by Christopher Diggins, 2005

Profiling the quick and dirty FFT against FFTW for N=1048577 samples yields 5.2 s against 1.8 s with FFTW_ESTIMATE, or 0.5 s if FFTW is using wisdom previously generated with FFTW_MEASURE. This is the order-of-magnitude speed-up reported in literature.

To compile, get fftw-3.0.1 and gmm-1.7.

paolo.greppi at nospam.tiscali.it
(remove nospam.)

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):
   
   
Comment:
   

Attached Files
file #3084:  domain.zip added by None (13KiB - application/x-zip-compressed)

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

Do you think this task is very important?
If so, you can click here to add your encouragement to it.
This task has 0 encouragements so far.

Only logged-in users can vote.

 

Please enter the title of George Orwell's famous dystopian book (it's a date):

 

 

Follows 1 latest change.

Date Changed By Updated Field Previous Value => Replaced By
Mon 13 Jun 2005 03:18:12 PM UTCNoneAttached File-=>Added domain.zip, #2602

Back to the top


Powered by Savane 3.1-cleanup1