How to use slepc4py - 10 common examples

# Now solve eigenvalue problem for $F_u(u, \lambda)\phi = r\phi$
    # Want eigenmode phi with minimal eigenvalue r
    B = derivative(residual(th, lm, TestFunction(V)), th, TrialFunction(V))

    petsc_M = assemble(inner(TestFunction(V), TrialFunction(V))*dx, bcs=bcs).petscmat
    petsc_B = assemble(B, bcs=bcs).petscmat

    num_eigenvalues = 1

    opts = PETSc.Options()
    opts.setValue("eps_target_magnitude", None)
    opts.setValue("eps_target", 0)
    opts.setValue("st_type", "sinvert")

    es = SLEPc.EPS().create(comm=COMM_WORLD)
    es.setDimensions(num_eigenvalues)
    es.setOperators(petsc_B, petsc_M)
    es.setProblemType(SLEPc.EPS.ProblemType.GHEP)
    es.setFromOptions()
    es.solve()

    ev_re, ev_im = petsc_B.getVecs()
    es.getEigenpair(0, ev_re, ev_im)
    eigenmode = Function(V)
    eigenmode.vector().set_local(ev_re)

    Z = MixedFunctionSpace([V, R, V])

    # Set initial guesses for state, parameter, null eigenmode
    z = Function(Z)
    z.split()[0].assign(th)

B = derivative(residual(th, lm, TestFunction(V)), th, TrialFunction(V))

    petsc_M = assemble(inner(TestFunction(V), TrialFunction(V))*dx, bcs=bcs).petscmat
    petsc_B = assemble(B, bcs=bcs).petscmat

    num_eigenvalues = 1

    opts = PETSc.Options()
    opts.setValue("eps_target_magnitude", None)
    opts.setValue("eps_target", 0)
    opts.setValue("st_type", "sinvert")

    es = SLEPc.EPS().create(comm=COMM_WORLD)
    es.setDimensions(num_eigenvalues)
    es.setOperators(petsc_B, petsc_M)
    es.setProblemType(SLEPc.EPS.ProblemType.GHEP)
    es.setFromOptions()
    es.solve()

    ev_re, ev_im = petsc_B.getVecs()
    es.getEigenpair(0, ev_re, ev_im)
    eigenmode = Function(V)
    eigenmode.vector().set_local(ev_re)

    Z = MixedFunctionSpace([V, R, V])

    # Set initial guesses for state, parameter, null eigenmode
    z = Function(Z)
    z.split()[0].assign(th)
    z.split()[1].assign(lm)
    z.split()[2].assign(eigenmode)

# We just need the Stiffness and Mass matrix
    a = inner(grad(u), grad(v))*dx
    m = inner(u, v)*dx

    A = topetsc(assemble(a, bcs=[bc]))
    M = topetsc(assemble(m, bcs=[bc]))

    # This shifts the "1.0" Eigenvalues out of the spectrum
    # of interest (which is around 0.0 in this case).
    vals = np.repeat(1E8, len(bc.nodes))
    A.setValuesLocalRCV(bc.nodes.reshape(-1, 1),
                        bc.nodes.reshape(-1, 1),
                        vals.reshape(-1, 1))
    A.assemble()

    E = SLEPc.EPS()
    E.create(comm=mesh.comm)
    E.setOperators(A, M)
    st = E.getST()
    st.setType('sinvert')
    kspE = st.getKSP()
    kspE.setType('fgmres')
    E.setDimensions(5, PETSc.DECIDE)
    E.solve()

    nconv = E.getConverged()
    assert(nconv > 0)

    # Create the results vectors
    vr, wr = A.getVecs()
    vi, wi = A.getVecs()
    ev = []

def eigenvalues(self):
        import slepc4py
        slepc4py.init()
        from slepc4py import SLEPc
        E = SLEPc.EPS(); E.create()

        E.setOperators(self.A)
        E.setProblemType(SLEPc.EPS.ProblemType.NHEP)
        E.setDimensions(1)
        E.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)
        E.setFromOptions()
        E.solve()
        n = E.getConverged()
        for i in range(0, n):
            print(E.getEigenvalue(0))
        return

def eigenvalues(self):
        import slepc4py
        slepc4py.init()
        from slepc4py import SLEPc
        E = SLEPc.EPS(); E.create()

        E.setOperators(self.A)
        E.setProblemType(SLEPc.EPS.ProblemType.NHEP)
        E.setDimensions(1)
        E.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)
        E.setFromOptions()
        E.solve()
        n = E.getConverged()
        for i in range(0, n):
            print(E.getEigenvalue(0))
        return

def init_petsc_slepc(self, comm):
        import os
        petsc_arch = os.environ.get("PETSC_ARCH", None)

        import petsc4py
        petsc4py.init(args=['-no_signal_handler'], arch=petsc_arch, comm=comm)
        import slepc4py
        slepc4py.init(args=['-no_signal_handler'], arch=petsc_arch)

        self._comm = comm
        self._PETSc = petsc4py.PETSc
        self._SLEPc = slepc4py.SLEPc

def init_petsc_slepc(self, comm):
        import os
        petsc_arch = os.environ.get("PETSC_ARCH", None)

        import petsc4py
        petsc4py.init(args=['-no_signal_handler'], arch=petsc_arch, comm=comm)
        import slepc4py
        slepc4py.init(args=['-no_signal_handler'], arch=petsc_arch)

        self._comm = comm
        self._PETSc = petsc4py.PETSc
        self._SLEPc = slepc4py.SLEPc

def init_slepc_args():
    try:
        import sys, slepc4py

    except ImportError:
        return

    argv = [arg for arg in sys.argv if arg not in ['-h', '--help']]
    slepc4py.init(argv)

conjunction with :class:`emopt.fdfd.FDFD` to simulated waveguide structures
which are particularly interesting for applications in silicon photonics, etc.

References
----------
[1] A. B. Fallahkhair, K. S. Li and T. E. Murphy, "Vector Finite Difference
Modesolver for Anisotropic Dielectric Waveguides", J. Lightwave Technol. 26(11),
1423-1431, (2008).
"""
from __future__ import absolute_import
# Initialize petsc first
from builtins import range
from builtins import object
import sys, slepc4py
from future.utils import with_metaclass
slepc4py.init(sys.argv)

from .defs import FieldComponent
from .misc import info_message, warning_message, error_message, \
NOT_PARALLEL, run_on_master, MathDummy

from . import grid

from abc import ABCMeta, abstractmethod
from petsc4py import PETSc
from slepc4py import SLEPc
from mpi4py import MPI
import numpy as np

__author__ = "Andrew Michaels"
__license__ = "GPL License, Version 3.0"
__version__ = "2019.5.6"

def eigenvalues(self):
        import slepc4py
        slepc4py.init()
        from slepc4py import SLEPc
        E = SLEPc.EPS(); E.create()

        E.setOperators(self.A)
        E.setProblemType(SLEPc.EPS.ProblemType.NHEP)
        E.setDimensions(1)
        E.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_REAL)
        E.setFromOptions()
        E.solve()
        n = E.getConverged()
        for i in range(0, n):
            print(E.getEigenvalue(0))
        return