Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

What would you say would be the advantages/disadvantages of two approaches to coding a general (finite volume, fem, dg) library for Computational Continuum Mechanics? This is how I see things right now, so please provide your own experiences and don't flame me for mine :) :

1) C++:

  • generic programming, virtual functions, overloading, speed... : all the genreic + OOP tools available to build whatever you want

  • low level libraries available mostly (no wide spread science&engineering library development such as the one for Python)

2) Python + wrappers for parallel computing (pyOpenCL and others)

  • huge amount of supporting libs of various kinds

  • code what you think: the implementation is done really fast

  • slower execution time

If you wanted to code a framework that would support various methods, work with complex geometries and problems, what would you choose and why?

share|improve this question
1  
I'm not very familiar with pyOpenCL, but generally speaking Python will be much too slow for even moderate-sized problems in 2D or 3D, unless your computational "kernels" are implemented in a low-level language (Fortran, C, etc.) –  David Ketcheson Dec 20 '11 at 17:01
add comment

5 Answers

up vote 12 down vote accepted

I would aim to get the best of both worlds and code the "user interface" (that is, the framework of functions that the user of your library will call to describe the geometry and other properties of the problem) in Python to get the quick turnaround time, then write the simulation run time in C++.

In fact, I would probably mock even the simulation run time up in Python first, then replace it by C++ code piece by piece. Eventually you could consider having your Python code generate C++ source, to be compiled and linked to your runtime online, so that the actual simulation doesn't need to call into Python at all - only return the results at the end. The nice thing about this setup is that it's inherently agile: you start with the quickest and easiest working solution, you'll find out quickly what works and doesn't work, and once you have something you like you can start speeding it up.

(This is how Maple's ODE/DAE solver works, except using Maple instead of Python. Full disclosure: I work for them.)

share|improve this answer
1  
+1. One of the nice bits of Python is the ability to step away from "Pure Python" if needed. –  Fomite Dec 20 '11 at 23:20
add comment

You can also use Cython for your algorithms. It's essentially Python with added type information for some variables that need to be "fast". It translates Python code to C code, which can subsequently be compiled by your favorite C compiler. Careful adding of this type information can make your code up to 150 times faster than naive Python code.

share|improve this answer
add comment

I think there are more to this question. First and foremost a developer will typically prefer what he/her is familiar with unless significant advantages (e.g. in productivity, development time and tools). Personally, I give priority to being productive (time is usually the most scarce resource!) and this favors choices which are close to my experience base.

Perhaps also relevant to take into account are

3) Development time

  • how much time is reserved for development
  • when is results of the work to be delivered? and how?
  • do a code already exist which can do the job? (uniqueness?)

4) Maintenance

  • how many (people) resources are devoted to maintenance?
  • how many people are to work on the code?
  • Is the code to be released at some point? (criteria?)
  • are the code going to rely on third-party libraries?

5) Licensing issue

  • is the code for research?
  • is the code for commercial applications?

6) Productivity and Fun factor (often overlooked!)

  • Where can one be most productive?
  • Where can one have the most fun developing?
  • Any opportunities to be a part of a (social) network?
share|improve this answer
add comment

This depends whether your code can be written as:

some_library_specific_type grid;

for t=0 to T do
    library_function_1(grid,...);
    library_function_2(grid,...);
end

or rather must be written as something like this:

some_home_made_mixture_of_native_types grid;

for t=0 to T do
    for all grid elements as g do
        some_function(g,...);
        library_function(g,...);
    end
end

In the first case choose what you like most to code in; in the second case don't use any scripting language or prepare to suffer from execution time.

share|improve this answer
add comment

As a corollary to Allan's answer (that your own developer time is the most valuable resource): Use what others have already done. You say that you want to develop a library for computational continuum mechanics but there are already several of those that are so large that they will almost invariably already have everything that you need. Take a look at deal.II for example for everything that can be written as a finite element problem, OpenFOAM for fluid dynamics, or PyCLAW/CLAWPACK for hyperbolic problems. deal.II, for example, asks you to program in C++ but in reality the level of programming is often so high that one could say it's like a domain-specific language for FEM codes using C++ syntax.

share|improve this answer
2  
I have never encountered a library that had everything I needed... –  Jack Poulson Dec 21 '11 at 19:57
    
Well, but you get the point I suppose. Some libraries have "almost everything" you may need. To just cite an example I'm particularly familiar with, a finite element solver on fully self-adaptive, 3d meshes running on 10,000+ processors using deal.II and PETSc 126 lines of code. That's clearly more than zero, but it's in fact a very small number given the complexity of what's under the hood. –  Wolfgang Bangerth Jan 2 '12 at 4:57
    
To play devil's advocate, it's trivial to run a code on 10,000 cores, but it's an entirely different matter to make it scalable. Not many parallel preconditioners for non-elliptic equations can even efficiently run on 300 cores. –  Jack Poulson Jan 2 '12 at 6:20
    
Sure. But the example I cite is scalable: math.tamu.edu/~bangerth/publications/2010-distributed.pdf . –  Wolfgang Bangerth Jan 5 '12 at 13:35
    
In the interest of full disclosure: I'm one of the authors of both the paper and the code referenced above, as well as of the deal.II library in general. –  Wolfgang Bangerth Nov 23 '12 at 23:51
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.