heat eq, implicit scheme using fort lib

Author: truong-dang-gft

src/fdm_heat_eq_implitcit_1D_use_fort_lib.py

@@ -3,35 +3,36 @@
 import numpy
 
 from src.utils import visu1D
-from ctypes import POINTER, c_int, c_float, cdll
+from ctypes import POINTER, c_int, c_double, cdll
 
 """
 MacOS: Before running this script, please build fortran library [1]:
-    gfortran -dynamiclib src/utils/tridiag.f90 -o src/tridiag.dylib
+    src$ gfortran -shared -fPIC libs/tridiag.f90 -o tridiag.dylib
 There will be a file called 'tridiag.dylib' in the project src folder.
 Reference:
 [1] http://jean-pierre.moreau.pagesperso-orange.fr/Fortran/tridiag_f90.txt
 [2] https://stackoverflow.com/questions/19263879/speeding-up-element-wise-array-multiplication-in-python/19458585#19458585
 """
 ROOT_DIR = os.path.dirname(os.path.abspath(__file__))
 FORTRAN_LIB_PATH = os.path.join(ROOT_DIR, 'tridiag.dylib')
-fortran = cdll.LoadLibrary(FORTRAN_LIB_PATH)
-fortran.tridiag.argtypes = [POINTER(c_float),
-                            POINTER(c_float),
-                            POINTER(c_float),
-                            POINTER(c_float),
-                            POINTER(c_float),
-                            POINTER(c_int),
-                            POINTER(c_int)]
-fortran.tridiag.restype = None
+fort_lib = cdll.LoadLibrary(FORTRAN_LIB_PATH)
+tri_diag = fort_lib.__getattr__("tridiag")
+tri_diag.argtypes = [POINTER(c_double),
+                     POINTER(c_double),
+                     POINTER(c_double),
+                     POINTER(c_double),
+                     POINTER(c_double),
+                     POINTER(c_int),
+                     POINTER(c_int)]
+tri_diag.restype = None
 
 if __name__ == "__main__":
     # FDM, implicit scheme for the heat equation:
     # INPUT
     # Advection velocity:
     mu = 1 / 16
     # Grid size; Use periodic boundary conditions
-    X, M, Tend = 1, 10, 0.5
+    X, M, Tend = 1, 5000, 1
     h = 1 / (M + 1)
     Dt = 0.02
     x = numpy.linspace(0, X, M + 2)  # x(1) =0, x(2) = h, , ..., x(M) = X -h; x(M+1) = X;
@@ -51,15 +52,17 @@
     B = numpy.zeros(M + 2)
     B[:] = 1.0 + 2.0 * mu
 
+    code = c_int(-1)
+    n = c_int(M + 2)
+
     while time < Tend:
         a, b, c = A, B, C
-        code = -1
-        fortran.tridiag(a.ctypes.data_as(POINTER(c_float)),
-                        b.ctypes.data_as(POINTER(c_float)),
-                        b.ctypes.data_as(POINTER(c_float)),
-                        wm.ctypes.data_as(POINTER(c_float)),
-                        wn.ctypes.data_as(POINTER(c_float)),
-                        c_int(M + 2), c_int(code))
+        tri_diag(c.ctypes.data_as(POINTER(c_double)),
+                 b.ctypes.data_as(POINTER(c_double)),
+                 a.ctypes.data_as(POINTER(c_double)),
+                 wm.ctypes.data_as(POINTER(c_double)),
+                 wn.ctypes.data_as(POINTER(c_double)),
+                 n, code)
         wn = wm
 
         time += Dt

src/libs/tridiag.f90

@@ -4,11 +4,10 @@ SUBROUTINE TRIDIAG(A, B, C, R, U, N, CODE) BIND(C, NAME='tridiag')
 ! M U = R, where A, B and C are the three main diagonals of matrix
 ! M(N,N), the other terms are 0. R is the right side vector.
 !*****************************************************************
-    USE ISO_C_BINDING, ONLY: C_FLOAT, C_INT
+    USE ISO_C_BINDING, ONLY: C_DOUBLE, C_INT
     IMPLICIT NONE
-    INTEGER, PARAMETER :: NMAX = 100
-    REAL(C_FLOAT), PARAMETER :: EPS = 1.0D-9
-    REAL(C_FLOAT) :: BET, GAM(NMAX), A(N), B(N), C(N), R(N), U(N)
+    REAL(C_DOUBLE), PARAMETER :: EPS = 1.0D-9
+    REAL(C_DOUBLE) :: BET, GAM(N+1), A(N), B(N), C(N), R(N), U(N)
     INTEGER(C_INT) :: J, N, CODE
 
     IF (B(1) <= EPS) THEN