From 689d7c4ed6311108a72622cc49c031e2f8b15b38 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 08:18:31 +0200
Subject: [PATCH 01/11] Added legacy Jacobi solver in fortran.

---
 poisson2d/legacy.f90 | 97 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 97 insertions(+)
 create mode 100644 poisson2d/legacy.f90

diff --git a/poisson2d/legacy.f90 b/poisson2d/legacy.f90
new file mode 100644
index 0000000..b0ae750
--- /dev/null
+++ b/poisson2d/legacy.f90
@@ -0,0 +1,97 @@
+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_dp) 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)
+
+    !----- Physical parameters -----
+    integer , parameter :: m = 300
+    !! Number of grid points in each direction. 
+    real(dp), parameter :: dx = 0.01_dp !1.0_dp / (m-1)
+    !! Uniform grid size in each direction.
+    real(dp) :: rhs(m, m)
+    !! Right-hand side of the Poisson equation.
+    real(dp) :: phi(m, m)
+    !! Solution of the Poisson equation.
+
+    !----- Jacobi solver -----
+    real(dp), parameter :: tolerance = 1e-6_dp
+    !! Tolerance for the iterative Jacobi solver.
+    real(dp) :: phi_prime(m, m)
+    !! Working array.
+    real(dp) :: residual
+    !! Residual is defined as the infinity-norm of the update, i.e. max | phi - phi' |
+    integer  :: iteration
+
+    !----- Miscellaneous -----
+    real(dp) :: start_time, end_time
+    integer             :: i, j
+    real(dp), parameter :: epsilon0 = 8.85E-12_dp
+
+    ! Tic.
+    call cpu_time(start_time)
+
+    ! Initialize right-hand side.
+    do j = 1, m
+        do i = 1, m
+            rhs(i,j) = rho(i*dx,j*dx)
+        end do
+    end do
+
+    !---------------------------------
+    !-----     JACOBI SOLVER     -----
+    !---------------------------------
+
+    ! Iteration counter.
+    iteration = 0
+    ! Working arrays.
+    phi = 0.0_dp ; phi_prime = 0.0_dp
+    ! Residual (ensuring at least one iteration is performed).
+    residual = 1.0_dp
+    ! Rescale rhs.
+    rhs = dx**2/4/epsilon0 * rhs
+
+    ! Iterative solver.
+    do while (residual > tolerance )
+        ! Update iteration counter.
+        iteration = iteration + 1
+        ! Reset residual.
+        residual = 0.0_dp
+        ! Jacobi iteration.
+        do j = 2, m-1
+            do i = 2, m-1
+                ! Jacobi update.
+                phi_prime(i,j) = (phi(i+1,j)+phi(i-1,j)+phi(i,j+1)+phi(i,j-1))/4 + rhs(i, j)
+                ! On-the-fly residual computation.
+                residual = max(residual, abs(phi_prime(i,j) - phi(i,j)))
+            end do
+        end do
+        ! Updated solution.
+        phi = phi_prime 
+    end do
+
+    write(*, *) "Solution converged to the desired tolerance after", iteration, "iterations."
+    
+    ! Toc.
+    call cpu_time(end_time)
+
+    write(*, *) "Computation required", end_time-start_time, "seconds (wall clock time)."
+
+end program

From 2847d636a06e21971a27161c7506406426d8fc6f Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 08:25:57 +0200
Subject: [PATCH 02/11] Run linter on poisson.py

---
 poisson2d/poisson.py | 29 +++++++++++++++--------------
 1 file changed, 15 insertions(+), 14 deletions(-)

diff --git a/poisson2d/poisson.py b/poisson2d/poisson.py
index e1659c1..6d5f051 100644
--- a/poisson2d/poisson.py
+++ b/poisson2d/poisson.py
@@ -13,14 +13,15 @@
 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)
+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:
+
+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:
+    elif x > 0.2 and x < 0.4 and y > 0.2 and y < 0.4:
         return -1
     else:
         return 0
@@ -29,25 +30,25 @@ def ro(x,y):
 # Main loop
 begin = datetime.datetime.now()
 delta = 1.0
-while delta>target:
+while delta > target:
 
     # Calculate new values of the potential
     a2 = a**2
-    for i in range(1,M):
-        for j in range(1,M):
+    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)																												
+            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
+    phi, phiprime = phiprime, phi
 end = datetime.datetime.now()
 dif = end-begin
 print(dif.total_seconds())
 # Make a plot
-plt.imshow(phi,origin='lower')
+plt.imshow(phi, origin='lower')
 plt.savefig("purepython.png")

From d30462dac2cd1b61bba6f58241cb327e4c7aeaca Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 08:27:18 +0200
Subject: [PATCH 03/11] Run linter on setup.py

---
 poisson2d/setup.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/poisson2d/setup.py b/poisson2d/setup.py
index d732f30..6b0104f 100644
--- a/poisson2d/setup.py
+++ b/poisson2d/setup.py
@@ -10,7 +10,7 @@
 from Cython.Build import cythonize
 import numpy
 setup(
-  name = 'yoo',
-  ext_modules = cythonize("poisson.pyx", annotate=True, language_level=3),
-  include_dirs=[numpy.get_include()],
+    name='yoo',
+    ext_modules=cythonize("poisson.pyx", annotate=True, language_level=3),
+    include_dirs=[numpy.get_include()],
 )

From 96665e69fc13300bb251f23a03b5b3c44ced8b85 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 09:37:32 +0200
Subject: [PATCH 04/11] Removed unefficient fortran implementations.

---
 poisson2d/naive.f90     | 51 --------------------------------------
 poisson2d/optimized.f90 | 55 -----------------------------------------
 2 files changed, 106 deletions(-)
 delete mode 100644 poisson2d/naive.f90
 delete mode 100644 poisson2d/optimized.f90

diff --git a/poisson2d/naive.f90 b/poisson2d/naive.f90
deleted file mode 100644
index ecbc75a..0000000
--- a/poisson2d/naive.f90
+++ /dev/null
@@ -1,51 +0,0 @@
-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_dp) 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_dp
-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
diff --git a/poisson2d/optimized.f90 b/poisson2d/optimized.f90
deleted file mode 100644
index 38dddce..0000000
--- a/poisson2d/optimized.f90
+++ /dev/null
@@ -1,55 +0,0 @@
-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_dp) 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_dp
-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

From dd7eac42fc928b115d58035c161f806bcdc6e451 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 09:45:49 +0200
Subject: [PATCH 05/11] Updated rhs module.

---
 poisson2d/legacy.f90 | 42 +++++++++++++++++++++++++-----------------
 1 file changed, 25 insertions(+), 17 deletions(-)

diff --git a/poisson2d/legacy.f90 b/poisson2d/legacy.f90
index b0ae750..a7d3b41 100644
--- a/poisson2d/legacy.f90
+++ b/poisson2d/legacy.f90
@@ -1,25 +1,33 @@
-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_dp) 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
+module rhs_module
+    implicit none
+    public
+
+    integer, parameter :: precision = 15, range = 307
+    integer, parameter :: dp = selected_real_kind(precision, range)
 
+    interface
+        pure real(dp) module function rho(x, y)
+            implicit none
+            real(dp), intent(in) :: x, y
+        end function
+    end interface
 end module
 
+submodule (rhs_module) rho_definition
+    implicit none
+    contains
+    module procedure rho
+        if (all([x, y] > 0.6_dp .and. [x, y] < 0.8_dp)) then
+            rho = 1.0_dp
+        else
+            rho = merge(-1.0_dp, 0.0_dp, all([x, y]>0.2_dp .and. [x, y] < 0.4_dp))
+        endif
+    end procedure
+end submodule
+
 program poisson
-    use rhofunc, only: rho
+    use rhs_module, only: rho, dp
     implicit none
-    integer, parameter :: dp = kind(0.d0)
 
     !----- Physical parameters -----
     integer , parameter :: m = 300

From a57c8dc8aef3239baf082f14d5023c34df824ac5 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 09:55:00 +0200
Subject: [PATCH 06/11] Remove rhs module.

---
 poisson2d/legacy.f90 | 49 +++++++++++++++-----------------------------
 1 file changed, 17 insertions(+), 32 deletions(-)

diff --git a/poisson2d/legacy.f90 b/poisson2d/legacy.f90
index a7d3b41..7dba9ec 100644
--- a/poisson2d/legacy.f90
+++ b/poisson2d/legacy.f90
@@ -1,38 +1,13 @@
-module rhs_module
+program poisson
     implicit none
-    public
-
+    !----- Sets the double precision kind -----
     integer, parameter :: precision = 15, range = 307
     integer, parameter :: dp = selected_real_kind(precision, range)
-
-    interface
-        pure real(dp) module function rho(x, y)
-            implicit none
-            real(dp), intent(in) :: x, y
-        end function
-    end interface
-end module
-
-submodule (rhs_module) rho_definition
-    implicit none
-    contains
-    module procedure rho
-        if (all([x, y] > 0.6_dp .and. [x, y] < 0.8_dp)) then
-            rho = 1.0_dp
-        else
-            rho = merge(-1.0_dp, 0.0_dp, all([x, y]>0.2_dp .and. [x, y] < 0.4_dp))
-        endif
-    end procedure
-end submodule
-
-program poisson
-    use rhs_module, only: rho, dp
-    implicit none
-
+ 
     !----- Physical parameters -----
     integer , parameter :: m = 300
     !! Number of grid points in each direction. 
-    real(dp), parameter :: dx = 0.01_dp !1.0_dp / (m-1)
+    real(dp), parameter :: dx = 1.0_dp / (m-1)
     !! Uniform grid size in each direction.
     real(dp) :: rhs(m, m)
     !! Right-hand side of the Poisson equation.
@@ -74,7 +49,7 @@ program poisson
     ! Residual (ensuring at least one iteration is performed).
     residual = 1.0_dp
     ! Rescale rhs.
-    rhs = dx**2/4/epsilon0 * rhs
+    rhs = dx**2/(4*epsilon0) * rhs
 
     ! Iterative solver.
     do while (residual > tolerance )
@@ -95,11 +70,21 @@ program poisson
         phi = phi_prime 
     end do
 
-    write(*, *) "Solution converged to the desired tolerance after", iteration, "iterations."
+    write(*, *) "Number of iterations  :", iteration
     
     ! Toc.
     call cpu_time(end_time)
 
-    write(*, *) "Computation required", end_time-start_time, "seconds (wall clock time)."
+    write(*, *) "Wall clock time (sec) :", end_time-start_time
 
+contains
+     pure real(dp) function rho(x, y)
+        implicit none
+        real(dp), intent(in) :: x, y
+        if (all([x, y] > 0.6_dp .and. [x, y] < 0.8_dp)) then
+            rho = 1.0_dp
+        else
+            rho = merge(-1.0_dp, 0.0_dp, all([x, y]>0.2_dp .and. [x, y] < 0.4_dp))
+        endif
+   end function
 end program

From 9e5ff2c037b86490ea1a52c7897aa4cc48e36a69 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 09:57:52 +0200
Subject: [PATCH 07/11] Moved initialization of rhs outside the timed section.

---
 poisson2d/legacy.f90 | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/poisson2d/legacy.f90 b/poisson2d/legacy.f90
index 7dba9ec..71bdc50 100644
--- a/poisson2d/legacy.f90
+++ b/poisson2d/legacy.f90
@@ -28,9 +28,6 @@ program poisson
     integer             :: i, j
     real(dp), parameter :: epsilon0 = 8.85E-12_dp
 
-    ! Tic.
-    call cpu_time(start_time)
-
     ! Initialize right-hand side.
     do j = 1, m
         do i = 1, m
@@ -38,6 +35,9 @@ program poisson
         end do
     end do
 
+    ! Tic.
+    call cpu_time(start_time)
+
     !---------------------------------
     !-----     JACOBI SOLVER     -----
     !---------------------------------
@@ -69,12 +69,11 @@ program poisson
         ! Updated solution.
         phi = phi_prime 
     end do
-
-    write(*, *) "Number of iterations  :", iteration
-    
+   
     ! Toc.
     call cpu_time(end_time)
 
+    write(*, *) "Number of iterations  :", iteration
     write(*, *) "Wall clock time (sec) :", end_time-start_time
 
 contains

From 37c198c1335bc36e4a98133a648886333ad9da13 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 10:07:18 +0200
Subject: [PATCH 08/11] Added documentation.

---
 poisson2d/legacy.f90 | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/poisson2d/legacy.f90 b/poisson2d/legacy.f90
index 71bdc50..4157684 100644
--- a/poisson2d/legacy.f90
+++ b/poisson2d/legacy.f90
@@ -1,4 +1,11 @@
 program poisson
+    !! This program solves the Poisson equation on the unit square with
+    !! homogeneous Dirichlet boundary conditions. The Laplace operator
+    !! is discretized with the standard second-order accurate central
+    !! finite difference scheme. The resulting discrete problem is solved
+    !! using the Jacobi iterative solver. The numerical implementation
+    !! relies on relatively standard Fortran constructs and is very similar
+    !! to what would be asked from students in a numerical analysis class.
     implicit none
     !----- Sets the double precision kind -----
     integer, parameter :: precision = 15, range = 307

From bbc2e64b3e0fe8226d5500cc2381e2df2fce436d Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 10:09:28 +0200
Subject: [PATCH 09/11] Created Jacobi folder if other solvers are implemented
 later on.

---
 poisson2d/{ => jacobi}/README.md   | 0
 poisson2d/{ => jacobi}/legacy.f90  | 0
 poisson2d/{ => jacobi}/naive.c     | 0
 poisson2d/{ => jacobi}/optimized.c | 0
 poisson2d/{ => jacobi}/poisson.py  | 0
 poisson2d/{ => jacobi}/poisson.pyx | 0
 poisson2d/{ => jacobi}/setup.py    | 0
 7 files changed, 0 insertions(+), 0 deletions(-)
 rename poisson2d/{ => jacobi}/README.md (100%)
 rename poisson2d/{ => jacobi}/legacy.f90 (100%)
 rename poisson2d/{ => jacobi}/naive.c (100%)
 rename poisson2d/{ => jacobi}/optimized.c (100%)
 rename poisson2d/{ => jacobi}/poisson.py (100%)
 rename poisson2d/{ => jacobi}/poisson.pyx (100%)
 rename poisson2d/{ => jacobi}/setup.py (100%)

diff --git a/poisson2d/README.md b/poisson2d/jacobi/README.md
similarity index 100%
rename from poisson2d/README.md
rename to poisson2d/jacobi/README.md
diff --git a/poisson2d/legacy.f90 b/poisson2d/jacobi/legacy.f90
similarity index 100%
rename from poisson2d/legacy.f90
rename to poisson2d/jacobi/legacy.f90
diff --git a/poisson2d/naive.c b/poisson2d/jacobi/naive.c
similarity index 100%
rename from poisson2d/naive.c
rename to poisson2d/jacobi/naive.c
diff --git a/poisson2d/optimized.c b/poisson2d/jacobi/optimized.c
similarity index 100%
rename from poisson2d/optimized.c
rename to poisson2d/jacobi/optimized.c
diff --git a/poisson2d/poisson.py b/poisson2d/jacobi/poisson.py
similarity index 100%
rename from poisson2d/poisson.py
rename to poisson2d/jacobi/poisson.py
diff --git a/poisson2d/poisson.pyx b/poisson2d/jacobi/poisson.pyx
similarity index 100%
rename from poisson2d/poisson.pyx
rename to poisson2d/jacobi/poisson.pyx
diff --git a/poisson2d/setup.py b/poisson2d/jacobi/setup.py
similarity index 100%
rename from poisson2d/setup.py
rename to poisson2d/jacobi/setup.py

From dcd69d4524b823ff5536a71db78a576c7440bba3 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 10:12:26 +0200
Subject: [PATCH 10/11] Removed  from the definition of rho.

---
 poisson2d/jacobi/legacy.f90 | 1 -
 1 file changed, 1 deletion(-)

diff --git a/poisson2d/jacobi/legacy.f90 b/poisson2d/jacobi/legacy.f90
index 4157684..2781bb5 100644
--- a/poisson2d/jacobi/legacy.f90
+++ b/poisson2d/jacobi/legacy.f90
@@ -85,7 +85,6 @@ program poisson
 
 contains
      pure real(dp) function rho(x, y)
-        implicit none
         real(dp), intent(in) :: x, y
         if (all([x, y] > 0.6_dp .and. [x, y] < 0.8_dp)) then
             rho = 1.0_dp

From 0efb3570a62eb66083abf73e838fae50c32d7da7 Mon Sep 17 00:00:00 2001
From: Jean-Christophe <loiseau.jc@gmail.com>
Date: Tue, 21 May 2024 15:07:48 +0200
Subject: [PATCH 11/11] Removed trailing whitespace.

---
 poisson2d/jacobi/legacy.f90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/poisson2d/jacobi/legacy.f90 b/poisson2d/jacobi/legacy.f90
index 2781bb5..d12c31d 100644
--- a/poisson2d/jacobi/legacy.f90
+++ b/poisson2d/jacobi/legacy.f90
@@ -59,7 +59,7 @@ program poisson
     rhs = dx**2/(4*epsilon0) * rhs
 
     ! Iterative solver.
-    do while (residual > tolerance )
+    do while (residual > tolerance)
         ! Update iteration counter.
         iteration = iteration + 1
         ! Reset residual.