Climate model code is so outdated, MIT starts from scratch

It is beyond readability. It is the ability to express ideas with as little code as possible. There was a paper (that I read) about 10 years ago that showed people can effectively write about 300 lines of code per day, regardless of the programming language. So, the more ideas a language can express in 300 lines, the more powerful the language and more productive the programmer will be. Fortran is quite powerful in this regard. Arjen gives the example of one-line sorting algorithm in his Fortran book I believe.
Metaprpgramming and preprocessing would improve this expressivity even further by orders of magnitude.

For a rather extreme example of expressiveness versus readability, consider a program to sieve and print primes up to 100. In APL, the program is

(2=+āŒæ0=(ā³X)āˆ˜.|ā³X)/ā³Xā†100

The corresponding C program is about forty lines, and you can see an explanation of the APL code at SieveOfEratosthenes - APL Wiki .

The APL program above may seem incomprehensible or absurd, but if you wish to see if it actually works, copy and paste the source code (just one line!) into the online Dyalog APL interpreter and press the ENTER key.

1 Like

Is writing a sort in 1 line here meant to be difficult? Something similar to
mysort!(v) = ([(v[i], v[j]) = (v[j], v[i]) for i = eachindex(v) for j = eachindex(v) if v[i] < v[j]]; v)
should work in most languages.
If you want a custom comparison operation, you could use
mysort!(v, o) = ([(v[i], v[j]) = (v[j], v[i]) for i = eachindex(v) for j = eachindex(v) if o(v[i], v[j])]; v)
which can be called as mysort([1,3,2], <) or mysort([1,3,2], >) to sort in forward or reverse order (this accepts an arbitrary ordering function).

1 Like

Students canā€™t read model code

The reason why students can not read the code likely has nothing to do with quality of the language (fortran vs. Julia). It is probably because the old code is terrible for these two main reasons:

  1. The climate model was written in the programming/computational paradigm of 50 years ago
  2. Graduate students have hacked away at the code for 50 years.

I have personally experienced this. In the past, for my work I used a chemical model of the atmosphere written the 80s (link here). Itā€™s honestly a disaster because a mixture of reasons (1) and (2) above.

I decided to re-write the chemical model (link here). I settled on Modern Fortran over Julia. I did this for a few reasons, including:

  1. None of the people I work with know Julia. They know Fortran.
  2. Making Julia code as fast as Fortran is sometimes hard.
  3. I find it easier to write in a strongly typed language.
  4. I knew Fortran better when I started the project.
6 Likes

Here is a Fortran equivalent, which I showed earlier. After defining some general-purpose operators as in APL, it is a two-liner.

module apl_mod
implicit none
interface operator(.i.)
   module procedure irange
end interface
!
interface operator(.d.)
   module procedure drop
end interface
!
interface operator(.o.)
   module procedure outer_product
end interface
!
interface operator (.ni.)
   module procedure not_in
end interface
!
contains
pure function irange(i) result(v)
! return integers from 1 to i
integer, intent(in) :: i
integer             :: v(i)
integer             :: j
do j=1,i
   v(j) = j
end do
end function irange
!
pure function drop(v) result(dv)
! return v(:) without the first element
integer, intent(in) :: v(:)
integer             :: dv(size(v)-1)
if (size(v) > 1) dv = v(2:)
end function drop
!
pure function outer_product(v) result(m)
! multiplication outer product
integer, intent(in) :: v(:)
integer             :: m(size(v),size(v))
integer             :: i,j
do i=1,size(v)
   do j=1,size(v)
      m(i,j) = v(i)*v(j)
   end do
end do
end function outer_product
!
pure function not_in(a,b) result(c)
! return in c the elements of a(:) not in b(:)
integer, intent(in)  :: a(:),b(:)
integer, allocatable :: c(:)
logical              :: amask(size(a))
integer              :: i
do i=1,size(a)
   amask(i) = .not. any(b == a(i)) 
end do
c = pack(a,amask)
end function not_in
end module apl_mod
!
program main
use apl_mod
implicit none
associate (t => .d. (.i. 100) )
   print "(20(1x,i0))",t .ni. [.o. t]
end associate
end program main

output:

 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71
 73 79 83 89 97

Fortran lets you define your own operators, unlike C++ or Python I believe. In C++ you can overload existing operators. For interactive use as in LFortran, abbreviations for transpose and matmul would be convenient.

There is a free 818-page book Mastering Dyalog APL. I wrote a

Fortran program -- compiles with ifort
module m
implicit none
interface operator (.plus.)
   module procedure plus
end interface operator (.plus.)
interface operator (.times.)
   module procedure times
end interface operator (.times.)
interface operator (.in.)
   module procedure in_int
end interface operator (.in.)
interface
   elemental function f_int_int(i,j) result(k)
   integer, intent(in) :: i,j
   integer             :: k
   end function f_int_int
end interface
contains
!
function in_int(x,y) result(x_in_y)
integer, intent(in) :: x(:),y(:)
logical             :: x_in_y(size(x))
integer             :: i
do i=1,size(x)
   x_in_y(i) = any(y == x(i))
end do
end function in_int
!
pure function pos(text,char) result(ipos)
! find the postions of char in text
character (len=*), intent(in) :: text
character (len=*), intent(in) :: char
integer, allocatable          :: ipos(:)
integer                       :: i,nfound,nlen
nlen = len(text)
if (nlen == 0 .or. len(char) /= 1) then
   allocate (ipos(0))
   return
end if
allocate (ipos(nlen),source=0)
nfound = 0
do i=1,nlen
   if (text(i:i) == char) then
      nfound = nfound + 1
      ipos(nfound) = i
   end if
end do
ipos = ipos(:nfound)
end function pos
!
subroutine disp(mat)
integer, intent(in) :: mat(:,:)
integer             :: i
do i=1,size(mat,1)
   print "(*(1x,i6))",mat(i,:)
end do
end subroutine disp  
!
elemental function mult(i,j) result(k)
integer, intent(in) :: i,j
integer             :: k
k = i*j
end function mult
!
function outerprod(i,j,func) result(k)
integer, intent(in)  :: i(:) ! (ni)
integer, intent(in)  :: j(:) ! (nj)
procedure(f_int_int) :: func
integer              :: k(size(i),size(j))
integer              :: ii,jj
do ii=1,size(i)
   do jj=1,size(j)
      k(ii,jj) = func(i(ii),j(jj))
   end do
end do
end function outerprod
!
function tpos(tf) result(ipos)
! return in ipos the positions of the true elements in tf(:)
logical, intent(in) :: tf(:)
integer             :: ipos(count(tf))
integer             :: i,j
j = 0
do i=1,size(tf)
   if (tf(i)) then
      j = j + 1
      ipos(j) = i
   end if
end do
end function tpos
!
pure integer function minus(v) result(z)
! subtract from right to left
integer, intent(in) :: v(:)
integer             :: i,n,w
n = size(v)
if (n == 0) then
   z = 0
   return
end if
z = v(n)
do i=n-1,1,-1
   z = v(i) - z
end do
end function minus
!
pure integer function plus(i,j)
integer, intent(in) :: i,j
plus = i+j
end function plus
!
pure integer function times(i,j)
integer, intent(in) :: i,j
times = i*j
end function times
!
pure function iota(n) result(v) ! p12 built-in APL function
integer, intent(in) :: n
integer             :: v(n)
integer             :: i
v = [(i,i=1,n)]
end function iota
!
elemental function int(tf) result(i)
logical, intent(in) :: tf
integer :: i
i = merge(1,0,tf)
end function int
!
end module m

program main
use m
implicit none
integer, parameter :: n1 = 4, n2 = 2
integer, allocatable :: val(:)
character (len=*), parameter :: fmt_ci = "(a,*(1x,i0))", fmt_cr = "(a,*(1x,f0.2))"
! p5
associate (price => [5.2, 11.5, 3.6, 4.0, 8.45], &
           qty   => [2, 1, 3, 6, 2], vat => 19.6)
associate (costs => price * qty)
print fmt_cr,"costs",costs
print fmt_cr,"sum(costs)",sum(costs)
associate (forecast => reshape( &
   [150,300,100,50,200,330,250,120],[n1,n2]), &
   actual => reshape([141,321,118,43,188,306,283,91],[n1,n2])) ! p8
print fmt_ci,"error",actual-forecast
print "(a,*(1x,f0.4))","VAT tax",price*vat/100 ! p7
print*,"max(75.6,87.3) =",max(75.6,87.3) ! p7
print fmt_ci,"max([11,28,52,14],[30,10,50,20]) =",max([11,28,52,14],[30,10,50,20])
print fmt_ci,"max([11,28,52,14],20) =",max([11,28,52,14],20)
associate (vec => [21,45,18,27,11]) ! p9-10
print fmt_ci,"vec",vec
print fmt_ci,"product(vec)",product(vec)
print fmt_ci,"maxval(vec)",maxval(vec)
val = [22,37,41,19,54,11,34]
print fmt_ci,"val",val
print fmt_cr,"mean(val)",sum(val)/real(size(val))
print fmt_cr,"mean(val)",mean(val) ! using a function
print fmt_ci,"(3 .plus. 6) .times. (5 .plus. 2)",(3 .plus. 6) .times. (5 .plus. 2) ! p11
print fmt_ci,"reduce(val,plus)",reduce(val,plus) ! p11
print fmt_ci,"val(4)=",val(4) ! p11
print fmt_ci,"val([2,4,7,1,4])=",val([2,4,7,1,4]) ! p11
val([3,5,1]) = 0 ! p12
print fmt_ci,"val",val ! p12
print fmt_ci,"val(1:3)",val(1:3) ! p12
print fmt_ci,"val(iota(3))",val(iota(3)) ! p12
associate (salaries => [4225,1619,3706,2240,2076,1389,3916,3918,4939,2735], &
           categories => [3,1,3,2,2,1,3,3,3,2], rates => 0.01*[8,5,2])
associate (salary_inc => salaries*rates(categories)) 
print fmt_ci,"salaries:",salaries
print fmt_cr,"salary increases:",salary_inc
print fmt_cr,"sum of salary increases:",sum(salary_inc) ! p13
print fmt_ci,"salaries > 3000:",int(salaries > 3000) ! p14
print fmt_ci,"# salaries < 2500:",count(salaries < 2500) ! p15
print fmt_ci,"categories == 3 .and. salaries < 4000:",int(categories == 3 .and. salaries < 4000) ! p14
print fmt_ci,"[0,0,1,1] /= [0,1,0,1]:",int([0,0,1,1] /= [0,1,0,1]) ! p14
print fmt_ci,"pack([23,55,17,46,81,82,83],[1,1,0,1,0,0,1]==1):",pack([23,55,17,46,81,82,83],[1,1,0,1,0,0,1]==1) ! p15
print fmt_ci,"pack(salaries, categories == 2):",pack(salaries, categories == 2) ! p15
associate (val => [22,37,41,19,54,11,34])
print fmt_ci,"val:",val
print fmt_ci,"positions with val > 35:",tpos(val > 35) ! p16
print "(a)","outer mult of [5,4,10,3] and [8,5,15,9,11,40]:"
print fmt_ci,"pos('Panama is a canal between Atlantic and Pacific','a'):",pos("Panama is a canal between Atlantic and Pacific","a") ! p16
print fmt_ci,"int([5,7,2,8,4,9] .in. [3,4,5,6]):",int([5,7,2,8,4,9] .in. [3,4,5,6]) ! p16
call disp(outerprod([5,4,10,3],[8,5,15,9,11,40],mult)) ! p23
associate (contents => [12,56,78,74,85,96,30,22,44,66,82,27]) 
print fmt_ci,"contents:",contents
print*,"all(contents > 20), any(contents < 30), count(contents < 30) =", &
        all(contents > 20), any(contents < 30), count(contents < 30) ! p106
print fmt_ci,"minus([45,9,11,2,5]):",minus([45,9,11,2,5]) ! p106
end associate ; end associate ; end associate ; end associate ; end associate 
end associate ; end associate ; end associate
!
contains
!
pure real function mean(v)
integer, intent(in) :: v(:)
mean = sum(v)/real(size(v))
end function mean
!
end program main

to do the calculations found in the first 106 pages with similar brevity and intend to continue.
Q is an APL descendant which can be used with a regular keyboard, and the 32-bit version is free. I want to see how much of Q for Mortals (also videos) can be translated.

2 Likes

I donā€™t think this is an example that will convince people that fortran is ā€œsimpleā€. The main thing that it shows to me is that a 2 line function needs 5 lines of boilerplate.

Only 3 of the boilerplate lines are required, so the main program could be

use apl_mod
associate (t => .d. (.i. 100) )
   print "(20(1x,i0))",t .ni. [.o. t]
end associate
end

Modern Fortran code does have boilerplate. A module may start with

module foo_mod
implicit none
private
public :: foo, bar
contains

before procedures foo and bar are defined, whereas in Fortran 77 you would define foo and bar at the top of the source file. But if the definitions of foo and bar are short, I donā€™t think the earlier boilerplate is a big problem.

5 Likes

My concern about the evolution of modern Fortran is hereā€¦ I read various forum discussions/posts about such user needs and response, but from those pages I am not very sure whether ā€œkey personsā€ (e.g. from the committee side) really understand their need and utility. I guess it might be intuitively difficult for them to understand such needs if they are no longer doing real application programming as daily works and have no opportunity of experiencing how powerful preprocessing / metaprogramming capability helps reduce the amount of coding for usersā€¦

3 Likes

The thing is there are way too many Julia PR pieces floating around the internet to warrant rebuttal one by one. Julia core team runs a for-profit front that needs this kind of exposure, while Fortran is at the other end of the hype spectrum. By now we should all be aware of Juliaā€™s strength and which direction Fortran needs to catch up.

Disclaimer: My first CUDA training is from a group leader in CLiMA project, in C++.

3 Likes

One would think with so much at stake with climate models, more than anything everyone should want

  • clarity and also
  • consistency, especially with the development and use of all the models with all the computational science that has led us here, with everything going at present, and whatā€™s going to come in the future.

Fortran shines in both these aspects. Modern Fortran code authored with a bit of care and discipline strikes about the right balance between readability and expressiveness that leads to healthy clarity. Fads have come and gone but Fortran has persisted and evolved, and it can survive and thrive.

And a few, a handful really, fundamental advances in the language standard as well as a leap in the ecosystem can catapult Fortran to be the lingua franca of core computing needed toward not only climate modeling but also most of other scientific and technical computing.

4 Likes

5 posts were split to a new topic: Language Marketing for Julia compared to Fortran

I think one of the issues with marketing Modern Fortran is there arenā€™t a lot (as compared to Python and C++) of new production level codes that make heavy (or exclusive) use of the new features Modern Fortran brings to the table. One example I point people to when asked about what can Modern Fortran do is the FEMPAR OO Finite Element framework. See

and the source code at

Projects like this should be highlighted in any attempt to illustrate the power of Modern Fortran to a larger audience. Iā€™m sure there are a few others but FEMPAR is the one Iā€™m most familiar with. This community needs to keep a running list of new projects that are using Modern Fortran (almost) exclusively

1 Like

Climate models are so important that it will be great if a set of Fortranners could band together and form a Fortran taskforce whose aim will be to offer modern Fortran translations of such efforts in academia that inevitably adopt programming paradigms du jour. For example, such a taskforce can take copies of the Julia code as they get developed and refactor them into modern Fortran.

First and foremost, this will help cross-check the numerical computing results from such models. There is simply too much at stake with climate modeling, the risk of bad policy decisions that can have tremendous adverse impact is far too great and is a gathering danger. Those who will be affected the most from bad policy decisions will be the weakest life forms on the planet. The cognoscenti - the scientists and engineers - have a profound obligation to check everything and avoid the kind of issues that arose with covid-19 policy making: The Most Devastating Software Mistake Of All Time. Why Is the Imperial Model Under Criticism? Code Review of Fergusonā€™s Model Note here the hubris due to computational advances has reached a stage the very first aspect of do no harm in the Hippocratic oath is often the immediate casualty.

Secondly, it will provide immediate feedback into the computational efficiency of modern Fortran, both its weaknesses and strengths relative to the latest attention grabbers.

Thirdly, it will provide immediate options using modern Fortran for many others in this climate modeling domain who, for whatever reason, are currently in the Fortran space and given their investments and other constraints cannot migrate away from Fortran but would like to evaluate / integrate latest developments in climate models in academia using Julia, etc.

2 Likes

I just added FEMPAR to my list, by request. But it has not been updated in 2 years, and the site http://www.fempar.org/ is gone. There was a 2020 paper by Santiago Badia and Alberto F. MartĆ­n. It may be that the successor project is in Julia (gasp!) GitHub - gridap/Gridap.jl: Grid-based approximation of partial differential equations in Julia. That project is currently active, and there is a paper by Santiago Badia and Francesc Verdugo, both of whom previously worked on FEMPAR. Verdugoā€™s site discusses his projects. It could be interesting to ask Gridap.jl contributors why they migrated from Fortran to Julia.

@Beliavsky, Thanks for pointing out the Gridap project. Actually looks kinda neat. I found the following presentation by Badia that hints at why they moved to Julia.

https://users.monash.edu/~jdroniou/MWNDEA/slides/slides-badia.pdf

A recurring theme in all projects of this type is the claim that programmers can be more productive in Julia (or Python or ā€¦insert your favorite Fortran alternative here). Sadly, thats probably true. I still question if Julia can scale to 400K plus cores and 60 billion unknowns like FEMPAR. Hopefully, the interactive version of LFortran can help with the productivity problem but until there is the wealth of libraries, packages etc that Julia has and something like numpy or the C++ STL for Fortran, Iā€™m afraid Fortran will always lose the productivity argument.

1 Like

A post was split to a new topic: Found contest for articles on Fortran?

Note that Julia has already been successfully scaled to over 600k cores 5 years ago (code: GitHub - jeff-regier/Celeste.jl: Scalable inference for a generative model of astronomical images, paper: Approximate inference for constructing astronomical catalogs from images). The type of communication problems in finite element analysis are almost certainly harder than those Celeste had to deal with, but Julia is a lot more mature than it was in 2018