Merge pull request #76 from devitocodes/jakubbober-advec

Advec completed

Author: rbanerjee20

.github/workflows/jupyter-notebooks.yml

@@ -82,6 +82,9 @@ jobs:
     - name: Diffusion notebooks (3.7)
       run: |
         $RUN_CMD python -m pytest -W ignore::DeprecationWarning --nbval --cov . --cov-config=.coveragerc --cov-report=xml:diffu_coverage.xml $SKIP fdm-devito-notebooks/03_diffu/diffu_rw.ipynb
+    - name: Advection notebook (4)
+      run: |
+        $RUN_CMD python -m pytest -W ignore::DeprecationWarning --nbval --cov . --cov-config=.coveragerc --cov-report=xml:advec_coverage.xml $SKIP fdm-devito-notebooks/04_advec/advec.ipynb
     - name: Upload coverage to Codecov
       uses: codecov/[email protected]
       with:

fdm-devito-notebooks/04_advec/advec.ipynb

@@ -1,138 +1,143 @@
 import numpy as np
 import matplotlib.pyplot as plt
+from devito import Grid, Eq, solve, TimeFunction, Operator
+
 
 def solver_FECS(I, U0, v, L, dt, C, T, user_action=None):
     Nt = int(round(T/float(dt)))
-    t = np.linspace(0, Nt*dt, Nt+1)   # Mesh points in time
+    t = np.linspace(0, Nt*dt, Nt+1)  # Mesh points in time
     dx = v*dt/C
     Nx = int(round(L/dx))
-    x = np.linspace(0, L, Nx+1)       # Mesh points in space
+    x = np.linspace(0, L, Nx+1)      # Mesh points in space
+    
     # Make sure dx and dt are compatible with x and t
-    dx = x[1] - x[0]
-    dt = t[1] - t[0]
+    dx = float(x[1] - x[0])
+    dt = float(t[1] - t[0])
     C = v*dt/dx
 
-    u   = np.zeros(Nx+1)
-    u_n = np.zeros(Nx+1)
+    grid = Grid(shape=(Nx+1,), extent=(L,))
+    t_s=grid.time_dim
 
-    # Set initial condition u(x,0) = I(x)
-    for i in range(0, Nx+1):
-        u_n[i] = I(x[i])
+    u = TimeFunction(name='u', grid=grid, space_order=2, save=Nt+1)
 
+    pde = u.dtr + v*u.dxc
+    
+    stencil = solve(pde, u.forward)
+    eq = Eq(u.forward, stencil)
+    
+    # Set initial condition u(x,0) = I(x)
+    u.data[1, :] = [I(xi) for xi in x]
+    
+    # Insert boundary condition
+    bc = [Eq(u[t_s+1, 0], U0)]
+    
+    op = Operator([eq] + bc)
+    op.apply(dt=dt, x_m=1, x_M=Nx-1)
+    
     if user_action is not None:
-        user_action(u_n, x, t, 0)
-
-    for n in range(0, Nt):
-        # Compute u at inner mesh points
-        for i in range(1, Nx):
-            u[i] = u_n[i] - 0.5*C*(u_n[i+1] - u_n[i-1])
-
-        # Insert boundary condition
-        u[0] = U0
-
-        if user_action is not None:
-            user_action(u, x, t, n+1)
-
-        # Switch variables before next step
-        u_n, u = u, u_n
+        for n in range(0, Nt + 1):
+            user_action(u.data[n], x, t, n)
 
 
 def solver(I, U0, v, L, dt, C, T, user_action=None,
            scheme='FE', periodic_bc=True):
-    Nt = int(round(T/float(dt)))
+    Nt = int(round(T/np.float64(dt)))
     t = np.linspace(0, Nt*dt, Nt+1)   # Mesh points in time
     dx = v*dt/C
     Nx = int(round(L/dx))
     x = np.linspace(0, L, Nx+1)       # Mesh points in space
+    
     # Make sure dx and dt are compatible with x and t
     dx = x[1] - x[0]
     dt = t[1] - t[0]
     C = v*dt/dx
-    print 'dt=%g, dx=%g, Nx=%d, C=%g' % (dt, dx, Nx, C)
+    print('dt=%g, dx=%g, Nx=%d, C=%g' % (dt, dx, Nx, C))
 
-    u   = np.zeros(Nx+1)
-    u_n = np.zeros(Nx+1)
-    u_nm1 = np.zeros(Nx+1)
     integral = np.zeros(Nt+1)
+    
+    grid = Grid(shape=(Nx+1,), extent=(L,), dtype=np.float64)
+
+    t_s=grid.time_dim
+    
+    def u(to=1, so=1):
+        u = TimeFunction(name='u', grid=grid, time_order=to, space_order=so, save=Nt+1)
+        return u
+        
+    if scheme == 'FE':
+        u   = u(so=2)
+        pde = u.dtr + v*u.dxc
+        
+        pbc = [Eq(u[t_s+1, 0], u[t_s, 0] - 0.5*C*(u[t_s, 1] - u[t_s, Nx]))]
+        pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]
+        
+    elif scheme == 'LF':
+        # Use UP scheme for the first timestep
+        u = u(to=2, so=2)
+        pde0 = u.dtr(fd_order=1) + v*u.dxl(fd_order=1)
+        
+        stencil0 = solve(pde0, u.forward)
+        eq0      = Eq(u.forward, stencil0).subs(t_s, 0)
+
+        pbc0 = [Eq(u[t_s, 0], u[t_s, Nx]).subs(t_s, 0)]
+            
+        # Now continue with LF scheme
+        pde = u.dtc + v*u.dxc
+        
+        pbc = [Eq(u[t_s+1, 0], u[t_s-1, 0] - C*(u[t_s, 1] - u[t_s, Nx-1]))]
+        pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]
+    
+    elif scheme == 'UP':
+        u   = u()
+        pde = u.dtr + v*u.dxl
+        
+        pbc = [Eq(u[t_s, 0], u[t_s, Nx])]
+    
+    elif scheme == 'LW':
+        u = u(so=2)
+        pde = u.dtr + v*u.dxc - 0.5*dt*v**2*u.dx2
+        
+        pbc = [Eq(u[t_s+1, 0], u[t_s, 0] - 0.5*C*(u[t_s, 1] - u[t_s, Nx-1]) + \
+                  0.5*C**2*(u[t_s, 1] - 2*u[t_s, 0] + u[t_s, Nx-1]))]
+        pbc += [Eq(u[t_s+1, Nx], u[t_s+1, 0])]
+    
+    else:
+        raise ValueError('scheme="%s" not implemented' % scheme)
 
+    stencil = solve(pde, u.forward)
+    eq = Eq(u.forward, stencil)
+    
+    bc_init = [Eq(u[t_s+1, 0], U0).subs(t_s, 0)]
+    
     # Set initial condition u(x,0) = I(x)
-    for i in range(0, Nx+1):
-        u_n[i] = I(x[i])
-
-    # Insert boundary condition
-    u[0] = U0
-
+    u.data[0, :] = [I(xi) for xi in x]
+        
     # Compute the integral under the curve
-    integral[0] = dx*(0.5*u_n[0] + 0.5*u_n[Nx] + np.sum(u_n[1:-1]))
-
+    integral[0] = dx*(0.5*u.data[0][0] + 0.5*u.data[0][Nx] + np.sum(u.data[0][1:Nx]))
+    
     if user_action is not None:
-        user_action(u_n, x, t, 0)
-
-    for n in range(0, Nt):
-        if scheme == 'FE':
-            if periodic_bc:
-                i = 0
-                u[i] = u_n[i] - 0.5*C*(u_n[i+1] - u_n[Nx])
-                u[Nx] = u[0]
-                #u[i] = u_n[i] - 0.5*C*(u_n[1] - u_n[Nx])
-            for i in range(1, Nx):
-                u[i] = u_n[i] - 0.5*C*(u_n[i+1] - u_n[i-1])
-        elif scheme == 'LF':
-            if n == 0:
-                # Use upwind for first step
-                if periodic_bc:
-                    i = 0
-                    #u[i] = u_n[i] - C*(u_n[i] - u_n[Nx-1])
-                    u_n[i] = u_n[Nx]
-                for i in range(1, Nx+1):
-                    u[i] = u_n[i] - C*(u_n[i] - u_n[i-1])
-            else:
-                if periodic_bc:
-                    i = 0
-                    # Must have this,
-                    u[i] = u_nm1[i] - C*(u_n[i+1] - u_n[Nx-1])
-                    # not this:
-                    #u_n[i] = u_n[Nx]
-                for i in range(1, Nx):
-                    u[i] = u_nm1[i] - C*(u_n[i+1] - u_n[i-1])
-                if periodic_bc:
-                    u[Nx] = u[0]
-        elif scheme == 'UP':
-            if periodic_bc:
-                u_n[0] = u_n[Nx]
-            for i in range(1, Nx+1):
-                u[i] = u_n[i] - C*(u_n[i] - u_n[i-1])
-        elif scheme == 'LW':
-            if periodic_bc:
-                i = 0
-                # Must have this,
-                u[i] = u_n[i] - 0.5*C*(u_n[i+1] - u_n[Nx-1]) + \
-                       0.5*C*(u_n[i+1] - 2*u_n[i] + u_n[Nx-1])
-                # not this:
-                #u_n[i] = u_n[Nx]
-            for i in range(1, Nx):
-                u[i] = u_n[i] - 0.5*C*(u_n[i+1] - u_n[i-1]) + \
-                       0.5*C*(u_n[i+1] - 2*u_n[i] + u_n[i-1])
-            if periodic_bc:
-                u[Nx] = u[0]
-        else:
-            raise ValueError('scheme="%s" not implemented' % scheme)
-
-        if not periodic_bc:
-            # Insert boundary condition
-            u[0] = U0
-
+        user_action(u.data[0], x, t, 0)
+
+    bc  = [Eq(u[t_s+1, 0], U0)]
+    
+    if scheme == 'LF':
+        op = Operator((pbc0 if periodic_bc else []) + [eq0] + (bc_init if not periodic_bc else []) \
+                      + (pbc if periodic_bc else []) + [eq] + (bc if not periodic_bc else []))
+    else:       
+        op = Operator(bc_init + (pbc if periodic_bc else []) + [eq] + (bc if not periodic_bc else []))
+        
+    op.apply(dt=dt, x_m=1, x_M=Nx if scheme == 'UP' else Nx-1)
+
+    for n in range(1, Nt+1):
         # Compute the integral under the curve
-        integral[n+1] = dx*(0.5*u[0] + 0.5*u[Nx] + np.sum(u[1:-1]))
+        integral[n] = dx*(0.5*u.data[n][0] + 0.5*u.data[n][Nx] + np.sum(u.data[n][1:Nx]))
 
         if user_action is not None:
-            user_action(u, x, t, n+1)
+            user_action(u.data[n], x, t, n)
 
-        # Switch variables before next step
-        u_nm1, u_n, u = u_n, u, u_nm1
-        print 'I:', integral[n+1]
+        print('I:', integral[n])
     return integral
 
+
 def run_FECS(case):
     """Special function for the FECS case."""
     if case == 'gaussian':
@@ -159,14 +164,12 @@ def plot(u, x, t, n):
         m = 40
         if n % m != 0:
             return
-        print 't=%g, n=%d, u in [%g, %g] w/%d points' % \
-              (t[n], n, u.min(), u.max(), x.size)
+        print('t=%g, n=%d, u in [%g, %g] w/%d points' % \
+              (t[n], n, u.min(), u.max(), x.size))
         if np.abs(u).max() > 3:  # Instability?
             return
         plt.plot(x, u)
         legends.append('t=%g' % t[n])
-        if n > 0:
-            plt.hold('on')
 
     plt.ion()
     U0 = 0
@@ -180,6 +183,7 @@ def plot(u, x, t, n):
     plt.axis([0, L, -0.75, 1.1])
     plt.show()
 
+
 def run(scheme='UP', case='gaussian', C=1, dt=0.01):
     """General admin routine for explicit and implicit solvers."""
 
@@ -203,7 +207,6 @@ def plot(u, x, t, n):
             lines = plt.plot(x, u)
             plt.axis([x[0], x[-1], -0.5, 1.5])
             plt.xlabel('x'); plt.ylabel('u')
-            plt.axes().set_aspect(0.15)
             plt.savefig('tmp_%04d.png' % n)
             plt.savefig('tmp_%04d.pdf' % n)
         else:
@@ -219,12 +222,11 @@ def plot(u, x, t, n):
         eps = 1E-14
         if abs(t[n] - 0.6) > eps and abs(t[n] - 0) > eps:
             return
-        print 't=%g, n=%d, u in [%g, %g] w/%d points' % \
-              (t[n], n, u.min(), u.max(), x.size)
+        print('t=%g, n=%d, u in [%g, %g] w/%d points' % \
+              (t[n], n, u.min(), u.max(), x.size))
         if np.abs(u).max() > 3:  # Instability?
             return
         plt.plot(x, u)
-        plt.hold('on')
         plt.draw()
         if n > 0:
             y = [I(x_-v*t[n]) for x_ in x]
@@ -264,16 +266,17 @@ def plot(u, x, t, n):
     plt.figure(2)
     plt.axis([0, L, -0.5, 1.1])
     plt.xlabel('$x$');  plt.ylabel('$u$')
-    plt.axes().set_aspect(0.5)  # no effect
     plt.savefig('tmp1.png'); plt.savefig('tmp1.pdf')
     plt.show()
     # Make videos from figure(1) animation files
     for codec in codecs:
         cmd = 'ffmpeg -i tmp_%%04d.png -r 25 -vcodec %s movie.%s' % \
               (codecs[codec], codec)
         os.system(cmd)
-    print 'Integral of u:', integral.max(), integral.min()
+    print('Integral of u:', integral.max(), integral.min())
+
 
+# TODO: IMPLEMENT THIS IN DEVITO
 def solver_theta(I, v, L, dt, C, T, theta=0.5, user_action=None, FE=False):
     """
     Full solver for the model problem using the theta-rule
@@ -282,7 +285,7 @@ def solver_theta(I, v, L, dt, C, T, theta=0.5, user_action=None, FE=False):
     Vectorized implementation and sparse (tridiagonal)
     coefficient matrix.
     """
-    import time;  t0 = time.clock()  # for measuring the CPU time
+    import time;  t0 = time.process_time()  # for measuring the CPU time
     Nt = int(round(T/float(dt)))
     t = np.linspace(0, Nt*dt, Nt+1)   # Mesh points in time
     dx = v*dt/C
@@ -292,7 +295,7 @@ def solver_theta(I, v, L, dt, C, T, theta=0.5, user_action=None, FE=False):
     dx = x[1] - x[0]
     dt = t[1] - t[0]
     C = v*dt/dx
-    print 'dt=%g, dx=%g, Nx=%d, C=%g' % (dt, dx, Nx, C)
+    print('dt=%g, dx=%g, Nx=%d, C=%g' % (dt, dx, Nx, C))
 
     u   = np.zeros(Nx+1)
     u_n = np.zeros(Nx+1)
@@ -354,7 +357,7 @@ def solver_theta(I, v, L, dt, C, T, theta=0.5, user_action=None, FE=False):
         # Update u_n before next step
         u_n, u = u, u_n
 
-    t1 = time.clock()
+    t1 = time.process_time()
     return integral
 
 

fdm-devito-notebooks/04_advec/src-advec/advec1D.py

@@ -522,7 +522,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### mathbb{R} exponentials\n",
+    "### Real exponentials\n",
     "\n",
     "Let $u^n = \\exp{(\\omega n\\Delta t)} = e^{\\omega t_n}$."
    ]
@@ -1041,7 +1041,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.3"
+   "version": "3.8.2"
   },
   "latex_envs": {
    "LaTeX_envs_menu_present": true,

fdm-devito-notebooks/A_formulas/formulas.ipynb

@@ -18,8 +18,8 @@
 
 - chapter: Advection equations
   sections:
-  - file: notebooks/advec/advec_placeholder.md
+  - file: notebooks/advec/advec.ipynb
 
 - chapter: Nonlinear problems
   sections:
-  - file: notebooks/nonlin/nonlin_placeholder.md
+  - file: notebooks/nonlin/nonlin_placeholder.md