Add some examples (#76)

Author: ZoeLeibowitz

.github/workflows/pytest-petsc.yml

@@ -80,7 +80,11 @@ jobs:
 
     - name: Test examples
       run: |
-        ${{ env.RUN_CMD }} mpiexec -n 1 python3 --cov --cov-config=.coveragerc --cov-report=xml examples/petsc/seismic/staggered_acoustic.py
+        ${{ env.RUN_CMD }} mpiexec -n 1 python3 examples/petsc/seismic/01_staggered_acoustic.py
+        ${{ env.RUN_CMD }} mpiexec -n 1 python3 examples/petsc/cfd/01_navierstokes.py
+        ${{ env.RUN_CMD }} mpiexec -n 1 python3 examples/petsc/Poisson/01_poisson.py
+        ${{ env.RUN_CMD }} mpiexec -n 1 python3 examples/petsc/Poisson/02_laplace.py
+        ${{ env.RUN_CMD }} mpiexec -n 1 python3 examples/petsc/random/01_helmholtz.py
 
     - name: Upload coverage to Codecov
       if: "!contains(matrix.name, 'docker')"

devito/petsc/types/array.py

@@ -5,7 +5,6 @@
 from devito.types.array import ArrayBasic
 from devito.finite_differences import Differentiable
 from devito.types.basic import AbstractFunction
-from devito.finite_differences.tools import fd_weights_registry
 from devito.tools import dtype_to_ctype, as_tuple
 from devito.symbolics import FieldFromComposite
 
@@ -27,9 +26,6 @@ class PETScArray(ArrayBasic, Differentiable):
 
     _data_alignment = False
 
-    # Default method for the finite difference approximation weights computation.
-    _default_fd = 'taylor'
-
     __rkwargs__ = (AbstractFunction.__rkwargs__ +
                    ('target', 'liveness', 'coefficients', 'localinfo'))
 
@@ -38,15 +34,8 @@ def __init_finalize__(self, *args, **kwargs):
         self._target = kwargs.get('target')
         self._ndim = kwargs['ndim'] = len(self._target.space_dimensions)
         self._dimensions = kwargs['dimensions'] = self._target.space_dimensions
-
         super().__init_finalize__(*args, **kwargs)
-
-        # Symbolic (finite difference) coefficients
-        self._coefficients = kwargs.get('coefficients', self._default_fd)
-        if self._coefficients not in fd_weights_registry:
-            raise ValueError("coefficients must be one of %s"
-                             " not %s" % (str(fd_weights_registry), self._coefficients))
-
+        self._coefficients = self._target.coefficients
         self._localinfo = kwargs.get('localinfo', None)
 
     @property
@@ -93,6 +82,10 @@ def shape(self):
     def space_order(self):
         return self.target.space_order
 
+    @property
+    def staggered(self):
+        return self.target.staggered
+
     @property
     def is_Staggered(self):
         return self.target.staggered is not None
@@ -113,6 +106,10 @@ def shape_allocated(self):
     def _C_ctype(self):
         return POINTER(dtype_to_ctype(self.dtype))
 
+    @cached_property
+    def _fd_priority(self):
+        return self.target._fd_priority
+
     @property
     def symbolic_shape(self):
         field_from_composites = [

examples/petsc/Poisson/01_poisson.py

@@ -0,0 +1,112 @@
+import os
+import numpy as np
+
+from devito import (Grid, Function, Eq, Operator, switchconfig,
+                    configuration, SubDomain)
+
+from devito.petsc import PETScSolve, EssentialBC
+from devito.petsc.initialize import PetscInitialize
+configuration['compiler'] = 'custom'
+os.environ['CC'] = 'mpicc'
+
+# Solving pn.laplace = 2x(y - 1)(y - 2x + xy + 2)e^(x-y)
+# Constant zero Dirichlet BCs.
+
+PetscInitialize()
+
+
+# Subdomains to implement BCs
+class SubTop(SubDomain):
+    name = 'subtop'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: x, y: ('right', 1)}
+
+
+class SubBottom(SubDomain):
+    name = 'subbottom'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: x, y: ('left', 1)}
+
+
+class SubLeft(SubDomain):
+    name = 'subleft'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: ('left', 1), y: y}
+
+
+class SubRight(SubDomain):
+    name = 'subright'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: ('right', 1), y: y}
+
+
+sub1 = SubTop()
+sub2 = SubBottom()
+sub3 = SubLeft()
+sub4 = SubRight()
+
+subdomains = (sub1, sub2, sub3, sub4)
+
+
+def analytical(x, y):
+    return np.float64(np.exp(x-y) * x * (1-x) * y * (1-y))
+
+
+Lx = np.float64(1.)
+Ly = np.float64(1.)
+
+n_values = list(range(13, 174, 10))
+dx = np.array([Lx/(n-1) for n in n_values])
+errors = []
+
+
+for n in n_values:
+
+    grid = Grid(
+        shape=(n, n), extent=(Lx, Ly), subdomains=subdomains, dtype=np.float64
+    )
+
+    phi = Function(name='phi', grid=grid, space_order=2, dtype=np.float64)
+    rhs = Function(name='rhs', grid=grid, space_order=2, dtype=np.float64)
+
+    eqn = Eq(rhs, phi.laplace, subdomain=grid.interior)
+
+    tmpx = np.linspace(0, Lx, n).astype(np.float64)
+    tmpy = np.linspace(0, Ly, n).astype(np.float64)
+    Y, X = np.meshgrid(tmpx, tmpy)
+
+    rhs.data[:] = np.float64(
+        2.0*X*(Y-1.0)*(Y - 2.0*X + X*Y + 2.0)
+    ) * np.float64(np.exp(X-Y))
+
+    # # Create boundary condition expressions using subdomains
+    bcs = [EssentialBC(phi, np.float64(0.), subdomain=sub1)]
+    bcs += [EssentialBC(phi, np.float64(0.), subdomain=sub2)]
+    bcs += [EssentialBC(phi, np.float64(0.), subdomain=sub3)]
+    bcs += [EssentialBC(phi, np.float64(0.), subdomain=sub4)]
+
+    exprs = [eqn] + bcs
+    petsc = PETScSolve(exprs, target=phi, solver_parameters={'ksp_rtol': 1e-8})
+
+    with switchconfig(language='petsc'):
+        op = Operator(petsc)
+        op.apply()
+
+    phi_analytical = analytical(X, Y)
+
+    diff = phi_analytical[1:-1, 1:-1] - phi.data[1:-1, 1:-1]
+    error = np.linalg.norm(diff) / np.linalg.norm(phi_analytical[1:-1, 1:-1])
+    errors.append(error)
+
+slope, _ = np.polyfit(np.log(dx), np.log(errors), 1)
+
+assert slope > 1.9
+assert slope < 2.1

examples/petsc/Poisson/02_laplace.py

@@ -0,0 +1,125 @@
+import os
+import numpy as np
+
+from devito import (Grid, Function, Eq, Operator, SubDomain,
+                    configuration, switchconfig)
+from devito.petsc import PETScSolve, EssentialBC
+from devito.petsc.initialize import PetscInitialize
+configuration['compiler'] = 'custom'
+os.environ['CC'] = 'mpicc'
+
+PetscInitialize()
+
+# Laplace equation, solving phi.laplace = 0
+
+# Constant Dirichlet BCs:
+# phi(x, 0) = 0
+# phi(0, y) = 0
+# phi(1, y) = 0
+# phi(x, 1) = f(x) = sin(pi*x)
+
+# The analytical solution is:
+# phi(x, y) = sinh(pi*y)*sin(pi*x)/sinh(pi)
+
+
+# Subdomains to implement BCs
+class SubTop(SubDomain):
+    name = 'subtop'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: x, y: ('right', 1)}
+
+
+class SubBottom(SubDomain):
+    name = 'subbottom'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: x, y: ('left', 1)}
+
+
+class SubLeft(SubDomain):
+    name = 'subleft'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: ('left', 1), y: y}
+
+
+class SubRight(SubDomain):
+    name = 'subright'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: ('right', 1), y: y}
+
+
+sub1 = SubTop()
+sub2 = SubBottom()
+sub3 = SubLeft()
+sub4 = SubRight()
+
+subdomains = (sub1, sub2, sub3, sub4)
+
+
+Lx = np.float64(1.)
+Ly = np.float64(1.)
+
+
+def analytical(x, y, Lx, Ly):
+    tmp = np.float64(np.pi)/Lx
+    numerator = np.float64(np.sinh(tmp*y)) * np.float64(np.sin(tmp*x))
+    return numerator / np.float64(np.sinh(tmp*Ly))
+
+
+n_values = list(range(13, 174, 10))
+dx = np.array([Lx/(n-1) for n in n_values])
+errors = []
+
+
+for n in n_values:
+
+    grid = Grid(
+        shape=(n, n), extent=(Lx, Ly), subdomains=subdomains, dtype=np.float64
+    )
+
+    phi = Function(name='phi', grid=grid, space_order=2, dtype=np.float64)
+    rhs = Function(name='rhs', grid=grid, space_order=2, dtype=np.float64)
+
+    phi.data[:] = np.float64(0.0)
+    rhs.data[:] = np.float64(0.0)
+
+    eqn = Eq(rhs, phi.laplace, subdomain=grid.interior)
+
+    tmpx = np.linspace(0, Lx, n).astype(np.float64)
+    tmpy = np.linspace(0, Ly, n).astype(np.float64)
+    Y, X = np.meshgrid(tmpx, tmpy)
+
+    # Create boundary condition expressions using subdomains
+    bc_func = Function(name='bcs', grid=grid, space_order=2, dtype=np.float64)
+    bc_func.data[:] = np.float64(0.0)
+    bc_func.data[:, -1] = np.float64(np.sin(tmpx*np.pi))
+
+    bcs = [EssentialBC(phi, bc_func, subdomain=sub1)]  # top
+    bcs += [EssentialBC(phi, bc_func, subdomain=sub2)]  # bottom
+    bcs += [EssentialBC(phi, bc_func, subdomain=sub3)]  # left
+    bcs += [EssentialBC(phi, bc_func, subdomain=sub4)]  # right
+
+    exprs = [eqn] + bcs
+    petsc = PETScSolve(exprs, target=phi, solver_parameters={'ksp_rtol': 1e-8})
+
+    with switchconfig(language='petsc'):
+        op = Operator(petsc)
+        op.apply()
+
+    phi_analytical = analytical(X, Y, Lx, Ly)
+
+    diff = phi_analytical[1:-1, 1:-1] - phi.data[1:-1, 1:-1]
+    error = np.linalg.norm(diff) / np.linalg.norm(phi_analytical[1:-1, 1:-1])
+    errors.append(error)
+
+slope, _ = np.polyfit(np.log(dx), np.log(errors), 1)
+
+assert slope > 1.9
+assert slope < 2.1