Add a couple equations, test parallelization

OpenMP is still slow for a problem this size.

src/nonlinalg.f90

@@ -3,7 +3,7 @@ module nonlinalg
   use linalg, only: linsolve_quick
   implicit none
 
-  real(wp), parameter :: pi = 4._wp * datan(1._wp)
+  ! real(wp), parameter :: pi = 4._wp * datan(1._wp)
 
   interface
     function fun_interf(n, x) result(y)
@@ -52,7 +52,7 @@ function jacobian(fun, n, x) result(dydx)
     procedure(fun_interf)                 :: fun
 
     integer                   :: jj
-    real(wp), dimension(n)    :: xx
+    real(wp), dimension(n)    :: xx, fun_plus, fun_minus
     real(wp), parameter       :: eps = epsilon(1e0)
 
     ! print*,
@@ -66,15 +66,18 @@ function jacobian(fun, n, x) result(dydx)
     ! end do
     ! stop
 
-    !$omp parallel do private (jj, xx) shared(dydx, x, n)
+    !$omp parallel do &
+    !$omp& private (jj, xx, fun_plus, fun_minus) shared(dydx, x, n)
     do jj = 1,n
       xx = x
       xx(jj) = xx(jj) + eps
-      dydx(:,jj) = fun(n, xx)
+      fun_plus = fun(n, xx)
 
       xx = x
       xx(jj) = xx(jj) - eps
-      dydx(:,jj) = ( dydx(:,jj) - fun(n, xx) ) / (2._wp * eps)
+      fun_minus = fun(n, xx)
+
+      dydx(:,jj) = (fun_plus - fun_minus ) / (2._wp * eps)
     end do
     !$omp end parallel do
 

src/runge_kutta.f90

@@ -77,7 +77,7 @@ subroutine rk_wrapper(n, num_t, t, y, dysub)
     ! EDIT 10-11-2016: change dt to be the minimum value between t(2)-t(1) and
     ! 1d-2 because sometimes you only want a very sparse output of data, and
     ! that could overshoot a good initial guess. This should be smarter
-    dt = minval([t(2) - t(1), 1d-1])
+    dt = minval([t(2) - t(1), 1d-2])
 
     do ii = 1, num_t-1
       yy = y(ii,:)
@@ -112,8 +112,12 @@ subroutine rk_adaptive(n, tspan, dt, y, dysub, a, b, bstar, c, m, p, explicit)
 
     ! Local arguments
     real(wp)                              :: t, error
-    real(wp), parameter                   :: eps_abs = 1d-3 ! sqrt(epsilon(1._wp))
-    real(wp), parameter                   :: eps_rel = 1d-3
+    ! real(wp), parameter                   :: eps_abs = epsilon(1e0)
+    ! real(wp), parameter                   :: eps_abs = 1d-7
+    real(wp), parameter                   :: eps_abs = sqrt(epsilon(1._wp))
+    ! real(wp), parameter                   :: eps_rel = epsilon(1e0)
+    ! real(wp), parameter                   :: eps_rel = 1d-7
+    real(wp), parameter                   :: eps_rel = sqrt(epsilon(1._wp))
     real(wp), dimension(n)                :: dummy_y, y_star
 
     integer                               :: eval, num_dt

src/test/plotter.py

@@ -12,25 +12,27 @@
 
 df = pd.read_csv(data_file,
                  delim_whitespace=True,
-                 names=['t', 'y1', 'y2'])
+                 names=['t', 'y1', 'y2', 'y3'])
+                #  names=['t', 'y1', 'y2'])
                 #  names=['t', 'x', 'y', 'z'])
                 #  names=['t', 'x1', 'x2', 'x3'])
 
 df_raw = pd.read_csv(raw_file,
                      delim_whitespace=True,
-                     names=['t', 'dt', 'y1', 'y2'])
+                     names=['t', 'dt', 'y1', 'y2', 'y3'])
 
 df.sort_values(by='t', inplace=True)
 df_raw.sort_values(by='t', inplace=True)
 
 t   = df_raw.t.values
 y1  = df_raw.y1.values
 y2  = df_raw.y2.values
-# y3  = df.y3.values
+y3  = df_raw.y3.values
 
 ax1 = fig.add_subplot(211)
 ax1.plot(t, y1, marker='o', linestyle='-', label='y1')
 ax1.plot(t, y2, marker='o', linestyle='-', label='y2')
+ax1.plot(t, y3, marker='o', linestyle='-', label='y3')
 # ax1.plot(t, x3, marker='o', linestyle='None', label='x3')
 
 ax1.set_xlabel('t')

src/test/src/main.f90

@@ -6,7 +6,7 @@ program main
   implicit none
 
   integer                       :: ii, ios
-  integer, parameter            :: n = 2, num_t = 2
+  integer, parameter            :: n = 3, num_t = 2
   real(wp)                      :: t0, tend
   real(wp), dimension(num_t)    :: t
   real(wp), dimension(n)        :: y0
@@ -16,14 +16,15 @@ program main
   ! Set time vector `t`
   t0   = 0.0_wp
   ! tend    = 4_wp*pi
-  tend = 1000._wp
+  tend = 50._wp
   call linspace(t0, tend, t)
 
   ! Set y vector `y` using y0
   ! y0 = [0._wp, 0.1_wp]
-  y0 = [2._wp, 0._wp]
+  ! y0 = [2._wp, 0._wp]
+  ! y0 = [ 0._wp ]
   ! y0 = [0.1_wp, 10.0_wp, 10.0_wp]
-  ! y0      = [0._wp, 1._wp, 0._wp]
+  y0      = [1._wp, 0._wp, 0._wp]
   y(1,:)  = y0
 
 

src/test/src/misc.f90

@@ -14,10 +14,11 @@ subroutine mysub(n, t, y, dy)
 
     ! call simple_trig(n, t, y, dy)
     ! call simple_ode(n, t, y, dy)
-    call vanderpol(n, t, y, dy)
+    ! call vanderpol(n, t, y, dy)
     ! call lorenz(n, t, y, dy)
     ! call brusselator(n, t, y, dy)
-
+    ! call stiff(n, t, y, dy)
+    call robertson(n, t, y, dy)
 
     return
   end subroutine mysub
@@ -91,4 +92,36 @@ subroutine brusselator(n, t, y, dy)
     return
   end subroutine brusselator
 
+  subroutine stiff(n, t, y, dy)
+    integer, intent(in)                 :: n
+    real(wp), intent(in)                :: t
+    real(wp), intent(in), dimension(n)  :: y
+    real(wp), intent(out), dimension(n) :: dy
+
+    dy(1) = -50._wp * ( y(1) - dcos(t) )
+
+    return
+  end subroutine stiff
+
+  subroutine robertson(n, t, y, dy)
+    integer, intent(in)                 :: n
+    real(wp), intent(in)                :: t
+    real(wp), intent(in), dimension(n)  :: y
+    real(wp), intent(out), dimension(n) :: dy
+
+    real(wp), dimension(n) :: k
+
+    k = [ 4d-2, 3d7, 1d4 ]
+
+    dy(1) = - k(1) * y(1) + k(3) * y(2) * y(3)
+    dy(2) =   k(1) * y(1) - k(3) * y(2) * y(3) - k(2) * y(2) * y(2)
+    dy(3) =                                      k(2) * y(2) * y(2)
+
+    ! print*, y
+    ! print*, dy
+    ! stop
+
+    return
+  end subroutine robertson
+
 end module misc