33
Ease of learning
Python and Fortran are both relatively easy-to-learn languages. It's probably easier to find good Python learning materials than good Fortran learning materials because Python is used more widely, and Fortran is currently considered a "specialty" language for numerical computing.
I believe the transition from Python to Fortran ...
27
It's a bit of a popular misnomer that there is a "version" of Fortran to know. With rare exception, the latest Fortran standards (and compilers) retain excellent backwards compatibility with older standards. This is with good reason: not many people would use Fortran today if it weren't for the large amounts of legacy code still in use. That is to say, a ...
14
To make more robust comparisons (on linux), you can :
1) On Intel CPUs the turbo overclocks your CPU. This is controlled by the temperature of the CPU, so it can behave differently from one run to the other. On Linux, you can block the frequency of the CPU as follows. For example, for 2.4GHz:
echo 1 > /sys/module/processor/parameters/ignore_ppc
for ...
12
Fortran's allocatable variables are automatically deallocated when the variable goes out of scope (see http://www.fortran90.org/src/best-practices.html#allocatable-arrays). This means that it is not possible to create a memory leak by failing to deallocate an allocatable array. This is one of the big benefits of using allocatable arrays rather than pointers.
...
11
Python is a very slow, high level language. For fast number crunching you'll have to write the main compute kernels in low level languages like C/C++ which means that now you have to learn not one but at least two languages. You'll also have to deal with additional headache associated with debugging/installation/maintenance etc. Most people use Python as a ...
10
When one uses a low–level programming language, e.g. C++ or FORTRAN, one essentially controls lots of things: how parameters are passed, how data structures are aligned in memory, what is the most efficient way to loop over elements of a big sparse matrix (see cache thrashing) when one multiples it, and so on. In contrast, when one uses high–level software, ...
9
Java has been around for almost 20 years now as a major programming language, but it hasn't caught on in scientific computing so far. I think that's a good indicator for what's going to happen in the future.
My take is that the issue isn't speed. Most people are probably willing to give up 20% of performance (or even a factor of 2) if they would be vastly ...
9
There is no built-in Fortran functionality to do linear interpolation. You could either use a library or write your own routine. I haven't tried compiling or testing and my fortran may be a bit rusty, but something like the following should work.
subroutine interp1( xData, yData, xVal, yVal )
! Inputs: xData = a vector of the x-values of the data to be ...
9
I think that the problem is linked to the way in which f2py generates the fortran interface: the argument to fortranrun.f2py should be stored as a F_CONTIGUOUS array, otherwise the interface will create an internal copy with the correct storage order.
Python 3.6.2 (default, Jul 22 2017, 21:19:22)
[GCC 7.1.1 20170516] on linux
Type "help", "copyright", "...
9
You should consider giving Julia a try. Let me explain what's going on in the design space right now that would be of interest to you. Full disclosure I am the lead developer of JuliaDiffEq.
JuliaDiffEq and DifferentialEquations.jl has a large feature set dedicated to efficiently integrating computationally-difficult differential equations. It has a simple ...
8
I would stay away from Fortan, or if you must, use a reasonably new version (2003 rather than 77). A lot of physics software (Monte Carlo simulations in particular) is written in Fortran, simply because the projects were originally started in the 80s.
That being said, python and Fortran are two very different language, and what they should be used for is ...
8
In general, you will be much more productive writing software in a higher-level
language (e.g MATLAB) that has features useful in describing problems in your particular
domain (e.g. matrices in computational science). Often much more time is spent
writing the software rather than acutally running it so reducing programming
time by using the right higher ...
8
I think it's generally true that there are no advantages of Fortran 77 over either newer versions of Fortran or in fact any number of other programming languages that are widely used in scientific computing.
The reason it's still used is because there are millions of lines of code around that are written in Fortran 77. Now, recall that it takes a good ...
7
There is almost no good reason to write your own dense matrix manipulation routines until you have compared to fast libraries (MKL, FLAME, MAGMA, etc). Writing these libraries is challenging work that depends highly on the target architecture and requires deeper concepts than naive nested do loops. There's a whole literature on the subject that you should ...
7
Shoutout to Kyle Mandli and Endulum who each contributed to this answer in the comments.
First, I took Endulum's suggestion and removed the redundant reshapes. After this change the Fortran version was beating Python on small scale examples, but at scale the Python version was still faster. Then I implemented Kyle Mandli's suggestion and replaced all of the ...
6
Disclaimer: I am the author of the linked modules.
First I would recommend to post this question at Stackoverflow, many Fortran experts meet there.
There are more possible approaches to this problem. One similar to C void pointers (using transfer) or safer one using class(*). Both theoretically allow heterogeneous lists. You are then responsible to ...
6
This is probably not the answer you are looking for, but I wanted to state it anyway: Your choice of programming languages introduces two difficulties you will encounter as your program grows.
First, you will find that there is a cost to trying to couple different languages. They may interoperate, but there is always some quirk to it and as your software ...
6
It seems unlikely to me. The Java MPI APIs haven't been worked on in years (so you're wrong about #4), and the JVM's floating-point performance is notoriously poor. Java may out perform C/C++ or Fortran in some areas due to rapid thread creation and easy memory management, but these aren't the bottlenecks in typical scientific programs.
As to your #5, the ...
6
I would argue that Java will in fact REDUCE productivity when compared with modern c++, or even with modern Fortran for the purpose of scientific computing.
Writing
A = B*C+2*D
is just so much more readable than
A = B.mult(C).add(D.mult(2))
Assuming the code above deals with arrays, both C++ and fortran will also produce significantly more efficient ...
6
There are three issues that are likely to cause such problems in pseudospectral methods:
Gibbs oscillations
Aliasing
Time step too large
In any case you likely develop oscillations in the solution until some point ends up with a negative density, resulting in a NaN when computing the pressure or sound speed or some other term. The solution to 3 is obvious, ...
6
One issue causing the jagged spectra at high wavenumbers is under sampling there. For example consider the 2D analogue of your binning procedure:
You don't want to sample from the red zones as they will become increasingly under-sampled as you move past a radius of size $|k_{x}|=|k_{y}|$.
6
Typically, you would add "guard cells", that is (for u) u(-1) and u(n+1) with your notation. Before each integration step:
u(n+1) = u(0)
u(-1) = u(n)
and similarly for the other variables. If you use higher order derivatives, you could also define u(-2) and u(n+2), etc.
5
Python
is very practical for full simulation analysis with well-documented versatile packages: grid generation, array computation and data structure handling (numpy and pandas) as well as data visualization with matplotlib.
For complex simulations with big result files, it's even better to work with the VTK package which allows exporting data to be read by ...
5
An MPI_Barrier can be used to synchronize all processes in a communicator. Each of the processes has to wait till all other processes reach the barrier before all of them can proceed further:
MPI_BARRIER(COMM, IERROR)
Here COMM is the communicator handle and IERROR the error status.
Note that sometimes there are more clever ways to deal with such ...
5
The clear answers if you want to keep Fortran speed are to use a language which has proper code generation like Julia or C++. C++ templates have already been mentioned, so I'll mention Julia's tools here. Julia's generated functions let you use its metaprogramming to build functions on demand via type information. So essentially what you can do here is do
@...
5
For the record, you might try putting an OpenMP worksharing region around your array operation syntax:
!$omp workshare
C=A+B
Don't forget to build with OpenMP enabled (-openmp for the Intel compilers) and to set OMP_NUM_THREADS.
5
Honestly, the best ideas have already been said. I'll try to synthesize my thoughts anyway.
First, the best way to write a program is whatever way gives you the results you need, for the least development time. For some applications, you NEED bare-metal performance, but a lot of the time you don't. Since you've been working in Octave up until now, you may ...
5
For a sparse parallel solver, it's your own responsibility to provide a matrix vector product and a suitable preconditioner. The data for the vector itself should fit into main memory in any case. If the matrix has at most a fixed small number of non-zero elements (<20) per column or row, then the same is also true for the matrix itself. In this case, an ...
5
I'd say Gmsh. I used it for a few finite element projects, and it was mostly easy to work with. The mesh output formats are very parseable, and there's at least one third-party parser (MeshPy) that can parse the output. It also has a C++ API, and the mailing list gets enough traffic (probably 10-20 messages a week) that your questions might be answered (in ...
5
Notice that the summation in $U^{(r)}$ is incorrect. You want to sum over all the copies of the atoms in the lattice of periodic boxes, not just those whose indices satisfy $i > j$. In the original box, of course you want to avoid self interactions $i = j$ but just in that one box. In other words, do not discard the electrostatic interactions between ...
Only top voted, non community-wiki answers of a minimum length are eligible