# ex13.py: Nonlinear eigenproblem with contour integral # ===================================================== # # This example solves a nonlinear eigenvalue problem with the `NEP # ` module configured to apply a contour integral # method. # # The problem arises from the PDE # # .. math:: # # u_{xx}(x) + n_c^2 \lambda^2 u(x) + g(\lambda) D_0 \lambda^2 u(x) = 0 # # where # # .. math:: # # g(\lambda) &= g_t/(\lambda-k_a + i g_t), \\ # k_a &= 8.0, \\ # g_t &= 0.5, \\ # D_0 &= 0.5, \\ # n_c &= 1.2, # # and the boundary conditions are # # .. math:: # # u(0) &= 0 \\ # u_x(1) &= i \lambda u(1). # # For the discretization, :math:`n` grid points are used, :math:`x_1=0,\dots,x_n=1`, # with step size :math:`h = 1/(n-1)`. Hence, # # .. math:: # # u_{xx}(x_i) &= \frac{1}{h^2} (u_{i-1} - 2 u_i + u_{i+1}) \\ # &= \frac{1}{h^2} [1, -2, 1] [u_{i-1}, u_i, u_{i+1}]^T. # # The boundary condition at :math:`x=0` is :math:`u_1=0`, and at :math:`x=1`: # # .. math:: # # u'(1) &= 1/2 ( u'(1+h/2) + u'(1-h/2) ) \\ # &= 1/2 ( (u_{n+1}-u_n)/h + (u_n-u_{n-1})/h ) \\ # &= 1/(2h) (u_{n+1} - u_{n-1}) = i\lambda u_n, # # and therefore # # .. math:: # # u_{n+1} = 2i h \lambda u_n + u_{n-1}. # # The Laplace term for :math:`u_n` is # # .. math:: # # &= \frac{1}{h^2} (u_{n-1} - 2u_n + u_{n+1}) \\ # &= \frac{1}{h^2} (u_{n-1} - 2u_n + 2ih\lambda u_n + u_{n-1}) \\ # &= \frac{1}{h^2} (2 u_{n-1} + (2ih\lambda - 2) u_n) \\ # &= \frac{1}{h^2} [2, 2ih\lambda -2] [u_{n-1}, u_n]^T. # # The above discretization allows us to write the nonlinear PDE in the # following split-operator form # # .. math:: # # \{A + \lambda^2 n_c^2 I_d + g(\lambda) \lambda^2 D_0 I_d + 2i\lambda/h D\} u = 0 # # so :math:`f_1 = 1`, :math:`f_2 = n_c^2 \lambda^2`, :math:`f_3 = g(\lambda) # \lambda^2D_0`, :math:`f_4 = 2i\lambda/h`, with coefficient matrices # # .. math:: # # A = \begin{bmatrix} # 1 & 0 & 0 & \dots & 0 \\ # 0 & * & * & \dots & 0 \\ # 0 & * & * & \dots & 0 \\ # \vdots & \vdots & \vdots & \ddots & \vdots \\ # 0 & 0 & 0 & \dots & * # \end{bmatrix}, # I_d = \begin{bmatrix} # 0 & 0 & 0 & \dots & 0 \\ # 0 & 1 & 0 & \dots & 0 \\ # 0 & 0 & 1 & \dots & 0 \\ # \vdots & \vdots & \vdots & \ddots & \vdots \\ # 0 & 0 & 0 & \dots & 0 # \end{bmatrix}, # D = \begin{bmatrix} # 0 & 0 & 0 & \dots & 0 \\ # 0 & 0 & 0 & \dots & 0 \\ # 0 & 0 & 0 & \dots & 0 \\ # \vdots & \vdots & \vdots & \ddots & \vdots \\ # 0 & 0 & 0 & \dots & 1 # \end{bmatrix}. # # Contributed by: Thomas Hisch. # # The full source code for this demo can be `downloaded here # <../_static/ex13.py>`__. # Initialization by importing slepc4py, petsc4py, numpy and scipy. import sys import slepc4py slepc4py.init(sys.argv) # isort:skip import numpy as np try: import scipy import scipy.optimize except ImportError: scipy = None from petsc4py import PETSc from slepc4py import SLEPc Print = PETSc.Sys.Print # Check that the selected PETSc/SLEPc are built with complex scalars. if not np.issubdtype(PETSc.ScalarType, np.complexfloating): Print("Demo should only be executed with complex PETSc scalars") exit(0) # Long function that defines the nonlinear eigenproblem and solves it. def solve(n): L = 1.0 h = L / (n - 1) nc = 1.2 ka = 10.0 gt = 4.0 D0 = 0.5 A = PETSc.Mat().create() A.setSizes([n, n]) A.setFromOptions() A.setOption(PETSc.Mat.Option.HERMITIAN, False) rstart, rend = A.getOwnershipRange() d0, d1, d2 = ( 1 / h**2, -2 / h**2, 1 / h**2, ) Print(f"dterms={(d0, d1, d2)}") if rstart == 0: # dirichlet boundary condition at the left lead A[0, 0] = 1.0 A[0, 1] = 0.0 A[1, 0] = 0.0 A[1, 1] = d1 A[1, 2] = d2 rstart += 2 if rend == n: # at x=1.0 neumann boundary condition (not handled here but in a # different matrix (D)) A[n - 1, n - 2] = 2.0 / h**2 A[n - 1, n - 1] = (-2) / h**2 # + 2j*k*h / h**2 (neumann) rend -= 1 for i in range(rstart, rend): A[i, i - 1 : i + 2] = [d0, d1, d2] A.assemble() Id = PETSc.Mat().create() Id.setSizes([n, n]) Id.setFromOptions() Id.setOption(PETSc.Mat.Option.HERMITIAN, True) rstart, rend = Id.getOwnershipRange() if rstart == 0: # due to dirichlet BC rstart += 1 for i in range(rstart, rend): Id[i, i] = 1.0 Id.assemble() D = PETSc.Mat().create() D.setSizes([n, n]) D.setFromOptions() D.setOption(PETSc.Mat.Option.HERMITIAN, True) _, rend = D.getOwnershipRange() if rend == n: D[n - 1, n - 1] = 1 D.assemble() Print(f"DOF: {A.getInfo()['nz_used']}, MEM: {A.getInfo()['memory']}") f1 = SLEPc.FN().create() f1.setType(SLEPc.FN.Type.RATIONAL) f1.setRationalNumerator([1.0]) f2 = SLEPc.FN().create() f2.setType(SLEPc.FN.Type.RATIONAL) f2.setRationalNumerator([nc**2, 0.0, 0.0]) f3 = SLEPc.FN().create() f3.setType(SLEPc.FN.Type.RATIONAL) f3.setRationalNumerator([D0 * gt, 0.0, 0.0]) f3.setRationalDenominator([1.0, -ka + 1j * gt]) f4 = SLEPc.FN().create() f4.setType(SLEPc.FN.Type.RATIONAL) f4.setRationalNumerator([2j / h, 0]) # Setup the solver nep = SLEPc.NEP().create() nep.setSplitOperator( [A, Id, Id, D], [f1, f2, f3, f4], PETSc.Mat.Structure.SUBSET, ) # Customize options nep.setTolerances(tol=1e-7) nep.setDimensions(nev=24) nep.setType(SLEPc.NEP.Type.CISS) # the rg params are chosen s.t. the singularity at k = ka - 1j*gt is # outside of the contour. radius = 3 * gt vscale = 0.5 * gt / radius rg_params = (ka, 3 * gt, vscale) R = nep.getRG() R.setType(SLEPc.RG.Type.ELLIPSE) Print(f"RG params: {rg_params}") R.setEllipseParameters(*rg_params) nep.setFromOptions() # Solve the problem nep.solve() its = nep.getIterationNumber() Print("Number of iterations of the method: %i" % its) sol_type = nep.getType() Print("Solution method: %s" % sol_type) nev, ncv, mpd = nep.getDimensions() Print("") Print("Subspace dimension: %i" % ncv) tol, maxit = nep.getTolerances() Print("Stopping condition: tol=%.4g" % tol) Print("") nconv = nep.getConverged() Print("Number of converged eigenpairs %d" % nconv) x = A.createVecs("right") evals = [] modes = [] if nconv > 0: Print() Print(" lam ||T(lam)x|| |lam-lam_exact|/|lam_exact| ") Print("--------------------- ------------- -----------------------------") for i in range(nconv): lam = nep.getEigenpair(i, x) error = nep.computeError(i) def eigenvalue_error_term(k): gkmu = gt / (k - ka + 1j * gt) nceff = np.sqrt(nc**2 + gkmu * D0) return -1j / np.tan(nceff * k * L) - 1 / nceff # compute the expected_eigenvalue # we assume that the numerically calculated eigenvalue is close to # the exact one, which we can determine using a Newton-Raphson # method. if scipy: expected_lam = scipy.optimize.newton( eigenvalue_error_term, np.complex128(lam), rtol=1e-11 ) rel_err = abs(lam - expected_lam) / abs(expected_lam) rel_err = "%6g" % rel_err else: rel_err = "scipy not installed" Print(" %9f%+9f j %12g %s" % (lam.real, lam.imag, error, rel_err)) evals.append(lam) modes.append(x.getArray().copy()) Print() return np.asarray(evals), rg_params, ka, gt # The main function reads the problem size ``n`` from the command line, # solves the problem with the above function, and then plots the computed # eigenvalues. def main(): opts = PETSc.Options() n = opts.getInt("n", 256) Print(f"n={n}") evals, rg_params, ka, gt = solve(n) if not opts.getBool("ploteigs", True) or PETSc.COMM_WORLD.getRank(): return try: import matplotlib.pyplot as plt from matplotlib.patches import Ellipse except ImportError: print("plot is not shown, because matplotlib is not installed") else: fig, ax = plt.subplots() ax.plot(evals.real, evals.imag, "x") height = 2 * rg_params[1] * rg_params[2] ellipse = Ellipse( xy=(rg_params[0], 0.0), width=rg_params[1] * 2, height=height, edgecolor="r", fc="None", lw=2, ) ax.add_patch(ellipse) ax.grid() ax.legend() plt.show() if __name__ == "__main__": main()