[naffrancois]

Hello,

I would like to write a very simple implicit scheme to solve Euler equations, backward Euler. Have one of you ever used operator overloading in Fortran to build the Jacobian matrix ? This is the line I am following but:

Second question, how do you usually debug the Jacobian ? I mean a good clue is when you reach bigger cfl than explicit. But if it does not, or if the time integration simply crashes, how do you check you are accurately computing the Jacobian ?

By the way, I am testing the algorithm on Sod shock tube problem

Cheers

[sbaffini]

The book of Damian Rouson:

Scientific Software Design: The object Oriented Way

has specific examples on this fortran topic. Still, i didn't find the book very pleasant at the first read.

[naffrancois]

Thank you very much for the reference, I will try to find it and have a look.

Nevertheless, I had to make a try on an even simpler pde, but I am getting the idea and it looks quite powerful. My long term strong doubt is about performance...

[sbaffini]

This guy:

https://github.com/szaghi/OFF

has developed a basic solver trying to follow the most modern Fortran features. I don't remember the exact details but i remember he complaining about some performance issues when using the most advanced features (probably polymorphicity for type bound procedures). You can find the details here:

http://www.sciencedirect.com/science...10465514001283

Hope it helps

PS There is still an other FV code out there, made by someone who is in contact with Damian Rouson himself. I remember the code being based on several Fortran advanced features. I don't have the link now, but i can provide it in few days. Just to have a comparison.

PPS Here it is the code i was referring to:

http://people.uniroma2.it/salvatore....emo/index.html

and the underlying Ph.D. thesis:

http://www.ce.uniroma2.it/nemo/toninel_phd.pdf

[naffrancois]

Many thanks for your help and links. I was able to apply the concept to the construction of an implicit operator for the 1d euler equations. As I expected, once all the machinery of operator and function overloading is set, it is very easy to switch from a model equation to another. Way much simpler than hand differentiation.

[sbaffini]

Do you have any feel of a sensible overhead (some might be unavoidable)?

[naffrancois]

Well, I did not check the specific computing time spent in the propagation of the derivatives. But, for my goal application which is to solve a sequence of large linear systems (implicit time discretization), the overhead in the end seems to be negligible (<1%). But more tests are coming and extensions to other solvers, I am quite optimist

I am also convinced that it is a much more robust approach, I think so many implicit hand-differentiated codes do not work the way they should just because a small differential mistake was done in the Jacobian. In a Newton-like solver, even a small mistake can spoil the convergence very badly