fortran95 or c for a new stellar evolution program?

[Bill Paxton]

On Nov 24, 2005, at 2:58 AM, Ross Church wrote:

> Apologies for butting into this conversation, but Piet forwarded  
> your emails
> on to me and I thought of one or two things that might be relevent.

No apologies needed, of course.  I'm glad to get your input, and this  
is definitely an "open source" conversation!  So anyone who wants to  
join in is welcome.

>
> On Mon, Nov 21, 2005 , Piet Hut wrote:
>> Thanks for the many comments and strongly expressed opinions!  It's
>> interesting that my preference, c, got essentially no support,

Just to be clear, that quote is from me rather than Piet -- his  
opinion of c was somewhat different!

>
> I might counter that here - if I were wanting to write or use such a
> programme I would prefer the dominant language to be C.  My
> experience of talking to computer scientists and other professional
> programmers has been that they tend not to take Fortran seriously,  
> and were
> you wanting to interact with people outside theoretical physics /  
> astronomy
> a code written soley in Fortran (even f90/95) would possibly be  
> unwise.

Yes, my computer science friends do laugh when they hear I've been  
programming in fortran.  But then I just adopt a superior attitude  
and tell them that f95 is actually quite nice (I don't mention the  
lack of pointers to procedures).

>   The idea of a framework written in a high level language with  
> pluggable
> procedures written in a variety of low-level languages, on the  
> other hand,
> sounds like an excellent idea, and something that could be a very  
> useful
> tool for both teaching and research, regardless of the language(s)  
> chosen.

Frank has pointed out that there are already existing systems, such  
as FLASH, that are multi-lingual in this manner, so it's not a  
completely crazy idea.   There may be a small performance penalty  
(Frank has mentioned 10%-15% for example), but that doesn't worry me  
as much as the possible "complexity" penalty.  I really want to end  
up with something that is easy to download and run, as I mentioned in  
my first message.  And perhaps I should add, easy to understand at  
the level of overall structure -- I may never understand all the  
details of the various physics packages (eos, nuclear burning, etc.),  
but I should be able to grasp the way things fit together without  
having to spend months trying to figure it out.

> It strikes me that one of the big problems in stellar evolution is  
> that
> there are lots of codes that produce similar but slightly different  
> results
> for the same problem, and nobody really knows why this is.  It's  
> very hard
> to compare models between two codes because small differences can be
> magnified into large ones by the non-linearity of the problem.  It  
> strikes
> me that a pluggable code could help us start to address this  
> problem.  As
> such it would be wonderful not only to be able to easily change  
> things such
> as the equation of state or opacities, but also the details of the  
> solving
> procedure.  For instance, it would be interesting to see how the  
> Eggletonian
> (for want of a better term) scheme of arranging a mesh by means of  
> a spacing
> function compares with the more traditional scheme of having  
> Lagrangian
> meshpoints that can be inserted and removed as the code sees fit.

Perhaps this can be addressed by making sure that the stellar  
structure equations are treated as a "pluggable" module just like the  
rest of the physics.  Then one equation package could include the  
remeshing (Eggleton style) while another could leave the mesh  
refinement to a separate operation taking place between steps.  One  
might even hope that there would be an option to replace a quasi- 
static set of equations by a set including hydrodynamics with shocks  
(along the lines of KEPLER, say).

> I think that being able to easily change just one such detail in a  
> code and see what
> its effect is on the resultant models would be a great step forward  
> towards
> quantifying the uncertainties in stellar models that come from the  
> numerical
> method of solving the equations.

You've hit on a point that I think is a critical one -- how can we  
assign "error bars" for the results from our codes?   As a specific  
example, what is the helium core size at the break even point (3  
alpha energy production = neutrino losses) for a 1 Msun, solar  
metallicity star?  Is it 0.47?  Or 0.4 +/- 0.2?

> A note of caution also though - my experience in the past (and I'm  
> sure
> everybody else's) has been that very small tweaks to how a code  
> works can be
> enough to cause it to fail to converge.  How much of a problem this  
> would be
> for a pluggable code when one was trying to change large chunks of  
> code
> arround I'm not sure.  It would be excellent were we to find out!

Perhaps this leads back to the comments from Piet and Jun about the  
need for good diagnostics tools.  We tend to think about  
visualization in terms of the pretty plots of interesting physics for  
the next talk or paper, but an equally critical role is in developing  
and debugging.   The system needs to be built with this is mind.  For  
example, when the solver fails to converge for some puzzling reason,  
wouldn't it be nice to be able to easily get detailed plots of the  
stellar structure it is trying at each step of the iteration as a  
first step toward trying to figure things out?

>
> Best wishes,
>
> Ross (suffering the Eggleton stellar evolution code since 2003)

Ross, Ross -- please, you're supposed to say "enjoying" the Eggleton  
code since 2003!  ; - )
Maybe we need to have a special beer session for all of us who have  
"enjoyed" that particular experience!  I'm sure there are many  
interesting stories to be told.  But for all the pain it may cause,  
it still is truly an amazing piece of work, you must agree.  And it  
will be a real challenge for any new system to match it's robustness  
and performance, let alone exceed them.

Cheers,
Bill