examples: Add 04_poisson example

Author: ZoeLeibowitz

examples/petsc/Poisson/03_poisson.py

@@ -10,10 +10,6 @@
 os.environ['CC'] = 'mpicc'
 
 # 1D test
-# ref - An efficient implementation of fourth-order compact finite
-# difference scheme for Poisson equation with Dirichlet
-# boundary conditions
-# https://pdf.sciencedirectassets.com/271503/1-s2.0-S0898122116X00090/1-s2.0-S0898122116300761/main.pdf?X-Amz-Security-Token=IQoJb3JpZ2luX2VjENH%2F%2F%2F%2F%2F%2F%2F%2F%2F%2FwEaCXVzLWVhc3QtMSJGMEQCIH%2F3tmUOBoYeSj%2FmE47N2yF5b5IvdQohX7JWYGRZhNFTAiAX0P7PnDVBy2D3pTjsZRFgFEh7nVzV8oxoUFevx%2FucWSqzBQhZEAUaDDA1OTAwMzU0Njg2NSIMHBuo2q1ii2daxqQiKpAFdh12uILywsya5sCFqTsXwUdb%2F4lvp5uZIigTrb18tjM8cPru5xrZDgDVYVIlT4G6L1SE05FkjWKSQ7AO24wec3y2bNKGgpC%2FPFbEmTv6CgpR%2FjoSpToblGLiOqTkUgICSK3EoMdah%2BWl552nO0Ajdwuor0brGfDM7C2fgH1FqM%2BLyJ2do33rYFjGAswzZsQOGcf%2BChdkaKA9bxvfgE2%2Bukf7RcYDsGteSEX5Zb9XoyvMheiMUZoZk7KVPWjj3JORx9qetLs9LkpPO3IU%2BqPxtM7Vt3BnEnXR9gQ2bnL%2FtcT%2FcvsZ7a8AdiU1j8%2F%2Fxi9nBgPow0MQTmaoe9n67XRS0BVE7wAWldDb2qdZuOfwYl%2F2iG78mMTn%2FC4YcOCezc4nUT9fTcTcv3wKZzA%2Bkh8Z%2BXvdTcdADCKdVaIXLylqlhEmBlwua4cGjBG0RbpvGa%2FOBk6CbZLpn7%2FLawxsVZ1U1ksGd8HGJ%2FGMYDOauM%2FhRGNWRFsXnn%2BsrPhaJ3SoirVeV3q9JVrjGT6%2FUT3W9qIDtdPP4MJae5mp6TG5fusJjkCLxLTbeXF0%2FhbwEnAA54uj3jpTsh7rXVDB%2B8skGSdMhIITz3%2ByS%2BdMqt7iEgFOWqYXGwgXLGbOqyGGz2ikth4cs1FMT4sYrA066%2BcMkE9q3l3bsFZHQMw13UPgJQp2f69JIzgHbZ%2FoCkdDYNxUutRhZ6cMitSLrIGtcAa7p%2Fevtejnw5eTz20kLNAxjB3CMUuS1H5qhxb6cSmxneilYH1WINNPjCrDPCJ3FxlKtCJo4QzIfIKogegd%2B44T78fQzt8RP7LfA%2FzjITD9bdiCYW0f81Q3O8zzL7l7RtfnLfYXAuTFh9GtAdE8D6b4F2pnXkMwrfCCwAY6sgG4%2BnyhdUNH%2FhdcK7GZ56erHPDOYF04vpG2hZy26v7cSnA3Xb7zrqVzkLxPdyAViJnMjzV1c8itVIHgnkLuA0C%2FPJrp3RPy0ivl9dofnd%2FLtoBkoBadnTgx2f7x4SZ62bdbWk5DJ%2FavMuOajJ%2F4tl9%2F7%2FLWoyi92xH2ZCvnT4wIIakx9ODzn2dRwSYwP20omrw5oAHK8KfXr39zDhQcs6FZMnWqYVxGlKHy0XIqJY8mTLeE&X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Date=20250417T092301Z&X-Amz-SignedHeaders=host&X-Amz-Expires=300&X-Amz-Credential=ASIAQ3PHCVTYXYYPN3VY%2F20250417%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Signature=bfaba1511f86b81dca784d137decbc826e87f24f30b4589d02ca67033b886446&hash=5e1a53dc45c46e59555d516468c6e00966f056dd50abcad3b239f603507d92a7&host=68042c943591013ac2b2430a89b270f6af2c76d8dfd086a07176afe7c76c2c61&pii=S0898122116300761&tid=spdf-f29c66e1-0a20-4e85-99d8-adbf3bfe5f8e&sid=68bdc92a37ea6249ef2b6425bac44510ce06gxrqb&type=client&tsoh=d3d3LnNjaWVuY2VkaXJlY3QuY29t&rh=d3d3LnNjaWVuY2VkaXJlY3QuY29t&ua=1d035b54560505005402&rr=931adc7ff9ea1b7e&cc=gb
 
 # Solving -u.laplace = pi^2 * sin(pix), 0<x<2
 # Constant zero Dirichlet BCs.
@@ -94,3 +90,9 @@ def analytical(x):
 
 assert slope > 1.9
 assert slope < 2.1
+
+
+# ref - An efficient implementation of fourth-order compact finite
+# difference scheme for Poisson equation with Dirichlet
+# boundary conditions
+# https://pdf.sciencedirectassets.com/271503/1-s2.0-S0898122116X00090/1-s2.0-S0898122116300761/main.pdf?X-Amz-Security-Token=IQoJb3JpZ2luX2VjENH%2F%2F%2F%2F%2F%2F%2F%2F%2F%2FwEaCXVzLWVhc3QtMSJGMEQCIH%2F3tmUOBoYeSj%2FmE47N2yF5b5IvdQohX7JWYGRZhNFTAiAX0P7PnDVBy2D3pTjsZRFgFEh7nVzV8oxoUFevx%2FucWSqzBQhZEAUaDDA1OTAwMzU0Njg2NSIMHBuo2q1ii2daxqQiKpAFdh12uILywsya5sCFqTsXwUdb%2F4lvp5uZIigTrb18tjM8cPru5xrZDgDVYVIlT4G6L1SE05FkjWKSQ7AO24wec3y2bNKGgpC%2FPFbEmTv6CgpR%2FjoSpToblGLiOqTkUgICSK3EoMdah%2BWl552nO0Ajdwuor0brGfDM7C2fgH1FqM%2BLyJ2do33rYFjGAswzZsQOGcf%2BChdkaKA9bxvfgE2%2Bukf7RcYDsGteSEX5Zb9XoyvMheiMUZoZk7KVPWjj3JORx9qetLs9LkpPO3IU%2BqPxtM7Vt3BnEnXR9gQ2bnL%2FtcT%2FcvsZ7a8AdiU1j8%2F%2Fxi9nBgPow0MQTmaoe9n67XRS0BVE7wAWldDb2qdZuOfwYl%2F2iG78mMTn%2FC4YcOCezc4nUT9fTcTcv3wKZzA%2Bkh8Z%2BXvdTcdADCKdVaIXLylqlhEmBlwua4cGjBG0RbpvGa%2FOBk6CbZLpn7%2FLawxsVZ1U1ksGd8HGJ%2FGMYDOauM%2FhRGNWRFsXnn%2BsrPhaJ3SoirVeV3q9JVrjGT6%2FUT3W9qIDtdPP4MJae5mp6TG5fusJjkCLxLTbeXF0%2FhbwEnAA54uj3jpTsh7rXVDB%2B8skGSdMhIITz3%2ByS%2BdMqt7iEgFOWqYXGwgXLGbOqyGGz2ikth4cs1FMT4sYrA066%2BcMkE9q3l3bsFZHQMw13UPgJQp2f69JIzgHbZ%2FoCkdDYNxUutRhZ6cMitSLrIGtcAa7p%2Fevtejnw5eTz20kLNAxjB3CMUuS1H5qhxb6cSmxneilYH1WINNPjCrDPCJ3FxlKtCJo4QzIfIKogegd%2B44T78fQzt8RP7LfA%2FzjITD9bdiCYW0f81Q3O8zzL7l7RtfnLfYXAuTFh9GtAdE8D6b4F2pnXkMwrfCCwAY6sgG4%2BnyhdUNH%2FhdcK7GZ56erHPDOYF04vpG2hZy26v7cSnA3Xb7zrqVzkLxPdyAViJnMjzV1c8itVIHgnkLuA0C%2FPJrp3RPy0ivl9dofnd%2FLtoBkoBadnTgx2f7x4SZ62bdbWk5DJ%2FavMuOajJ%2F4tl9%2F7%2FLWoyi92xH2ZCvnT4wIIakx9ODzn2dRwSYwP20omrw5oAHK8KfXr39zDhQcs6FZMnWqYVxGlKHy0XIqJY8mTLeE&X-Amz-Algorithm=AWS4-HMAC-SHA256&X-Amz-Date=20250417T092301Z&X-Amz-SignedHeaders=host&X-Amz-Expires=300&X-Amz-Credential=ASIAQ3PHCVTYXYYPN3VY%2F20250417%2Fus-east-1%2Fs3%2Faws4_request&X-Amz-Signature=bfaba1511f86b81dca784d137decbc826e87f24f30b4589d02ca67033b886446&hash=5e1a53dc45c46e59555d516468c6e00966f056dd50abcad3b239f603507d92a7&host=68042c943591013ac2b2430a89b270f6af2c76d8dfd086a07176afe7c76c2c61&pii=S0898122116300761&tid=spdf-f29c66e1-0a20-4e85-99d8-adbf3bfe5f8e&sid=68bdc92a37ea6249ef2b6425bac44510ce06gxrqb&type=client&tsoh=d3d3LnNjaWVuY2VkaXJlY3QuY29t&rh=d3d3LnNjaWVuY2VkaXJlY3QuY29t&ua=1d035b54560505005402&rr=931adc7ff9ea1b7e&cc=gb

examples/petsc/Poisson/04_poisson.py

@@ -0,0 +1,123 @@
+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'
+
+
+# 2D test
+# Solving u.laplace = 0
+# Dirichlet BCs.
+# ref - https://www.scirp.org/journal/paperinformation?paperid=113731#f2
+# example 2 -> note they wrote u(x,1) bc wrong, it should be u(x,y) = e^-pi*sin(pix)
+
+
+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(-y*np.pi)) * np.float64(np.sin(np.pi*x))
+
+
+Lx = np.float64(1.)
+Ly = np.float64(1.)
+
+n_values = [10, 30, 50, 70, 90, 110]
+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
+    )
+
+    u = Function(name='u', grid=grid, space_order=2)
+    rhs = Function(name='rhs', grid=grid, space_order=2)
+
+    eqn = Eq(rhs, u.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[:] = 0.
+
+    bcs = Function(name='bcs', grid=grid, space_order=2)
+
+    bcs.data[:, 0] = np.sin(np.pi*tmpx)
+    bcs.data[:, -1] = np.exp(-np.pi)*np.sin(np.pi*tmpx)
+    bcs.data[0, :] = 0.
+    bcs.data[-1, :] = 0.
+    
+    # # Create boundary condition expressions using subdomains
+    bc_eqns = [EssentialBC(u, bcs, subdomain=sub1)]
+    bc_eqns += [EssentialBC(u, bcs, subdomain=sub2)]
+    bc_eqns += [EssentialBC(u, bcs, subdomain=sub3)]
+    bc_eqns += [EssentialBC(u, bcs, subdomain=sub4)]
+
+    exprs = [eqn]+bc_eqns
+    petsc = PETScSolve(exprs, target=u, solver_parameters={'ksp_rtol': 1e-6})
+
+    with switchconfig(language='petsc'):
+        op = Operator(petsc)
+        op.apply()
+
+    u_exact = analytical(X, Y)
+
+    diff = u_exact[1:-1] - u.data[1:-1]
+    error = np.linalg.norm(diff) / np.linalg.norm(u_exact[1:-1])
+    errors.append(error)
+
+slope, _ = np.polyfit(np.log(dx), np.log(errors), 1)
+
+
+assert slope > 1.9
+assert slope < 2.1