broken loops in 2d but working jitback for 2d ed bueler example

Author: ZoeLeibowitz

devito/petsc/iet/builder.py

@@ -361,7 +361,7 @@ def _create_dmda_calls(self, dmda):
         )
 
         set_point_bcs = petsc_call(
-            self.callback_builder._point_bc_efunc.name, [dmda, Byref(sobjs['numBC'])]
+            self.callback_builder._point_bc_efunc.name, [dmda, sobjs['numBC']]
         )
 
         get_local_section = petsc_call('DMGetLocalSection', [dmda, Byref(sobjs['lsection'])])

devito/petsc/types/equation.py

@@ -59,12 +59,23 @@ class ConstrainBC(EssentialBC):
 #         return Inc.__new__(Inc, *args, **kwargs)
 
 
-class NoOfEssentialBC(ConstrainBC, Inc):
+class NoOfEssentialBC(ConstrainBC):
     """Equation used count essential boundary condition nodes.
     This type of equation is generated inside
     petscsolve if the user sets `constrain_bcs=True`."""
-    def __new__(cls, *args, **kwargs):
-        return Inc.__new__(Inc, *args, **kwargs)
+    # def __new__(cls, *args, **kwargs):
+    #     return Inc.__new__(Inc, *args, **kwargs)
+    pass
+    
+
+# class NoOfEssentialBC(Inc, ConstrainBC):
+#     """
+#     Equation used to count essential boundary condition nodes.
+#     """
+
+#     def __new__(cls, *args, **kwargs):
+#         obj = super().__new__(cls, *args, **kwargs)
+#         return obj
 
 
 class PointEssentialBC(ConstrainBC):

devito/petsc/types/metadata.py

@@ -760,11 +760,13 @@ def _make_increment_expr(self, expr):
         Make the Eq that is used to increment the number of essential
         boundary nodes in the generated ccode.
         """
-        # rhssymb = PetscInt(name='rhssymb')
         if isinstance(expr, EssentialBC):
             assert expr.lhs == self.target
+            from math import prod
+
+            rhs = prod(expr.rhs.dimensions)
             return NoOfEssentialBC(
-                TempSymb, expr.rhs,
+                TempSymb, rhs,
                 subdomain=expr.subdomain
             )
         else:

examples/petsc/Poisson/ed_bueler_2d.py

@@ -0,0 +1,134 @@
+import os
+import numpy as np
+
+from devito import (Grid, Function, Eq, Operator, switchconfig,
+                    configuration, SubDomain, norm, mmax)
+
+from devito.petsc import petscsolve, EssentialBC
+from devito.petsc.initialize import PetscInitialize
+
+import matplotlib.pyplot as plt
+
+configuration['compiler'] = 'custom'
+os.environ['CC'] = 'mpicc'
+
+
+# 2D test
+# Solving -u_xx - u_yy = f(x,y)
+# Dirichlet BCs: u(0,y) = 0, u(1,y)=-e^y, u(x,0) = -x, u(x,1)=-xe
+# Manufactured solution: u(x,y) = -xe^(y), with corresponding RHS f(x,y) = xe^(y)
+# ref - https://github.com/bueler/p4pdes/blob/master/c/ch6/fish.c
+
+PetscInitialize()
+
+# Subdomains to implement BCs
+class SubTop(SubDomain):
+    name = 'subtop'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: ('middle', 1, 1), y: ('right', 1)}
+
+
+class SubBottom(SubDomain):
+    name = 'subbottom'
+
+    def define(self, dimensions):
+        x, y = dimensions
+        return {x: ('middle', 1, 1), 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 exact(x, y):
+    return -x*np.float64(np.exp(y))
+
+Lx = np.float64(1.)
+Ly = np.float64(1.)
+
+n = 17
+h = Lx/(n-1)
+
+
+grid = Grid(
+    shape=(n, n), extent=(Lx, Ly), subdomains=subdomains, dtype=np.float64
+)
+
+u = Function(name='u', grid=grid, space_order=2)
+f = Function(name='f', grid=grid, space_order=2)
+bc = Function(name='bc', grid=grid, space_order=2)
+
+eqn = Eq(-u.laplace, f, 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)
+
+f.data[:] = X*np.float64(np.exp(Y))
+
+bc.data[0, :] = 0.
+bc.data[-1, :] = -np.exp(tmpy)
+bc.data[:, 0] = -tmpx
+bc.data[:, -1] = -tmpx*np.exp(1)
+
+# # Create boundary condition expressions using subdomains
+bcs = [EssentialBC(u, bc, subdomain=sub1)]
+bcs += [EssentialBC(u, bc, subdomain=sub2)]
+bcs += [EssentialBC(u, bc, subdomain=sub3)]
+bcs += [EssentialBC(u, bc, subdomain=sub4)]
+
+exprs = [eqn] + bcs
+petsc = petscsolve(
+    exprs, target=u,
+    solver_parameters={'ksp_rtol': 1e-12, 'ksp_type': 'cg', 'pc_type': 'none'},
+    options_prefix='poisson_2d',
+    constrain_bcs=True
+)
+
+with switchconfig(log_level='DEBUG'):
+    op = Operator(petsc, language='petsc')
+    summary = op.apply()
+    # print(op.arguments())
+
+
+# print(op.ccode)
+# iters = summary.petsc[('section0', 'poisson_2d')].KSPGetIterationNumber
+
+u_exact = Function(name='u_exact', grid=grid, space_order=2)
+u_exact.data[:] = exact(X, Y)
+print(u_exact)
+
+diff = Function(name='diff', grid=grid, space_order=2)
+diff.data[:] = u_exact.data[:] - u.data[:]
+
+# # Compute infinity norm using numpy
+# # TODO: Figure out how to compute the infinity norm using Devito
+infinity_norm = np.linalg.norm(diff.data[:].ravel(), ord=np.inf)
+print(f"Infinity Norm={infinity_norm}")
+
+# # Compute discrete L2 norm (RMS error)
+n_interior = np.prod([s - 1 for s in grid.shape])
+discrete_l2_norm = norm(diff) / np.sqrt(n_interior)
+print(f"Discrete L2 Norm={discrete_l2_norm}")