CFD absractions in Modern Fortran

Question

I'm a novice in object-oriented programming, but this design technique seems very interesting for me. I have a Fortran 90 CFD package , not very complex, but it is already very hard to maintain and to add new physics or other features. As I learnt for now it can reduce complexity of program, but most important from my point of view - intrinsic object orientation of mathematical abstractions. Now I'm trying to draw some parallels between CFD mathematics and object-oriented concepts introduced in Fortran 2003 standard. CFD begins with introducing field abstraction. In most general case it is three-dimensional tensor field, which can be described as five-dimensional array tensor_field(:,:,:,:,:). That array stores nodal values of the tensor field. Similarly , scalar field scalar_field(:,:,:) and vector_field(:,:,:,:). For now everything is clear and straightforward. But when one starts to try to develop classes for this field types it becomes tricky. There are some questions that arise :

1. If that arrays are parts of corresponding classes than how to implement grad and div operations ? In elegant mathematical way grad should be procedure in scalar_field class , that gives vector_field as its result. But similarly , div should be in vector_field class and should result scalar_field. How can be achieved such interactions between objects in object oriented paradigm ?
2. From mathematical point of view tensor_field should be parent class to vector and scalar classes. But it seems a bit overwhelming to use only five-dimensional arrays. I thought about using scalar class as a parental, but it seems more like an element of vector class, not a parent.

To conclude, the main idea of my questions is : are there any object-oriented design pattens for CFD programming , that are wide acceptable and tested ? Maybe you can share here yours implementations in pseudo code, in Fortran or in C++ to discuss them ?

Best regards, Ivan Y.

Answer 1

I don't know if there are design patterns that are valid in the whole CFD field, because there are a lot of different methods available to solve this kind of problem (PDE based, particle based,...). To have concepts that make sense for all these methods, they would probably need to be extremely general, more than what you wrote: a scalar field would be something (probably a method of an object) taking coordinates in input and returning a scalar value (the same for vector/tensor field, with different return types). All this without speaking about "node" or anything related to the method used.

Now, if you restrict yourself to a particular kind of method (e.g. node based), you can more specify your interface in the way you do. If you manage to stick with this kind of interface, you will end up with a very powerful framework which can allow you to treat a wide variety of problems. The first thing I would however change in the interface that you propose is to hide the "array" part, use instead a virtual (abstract) method: you might want to use an analytical function, in which case you probably do not need to store the values of the function in all the nodes.

I coded a few years ago a similar framework for finite elements, which was quite flexible with the definition of scalar/vector/tensor fields: anything that could give a value (scalar/vector/tensor) for a given quadrature node (~position) could be used in the framework. Non-linearity, exact solutions, coulings between equations... could all be mixed and used very easily in the end. I wrote that in C++ and the code can still be found in the ETA module of the LifeV library. It is quite a tough machinery because it relies heavily on C++ template metaprogramming (for performance reasons), but the user interface was very simple and the same could be achieved with simpler code.

For your two particular questions:

1. I would create functions for div and grad outside the field classes, so that there are no problem with circular dependencies. You just write grad(my_field) instead of my_field.grad()
2. I would keep the classes scalar/vector/tensor field separated. There are interfaces that make sense only for scalars and not for vectors, and vice versa. If you need a base class for all the fields, do not specify the storage (array) in the mother class but keep it for particular implementations.

[Disclaimer: I used to be a developer of the LifeV library]

Answer 2

You can see the book by Damian Rousen (Scientific Software Design: The Object-Oriented Way). He uses Fortran 2003 extensively and is one of the author of the Trilinos project. You can also see section 3 here but the book is better and available here.