I recently found some old code I had used to solve the Poisson equation with a somewhat naive double loop relaxation method on a 2d grid. The C code is as follows:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
void swap(double a[101][101], double b[101][101], int n)
{
double swp[101][101];
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
swp[i][j] = a[i][j];
a[i][j] = b[i][j];
b[i][j] = swp[i][j];
}
}
}
double mmax(double phi[101][101], double phip[101][101], int n)
{
double max = 0.0;
double diff = 0.0;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
diff = abs(phi[i][j]-phip[i][j]);
if (diff > max)
max = diff;
}
}
return max;
}
double rho(double x, double y)
{
double s1 = 0.6;
double e1 = 0.8;
double s2 = 0.2;
double e2 = 0.4;
if (x > s1 && x < e1 && y > s1 && y < e1) {
return 1.0;
} else if (x > s2 && x < e2 && y > s2 && y < e2 ) {
return -1.0;
} else {
return 0.0;
}
}
void run(double toler, double a)
{
double epsilon0 = 8.85e-12;
double phi[101][101];
double phip[101][101];
memset(phip, 0, sizeof(phip[0][0])*101*101);
memset(phi, 0, sizeof(phi[0][0])*101*101);
double delta = 1.0;
while (delta > toler) {
for (int i=1; i < 100; i++) {
for (int j=1; j < 100; j++) {
phip[i][j] = (phi[i+1][j] + phi[i-1][j] +
phi[i][j+1] + phi[i][j-1])/4 +
pow(a,2.0)/(4.0 * epsilon0)*rho(i * a, j * a);
}
}
delta = mmax(phi, phip, 101);
swap(phi, phip, 101);
}
}
int main(int argc, char *argv[])
{
clock_t start = clock();
run(1e-6, 0.01);
clock_t end = clock();
double total = ((double)(end - start)) / CLOCKS_PER_SEC;
printf("Execution time: %f\n",total);
}
I decided to write it in Fortran to see how it does performance wise. And the Fortran code is as follows:
module rhofunc
implicit none
public
integer, parameter :: dp=kind(0.d0)
contains
pure real(dp) function rho(x,y)
real(dp), intent(in) :: x,y
if (x > 0.6_dp .and. x < 0.8_dp .and. y > 0.6_dp .and. y<0.8) then
rho = 1.0_dp
else if (x> 0.2_dp .and. x<0.4_dp .and. y>0.2_dp .and. y<0.4_dp) then
rho = -1.0_dp
else
rho = 0.0_dp
end if
end function
end module
program poisson
use rhofunc, only: rho
implicit none
integer, parameter :: dp=kind(0.d0), M=100
integer :: i,j
real(dp),parameter :: epsilon0=8.85e-12_dp, target=1e-6, a=0.01
real(dp) :: delta, start, end, phiprime(M,M), phi(M,M), temp(M,M), a2
delta = 1.0_dp
call cpu_time(start)
do while (delta > target )
a2 = a**2.0_dp
do i=2, M-1
do j=2, M-1
phiprime(i,j) = (phi(i+1,j) + phi(i-1,j) + phi(i,j+1) + phi(i,j-1))/4.0_dp &
+ a2/4.0_dp/epsilon0*rho(i*a,j*a)
end do
end do
delta = maxval(abs(phiprime - phi))
temp = phi
phi = phiprime
phiprime = temp
end do
call cpu_time(end)
print *, end-start
end program
On my machine, the C code runs in 0.12 seconds, but the Fortran code takes 0.37 seconds. Is there a particular reason why the C code is faster here? (I’m able to somewhat improve C performance by doing some unholy pointer tricks when swapping to 0.09 seconds)
Both are single-file, so I simply used gcc/gfortran with the -Ofast optimization flag. Edit: I also ran both with WSL2 on Windows.
After fixing this (plus some more modifications like declaring arrays as 