Merge pull request #9 from arunningcroc/master

Poisson 2D Jacobi solver benchmark with Python, Cython, C, and Fortran code

Author: milancurcic

poisson2d/README.md

@@ -0,0 +1,31 @@
+## Poisson solver in various languages
+
+This solver uses a simple Jacobi iteration to solve a Poisson equation. For details, see, for example,
+"Computational Physics" by Mark Newman, Chap 9. Compiler commands were:
+
+```gfortran sourcefile.f95 -Ofast -o program```
+
+```gcc sourcefile.c -Ofast -o program```
+
+For Cython module (import it to run):
+
+```python setup.py build_ext --inplace```
+
+The grid is a 2-dimensional array of size MxM. The timings for different values of M are, in seconds:
+
+| Language            | M=100           | M=200                 | M=300                  |
+|---------------------|-----------------|-----------------------|------------------------|
+| Python (pure)       | 276             | n/a                   | n/a                    |
+| Cython              | 1.02            | 32.8                  | 229                    |
+| Fortran (naive)     | 0.34            | 13.2                  | 69.7                   |
+| Fortran (optimized) | 0.18            | 6.25                  | 31.4                   |
+| C (naive)           | 0.42*           | 7.25                  | 33.7                   |
+| C (optimized)       | 0.37*           | 6.80                  | 32.8                   |
+
+* For all of these results, the amount of iterations performed by the respective codes was approximately
+the same, with the exception of the 100x100 grid C codes, which did nearly a double amount of iterations
+compared to the rest, for some reason. The timings are on AMD Ryzen 5 3600 @3.6GHz, using WSL2 on Windows 10.
+
+Some thoughts on the code at https://runningcrocodile.fi/articles/pythonperformance.html . Also, there was
+discussion about this problem at Discourse: https://fortran-lang.discourse.group/t/performance-c-vs-fortran/1461 
+with important contributions from numerous users (Mason, Beliavsky, septc, implicitall, nncarlson, han190 and pmk)

poisson2d/naive.c

@@ -0,0 +1,90 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+#define M 300
+void swap(double (*a)[M], double (*b)[M], int n)
+{
+  double swp;
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      swp = a[i][j];
+      a[i][j] = b[i][j];
+      b[i][j] = swp;
+    }
+  }
+}
+double mmax(double (*phi)[M], double (*phip)[M], 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 = fabs(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 a2;
+
+  double (*phi)[M];
+  double (*phip)[M]; 
+  double (*rhoa)[M]; 
+  phi = malloc(sizeof(double[M][M]));
+  phip = malloc(sizeof(double[M][M]));
+  rhoa = malloc(sizeof(double[M][M]));
+
+  int iter = 0;
+
+  memset(phip, 0, sizeof(phip[0][0])*M*M);
+  memset(phi, 0, sizeof(phi[0][0])*M*M);
+  
+  double delta = 1.0;
+  a2 = pow(a,2.0);
+  while (delta > toler) {
+    iter += 1;
+    
+    for (int i=1; i < M-1; i++) {
+      for (int j=1; j < M-1; j++) {
+        phip[i][j] = (phi[i+1][j] + phi[i-1][j] +
+                      phi[i][j+1] + phi[i][j-1])/4.0 +
+                      a2/(4.0 * epsilon0)*rho(i*a,j*a);
+      }
+    }
+    delta = mmax(phi, phip, M);
+    swap(phi, phip, M);
+
+  }
+  printf("iters %d", iter);
+
+}
+
+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);
+}

poisson2d/naive.f90

@@ -0,0 +1,51 @@
+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=300
+integer            :: i,j, iter
+real(dp),parameter :: epsilon0=8.85E-12_dp, target=1E-6_dp, a=0.01
+real(dp)           :: delta, b, e, phiprime(M,M), phi(M,M), a2, rhoarr(M,M), temp(M,M)
+
+
+delta = 1.0_dp
+iter = 0
+call cpu_time(b)
+phiprime(:,:) = 0.0_dp
+phi(:,:) = 0.0_dp
+
+do while (delta > target )
+    iter = iter + 1
+    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(e)
+print *, e-b, iter
+end program

poisson2d/optimized.c

@@ -0,0 +1,92 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <time.h>
+#define M 300
+void swap(double (*a)[M], double (*b)[M], int n)
+{
+  double swp;
+  for (int i = 0; i < n; i++) {
+    for (int j = 0; j < n; j++) {
+      a[i][j] = b[i][j];
+    }
+  }
+}
+double mmax(double (*phi)[M], double (*phip)[M], 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 = fabs(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 a2;
+
+  double (*phi)[M];
+  double (*phip)[M]; 
+  double (*rhoa)[M]; 
+  phi = malloc(sizeof(double[M][M]));
+  phip = malloc(sizeof(double[M][M]));
+  rhoa = malloc(sizeof(double[M][M]));
+
+  int iter = 0;
+
+  memset(phip, 0, sizeof(phip[0][0])*M*M);
+  memset(phi, 0, sizeof(phi[0][0])*M*M);
+  for (int i=1; i<M-1; i++)  {
+    for (int j=1; j<M-1; j++) {
+      rhoa[i][j] =  rho(i*a,j*a);
+    }
+  }
+  double delta = 1.0;
+  a2 = pow(a,2.0);
+  while (delta > toler) {
+    iter += 1;
+    
+    for (int i=1; i < M-1; i++) {
+      for (int j=1; j < M-1; j++) {
+        phip[i][j] = (phi[i+1][j] + phi[i-1][j] +
+                      phi[i][j+1] + phi[i][j-1])/4.0 +
+                      a2/(4.0 * epsilon0)*rhoa[i][j];
+      }
+    }
+    delta = mmax(phi, phip, M);
+    swap(phi, phip, M);
+
+  }
+  printf("iters %d", iter);
+
+}
+
+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);
+}

poisson2d/optimized.f90

@@ -0,0 +1,55 @@
+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=300
+integer            :: i,j, iter
+real(dp),parameter :: epsilon0=8.85E-12_dp, target=1E-6_dp, a=0.01
+real(dp)           :: delta, b, e, phiprime(M,M), phi(M,M), a2, rhoarr(M,M)
+
+
+delta = 1.0_dp
+iter = 0
+call cpu_time(b)
+phiprime(:,:) = 0.0_dp
+phi(:,:) = 0.0_dp
+do i=1, M
+    do j=1, M
+        rhoarr(i,j) = rho(i*a,j*a)
+    end do
+end do
+
+do while (delta > target )
+    iter = iter + 1
+    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*rhoarr(i,j)
+        end do
+    end do
+    delta = maxval(abs(phiprime - phi))
+    phi = phiprime 
+
+
+end do
+call cpu_time(e)
+print *, e-b, iter
+end program

poisson2d/poisson.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import matplotlib.pyplot as plt
+import datetime
+
+# Constants
+M = 100         # Grid squares on a side
+V = 1.0         # Voltage at top wall
+target = 1e-6   # Target accuracy
+a = 0.01
+epsilon0 = 8.85e-12
+
+# Create arrays to hold potential values
+phi = np.zeros([M+1,M+1],float)
+#phi[0,:] = V
+phiprime = np.zeros([M+1,M+1],float)
+
+def ro(x,y):
+    if x>0.6 and x<0.8 and y>0.6 and y<0.8:
+        return 1
+    elif x>0.2 and x<0.4 and y>0.2 and y<0.4:
+        return -1
+    else:
+        return 0
+
+
+# Main loop
+begin = datetime.datetime.now()
+delta = 1.0
+while delta>target:
+
+    # Calculate new values of the potential
+    a2 = a**2
+    for i in range(1,M):
+        for j in range(1,M):
+
+            phiprime[i,j] = (phi[i+1,j] + phi[i-1,j] \
+                             + phi[i,j+1] + phi[i,j-1])/4 \
+                    + a2/4/epsilon0*ro(i*a,j*a)																												
+
+    # Calculate maximum difference from old values
+    delta = np.max(abs(phi-phiprime))
+
+    # Swap the two arrays around
+    phi,phiprime = phiprime,phi
+end = datetime.datetime.now()
+dif = end-begin
+print(dif.total_seconds())
+# Make a plot
+plt.imshow(phi,origin='lower')
+plt.savefig("purepython.png")