Need Help converting a Fortran CFD code into C++ (Using MS Visual C++ 2008 Express)

[Thunder Chicken]

Hello all,

Back in university I wrote a flow solver as part of a graduate CFD (computational fluids dynamics) course. This was written in Fortran 90 on a Linux machine, I forget what compiler I used. It is simple main program that utilizes several subroutines (one to allocate the memory, another to create the mesh of control volumes, another to initialize the variables, calculate divergence, etc..). It was nicely modular, but it still had some bugs.

I want to resurrect this code and rewrite it in C++. I am trying to learn C++ and the class structure is giving me fits - every 'explanation' I read is thoroughly non-intuitive and I am left more confused. I have written modules for Orbiter by religiously copying everything in the SDK manual, but this is my first unguided foray into trying to make a C++ console app (it just needs to read some inputs from a file, run some calculations, and write out a data file, graphics and GUIs are not needed).

I *think* the way to go is to translate the subroutines into C++ code into individual header files, and then call them from a main cpp code (class?). Can anyone put me on the right path? I get along OK with the C++ language (formulas and such), it really is just the class structure and related code relating base and main classes and headers that I can't get straight in my head.

I am using MS Visual C++ 2008 Express Edition if that helps.

 

[tblaxland]

If it was written in Fortran, you shouldn't need classes to do a fairly straight translation, so the code will probably look more like C than C++.

It sounds like you are getting confused between classes and the file structure. A class is a special construct to that is used to define objects (often an abstraction of a real world object). Classes are generally declared in a header file and defined in a code file.

Header files (.h) and code files (.c/.cpp) have nothing to do with classes directly. C (which has no classes) also uses the same .h/.c file structure. The idea is that your function and type declarations are done in the header file, whereas the definitions are done in the code file. It is actually possible to write a program without a header file and just a code file but you need to be careful because you cannot call a function in your code before you have declared it. See section "Declarations and Definitions" in this document for the difference, understanding this is quite important to understanding the file structure. A common method would be to put your subroutine declarations in a header file (MySubroutines.h), their definitions in a code file (MySubroutines.c) and your main program in a different code file (MyProgram.c) that has a line:

#include MySubroutines.h

at the top.

Some features of C++ may be useful for your translation if you feel limited by C. I'm thinking of the Standard Template Library, in particular std::vector and std::string, and function overloading. Defining your own classes could also be useful but it would be beginning to look like a re-write of your code rather than a translation.

HTH...

 

[Thunder Chicken]

I've actually used C a fair amount (one commercial CFD code that I have used takes user defined functions written in C). The .cpp/.h arrangement is pretty much how I had the original program arranged, except that the header files each contained a subroutine, and so the include for these was done in-line.

As you say, I think I am confused by the class terminology and really had no feel or understanding as to what it was doing or whether it was needed.

Let me sketch up a template .cpp and .h and I'll get back to you.

Thanks for the lead!

 

[mjessick]

If you really just need to run the FORTRAN code, you could buy a new empty hard drive (unformatted), install Linux so that you can boot either operating system ("dual boot"), spend 6 months getting familiar with Linux and you are done. It shoudln't take too long to actually compile and link the FORTRAN progam. The "make"-ing programs (build scripts) and compilers seem to be more standardized than anything else in Linux land.

It comes down to whether you want to learn useful programming languages like C and C++. If you don't have enough example cases to test your new code, you might want to both anyway soyou can run the same problem in both FORTRAN and C.

The "Numerical Recipes in C" book used to have a good section on the challenges of writing numerical algorithms in C from a FORTRAN perspective. And also very valuable in general when converting FORTRAN codes because they provide systems that convert the array indices easily. (FORTRAN indices go 1,2,3 while C arrays start at zero.)
________________
Numerical Recipes in C, second edition. Available free online at: http://www.nr.com/oldverswitcher.html with some other older books.
direct link: http://www.nrbook.com/a/bookcpdf.php

Section 1.2 Some C Conventions for Scientific Computing

 

[Thunder Chicken]

OK, here is the main program and one of the subroutines (functions) that I need to call from the main program.

The main program is pretty much just this for now:

#include "Divergence.h"

int main ()
{

//constants defined here

const int ISIZE = 100;
const int JSIZE = 10;

//memory allocated in here

double  dx, dy;
double u [ISIZE] [JSIZE-1];
double v [ISIZE-1] [JSIZE];
double div_u [ISIZE-1] [JSIZE-1];

//bunch of loops to initialize variables here

.
.
.

return 0

}

Here is the Divergence.h - It just calculates the divergence of a velocity vector field:

double Divergence( int ISIZE, int JSIZE, double u, double v, double dx, double dy)
{

    int i, j;

    for (j = 1; j < JSIZE-1; j++)
    {
        for (i = 1; i < ISIZE-1; i++)
        {
            div[i][j] = (u[i+1][j] - u[i][j])/dx + (v[i][j+1] - v[i][j])/dy;
        }
    }
}

The loops in the divergence routine look correct to me. I'm clearly passing u and v incorrectly, but I haven't hit on how to pass the array variables. Also, the function Divergence needs to return an array of double, so that is clearly wrong as it is.

The good news is that the main program can see the header containing the function, so far so good. Can someone help me with the syntax required to pass the array variables?

 

[EliNaut]

Small off-topic comment here..

I have much sympathy for you, Thunder.
I attempted learning C++ from a Dummies book, as the series has helped me with other things before, and I was completely dumbstruck when it all of a sudden introduced classes. Thats what made me backburn it and eventually lost patience with it as I cycled around with my regular hobbies. =/
I want to make a more headstrong attempt sometime soon, hopefully.

 

[tblaxland]

The loops in the divergence routine look correct to me. I'm clearly passing u and v incorrectly, but I haven't hit on how to pass the array variables. Also, the function Divergence needs to return an array of double, so that is clearly wrong as it is.

This article probably explains it better than I could:
http://www.eskimo.com/~scs/cclass/int/sx9b.html

Another minor point: C/C++ arrays start at index 0, not 1.


-----Post Added-----


I've been thinking about this some more and I wonder if you might be better served by using a std::vector instead of C-style arrays. std::vector is a special template class for storing lists of objects - don't let that scare you, they are pretty easy to use. For a two dimensional array, you would use a vector of vectors. For example, I would re-write your code like this:

(I noticed in Divergence that div is not defined anywhere. I have assumed that it is meant to be the return value and that you want it stored in div_u)

 #include <vector>
#include "Divergence.h"

int main ()
{
  //constants defined here

  const int ISIZE = 100;
  const int JSIZE = 10;

  //memory allocated in here

  double  dx, dy;
  std::vector<std::vector<double>> u(ISIZE,std::vector<double>(JSIZE-1));
  std::vector<std::vector<double>> v(ISIZE-1,std::vector<double>(JSIZE));
  std::vector<std::vector<double>> div_u(ISIZE-1,std::vector<double>(JSIZE-1));

  //bunch of loops to initialize variables here
  // ...
  Divergence(div_u,u,v,dx,dy);
  // ...
  return 0;
}

And Divergence.h:

void Divergence(std::vector<std::vector<double>> &div, const std::vector<std::vector>> &u, const std::vector<std::vector>> v, double dx, double dy)
{

  int i, j;

  for (j = 0; j < (v[0].size() - 1); j++)
  {
    for (i = 1; i < (u[0].size() - 1); i++)
    {
      div[i][j] = (u[i+1][j] - u[i][j])/dx + (v[i][j+1] - v[i][j])/dy;
    }
  }
}

There a couple of interesting techniques thing used here:

First is "pass by reference". This means that the div,u,v in Divergence only have a local scope but they are aliases for the div_u,u,v in main(). This means that assignments to div actually occur in div_u and means that you don't have to copy large amounts of data on returning from the function. That is the function of the "&"s in the Divergence argument list.

Second thing is the use of "const" in the Divergence argument list. Because u,v are aliases of the u,v in main(), Divergence could potentially corrupt that data if you made a mistake in your code. "const" protects against that by making sure you can only read from u,v not write to them. This is known as "const correctness" and is not necessary but is definitely a Good Thing.

I don't have time to go into std::vector just now but there is a good tutorial here:
http://www.codeguru.com/cpp/cpp/cpp_mfc/stl/article.php/c4027

and an article that explains vectors of vectors here:
http://www.tek-tips.com/faqs.cfm?fid=5575

 

[Thunder Chicken]

Another minor point: C/C++ arrays start at index 0, not 1.

Yeah, I am aware of this. The FORTRAN code has a bunch of loops that index starting at 1, I'll have to sift through carefully and fix them all.

I've been thinking about this some more and I wonder if you might be better served by using a std::vector instead of C-style arrays.

Maybe, but I think one thing I dislike about C++ is the visual clutter you get in the code when you start using these classes - simple things start to look baffling after a very short while. Fortran may not be as powerful as C++, but a child who has coded in BASIC could understand what was going on in a Fortran code. I'd like to see if I can keep that simplicity, if I can help it.

There a couple of interesting techniques thing used here:

First is "pass by reference". This means that the div,u,v in Divergence only have a local scope but they are aliases for the div_u,u,v in main(). This means that assignments to div actually occur in div_u and means that you don't have to copy large amounts of data on returning from the function. That is the function of the "&"s in the Divergence argument list.

Second thing is the use of "const" in the Divergence argument list. Because u,v are aliases of the u,v in main(), Divergence could potentially corrupt that data if you made a mistake in your code. "const" protects against that by making sure you can only read from u,v not write to them. This is known as "const correctness" and is not necessary but is definitely a Good Thing.

Very helpful points. Thanks for your help!

 

[mjessick]

Yeah, I am aware of this. The FORTRAN code has a bunch of loops that index starting at 1, I'll have to sift through carefully and fix them all.

Or use the Numerical Recipes techniques and leave the arrays 1 based, at the cost of one unused [0] element per array.

Frankly, when I convert FORTRAN code I care more about not hand changing the cryptic indices than about the one unused memory space.

 

[movieman]

Maybe, but I think one thing I dislike about C++ is the visual clutter you get in the code when you start using these classes - simple things start to look baffling after a very short while.

Typedef is your friend.

 

[tblaxland]

Maybe, but I think one thing I dislike about C++ is the visual clutter you get in the code when you start using these classes - simple things start to look baffling after a very short while.

Typedef is your friend.

As are wrapper classes. Check out the link I posted earlier, sections titled "Vector-Vector with wrapper class" and "Custom class dynamic_2d_array":
http://www.tek-tips.com/faqs.cfm?fid=5575

 

[Thunder Chicken]

I guess my issue with the vector and other classes is that I don't see the advantages in using them. Renaming variables, using vectors instead of array variables, etc.. what does it buy me except more complex-looking code? What are the advantages over simply using an array variable (e.g. u [isize][jsize]) and simply calling it u all of the time?

 

[tblaxland]

I guess my issue with the vector and other classes is that I don't see the advantages in using them. Renaming variables, using vectors instead of array variables, etc.. what does it buy me except more complex-looking code? What are the advantages over simply using an array variable (e.g. u [isize][jsize]) and simply calling it u all of the time?

std::vector has a number of useful functions (maybe not useful for your particular application) such as the ability to easily insert, delete and sort its elements. They also are more secure in that the size of the vector is stored with the class and managed automatically so it is harder for you to accidentally run into "out of range" errors. The class also handles memory allocation/deallocation for you. Using a wrapper class makes thing pretty neat too. Compare:

dynamic_2d_array <double> MyMatrix(x, y);
double MyMatrix[x][y];

Not too bad is it? Also, resizing that array is effectively impossible. You really need to make it bigger than you ever expect to need, using unnecessary amounts of memory in the process, and then manually manage a size variable or end pointer, which is fraught with potential errors. Compare that with MyMatrix.resize(x,y). Easy!

If you don't need those features, then I see no problem use the pointer-pointer method I linked to earlier. The code for that is pretty clean except for the memory allocation-deallocation parts, but I don't see that it is any cleaner than the vector-vector-wrapper method, with none of the benefits.

BTW, what renaming of variables are you referring to?

I'm not too familiar with arrays in other languages, but do they offer features similar to a vector?

 

[Thunder Chicken]

I guess it boils down to application. In CFD, dynamically allocated arrays pretty much handle everything without a problem. We're working with fluids, 3 array indices at most. Maybe better heads than mine can utilize all these c++ features, but I am all for simple, dirt simple, code.

The one question that I still have is how to write an array function. My divergence function needs to return an 2D array of float. How do you specify that?

 

[tblaxland]

The one question that I still have is how to write an array function. My divergence function needs to return an 2D array of float. How do you specify that?

I did that in my earlier post (#7) already. The first argument to the Divergence function is a reference to the 2d array (vector-vector in that case) where the result is stored. This same principle could be applied to a pointer-pointer instead. I can do it for you if you have trouble (let me know), but vector-vector is better, IMHO.

The function could return a vector-vector also (instead of storing the result by reference), but you would then be just be copying it into a variable with main() local scope anyway. Save yourself the memory cycles

 

[Thunder Chicken]

I did that in my earlier post (#7) already. The first argument to the Divergence function is a reference to the 2d array (vector-vector in that case) where the result is stored. This same principle could be applied to a pointer-pointer instead. I can do it for you if you have trouble (let me know), but vector-vector is better, IMHO.

The function could return a vector-vector also (instead of storing the result by reference), but you would then be just be copying it into a variable with main() local scope anyway. Save yourself the memory cycles

Can't you reference arrays like this?:

void Divergence(int isize, int jsize, &div[isize][jsize-], u[isize-1][jsize-1], v[isize-1][jsize-1], double dx, double dy)
{

  int i, j;

  for (j = 0; j < (jsize - 1); j++)
  {
    for (i = 0; i < (isize - 1); i++)
    {
      div[i][j] = (u[i+1][j] - u[i][j])/dx + (v[i][j+1] - v[i][j])/dy;
    }
  }
}

The loop indices are all kludged up yet and I can't test this until that is fixed. Do you need a return statement, or is that handled by the '&'? This is almost the same way Fortran handles arrays and so it would not be such a stretch to convert the code if this does work.

Thanks for your help.

 

[tblaxland]

Can't you reference arrays like this?:

void Divergence(int isize, int jsize, &div[isize][jsize-], u[isize-1][jsize-1], v[isize-1][jsize-1], double dx, double dy)

No. The array indexes cannot be variables when used in a function argument list. You can use that syntax if the array indexes are macros, eg:

#define ISIZE 100
#define JSIZE 10

void Divergence(double div[ISIZE][JSIZE-1], const double u[ISIZE-1][JSIZE-1], const double v[ISIZE-1][JSIZE-1], double dx, dobule dy)

This is because C interprets div,u,v as pointers to the original array, it does not make copies of the arrays for local use. It then needs to know, at compile time, how much to offset those pointers before dereferencing them. Because div,u,v are pointers, you do not need the "&" in front of div in order to access the original array that was allocated in main().

 

[JamesG]

If you haven't already, it will probably be worth your while to read up on object oriented design/programming a bit. Just google OOP and there are quite a few online tutorials.

Once you understand the way C++ etc. is intended to be used, you won't keep trying to beat your square structural block into an object oriented hole.

 

[tblaxland]

Can't you reference arrays like this?:

void Divergence(int isize, int jsize, &div[isize][jsize-], u[isize-1][jsize-1], v[isize-1][jsize-1], double dx, double dy)

A clarification on my last post. The C99 standard supports the above syntax, but it is not implemented by Visual C++. GCC implements it though, so maybe you want to compile with that?
http://gcc.gnu.org/onlinedocs/gcc/Variable-Length.html

 

[Thunder Chicken]

Ahhh, now it is coming together. The code is converted and compiles, and the executable runs (but is giving me garbage numbers, just like the old Fortran version). Now it is a matter of sorting out the code and triple checking all of the indices.

It is a tough code to keep straight, as their are several flow-related variables (u and v velocities, pressure, and others) that do not share the same array bounds. The solver, like many incompressible codes, uses a staggered grid where velocities are defined at the faces of cells, while pressure is defined at the cell center.

Anyway, it is an algorithm problem now, I seem to have enough handle of C++ to run with this now. Thanks for your help!