An implemetation of the contour integral method using GetDP python bindings.

The python script uses GetDP to represent the operator considered in the CIM.
It seems OK.

beyn.py

+import solver as Solver
+import numpy  as np
+import numpy.matlib
+
+def simple(operator, radius, origin,
+           nodes=100, maxIt=10, lStart=1, lStep=1, rankTol=1e-4):
+    """Solve an eigenvalue problem using Beyn's algorithm (simple version)
+
+    Keyword arguments:
+    operator -- the solver defining the operator to use
+    radius -- the radius of the circular contour used to search the eigenvalues
+    origin -- the origin (in the complex plane) of the above circular contour
+    nodes -- the number of nodes for the trapezoidal integration rule (optional)
+    lStart -- the number of columns used for A0 when algorithm starts (optional)
+    lStep -- the step used for increasing the number of columns of A0 (optional)
+    rankTol -- the tolerance for the rank test
+    """
+
+    # Initialise A0 search
+    myPath = path(nodes, origin, radius)
+    hasK   = False
+    it     = 0
+    m      = operator.size()
+    l      = lStart
+    k      = -1
+
+    # Search A0
+    while(not hasK and it != maxIt):
+        print "Iteration: " + str(it+1)
+
+        vHat = randomMatrix(m ,l)                       # Take a random VHat
+        A0   = integrate(operator, myPath, 0, vHat)     # Compute A0
+        k    = np.linalg.matrix_rank(A0, tol=rankTol)   # Rank test
+
+        if(k == 0):                                     # Rank is zero?
+            raise RuntimeError(zeroRankErrorStr())      #  -> Stop
+        elif(k == m):                                   # Maximum rank reached?
+            raise RuntimeError(maxRankErrorStr())       #  -> Stop
+        elif(k < l):                                    # Found a null SV?
+            hasK = True                                 #  -> We have A0
+        else:                                           # Matrix is full rank?
+            l = l + lStep                               #  -> Increase A0 size
+
+        it += 1                                         # Keep on searching A0
+
+    # Check if maxIt was reached
+    if(it == maxIt):
+        raise RuntimeError(maxItErrorStr())
+    else:
+        print "Constructing linear EVP..."
+
+    # Compute V, S and Wh
+    #  NB: For SVD(A) = V*S*Wh, numpy computes {v, s, w}, such that:
+    #      v = V; diag(s) = S and w = Wh
+    V, S, Wh = np.linalg.svd(A0, full_matrices=False, compute_uv=1)
+
+    # Extract V0, W0 and S0Inv
+    V0    = np.delete(V,  l-1, 1)
+    W0    = np.delete(Wh, l-1, 0).H
+    S0Inv = np.matrix(np.diag(1/np.delete(S, l-1, 0)))
+
+    # Compute A1 and B
+    A1 = integrate(operator, myPath, 1, vHat)
+    B  = V0.H * A1 * W0 * S0Inv
+
+    # Eigenvalues of B
+    print "Solving linear EVP..."
+    myLambda, QHat = numpy.linalg.eig(B)
+
+    # Display
+    display(myLambda)
+
+    # Done
+    return
+
+
+## Import only simple (other functions are just helpers)
+__all__ = ['simple']
+
+
+## Helper functions
+def path(nodes, origin, radius):
+    """Returns a list with the coordinates of a circular contour
+
+    Keyword arguments:
+    nodes -- the number of nodes used to discretise the contour
+    radius -- the radius of the circular contour
+    origin -- the origin (in the complex plane) of the circular contour
+
+    The returned list contains nodes+1 elements,
+    such that the first and the last are identical (up the machine precision)
+    """
+
+    step         = 1.0 / nodes
+    nodesPlusOne = nodes + 1
+
+    path = list()
+    for i in range(nodesPlusOne):
+        path.append(origin + radius * np.exp(1j * 2 * np.pi * i * step))
+
+    return path
+
+
+def integrate(operator, path, order, vHat):
+    """Computes the countour integral of Beyn's method, that is matrix A_p
+
+    Keyword arguments:
+    operator -- the solver defining the operator to use
+    path -- the path to integrate on
+    order -- the order of Beyn's integral (that is the 'p' defining matrix A_p)
+    vHat -- the RHS matrix defing Beyn's integral
+    """
+
+    # Initialise I
+    I = np.matlib.zeros(vHat.shape, dtype=complex)
+
+    # Initialise integration loop
+    F1  = multiSolve(operator, vHat, path[0])
+    F1 *= np.power(path[0], order)
+
+    # Integration loop
+    pathSizeMinus = len(path) - 1;
+    for i in range(pathSizeMinus):
+        F2  = multiSolve(operator, vHat, path[i + 1])
+        F2 *= np.power(path[i + 1], order)
+
+        I += (F1 + F2) * (path[i + 1] - path[i])
+        F1 = F2
+
+    # Done
+    return I
+
+
+def multiSolve(solver, B, w):
+    """Solves for multiple RHS
+
+    Keyword arguments:
+    solver -- the solver to use
+    B -- the matrix of RHS
+    w -- the complex frequency to use
+    """
+
+    # Number of solve
+    nSolve = B.shape[1]
+
+    # Initialise X
+    size = solver.size()
+    X    = np.matlib.empty((size, nSolve), dtype=complex)
+
+    # Loop and solve
+    for i in range(nSolve):
+        b = B[:, i]
+        solver.solve(b, w)
+        X[:, i] = solver.solution()
+
+    # Done
+    return X
+
+
+def randomMatrix(n, m):
+    """Returns a random complex matrix of the given size (n x m)"""
+    return np.matlib.rand(n, m) + np.matlib.rand(n, m) * 1j
+
+
+def display(v):
+    """Displays the eigenvalues given in the array v"""
+    size = v.shape[0]
+    print "Eigenvalues:"
+    print "----------- "
+    for i in range(size):
+        print "(" + format(np.real(v[i]).tolist(), '+e') + ")",
+        print "+",
+        print "(" + format(np.imag(v[i]).tolist(), '+e') + ")*j"
+    print "----------- "
+    print "Found: " + str(size)
+
+
+def zeroRankErrorStr():
+    """Returns a string explaining the probable reason of a zero rank"""
+    return ("Found a rank of zero: " +
+            "the contour is probably enclosing no eigenvalues")
+
+
+def maxRankErrorStr():
+    """Returns a string explaining the probable reason of a maximal rank"""
+    return ("Maximal rank found: " +
+            "the contour is probably enclosing to many eigenvalues")
+
+
+def maxItErrorStr():
+    """Returns a string claiming: the maximum number of iterations is reached"""
+    return "Maximum number iterations is reached!"

main.py

+import solver as Solver
+import beyn   as Beyn
+import numpy  as np
+
+# Data
+pro        = "maxwell.pro"
+mesh       = "square.msh"
+resolution = "Maxwell"
+
+# Initialise GetDP
+operator = Solver.GetDPWave(pro, mesh, resolution)
+
+# Compute
+Beyn.simple(operator, 1e8, complex(9e8, 0))

maxwell.pro

+Group{
+  Omega = Region[1]; // Domain
+  Bndry = Region[2]; // Boundary
+}
+
+Function{
+  // Physical data //
+  DefineConstant[angularFreqRe = 5,  // [rad/s]
+                 angularFreqIm = 0]; // [rad/s]
+  C0    = 299792458;                 // [m/s]
+  aFC[] = Complex[angularFreqRe,
+                  angularFreqIm];    // [rad/m]
+  k[]   = aFC[]/C0;                  // [rad/m]
+
+  Printf["Angular frequency set to: %g + i*%g [rad/s]",
+         angularFreqRe, angularFreqIm];
+
+  // Algebraic data //
+  DefineConstant[x() = {},  // Solution
+                 b() = {},  // Right hand side
+                 d() = {}]; // DoFs
+
+  // Control data //
+  DefineConstant[imposeRHS = 0]; // Should I use an imposed RHS?
+}
+
+Jacobian{
+  { Name JVol;
+    Case{
+      { Region All; Jacobian Vol; }
+    }
+  }
+}
+
+Integration{
+  { Name IP2;
+    Case{
+      { Type Gauss;
+        Case{
+          { GeoElement Line;        NumberOfPoints 2; }
+          { GeoElement Triangle;    NumberOfPoints 3; }
+          { GeoElement Tetrahedron; NumberOfPoints 4; }
+        }
+      }
+    }
+  }
+}
+
+Constraint{
+  { Name Dirichlet; Type Assign;
+    Case{
+      { Region Bndry; Value 0; }
+    }
+  }
+}
+
+FunctionSpace{
+  { Name HCurl; Type Form1;
+    BasisFunction{ // Nedelec
+      { Name se; NameOfCoef ee; Function BF_Edge;
+        Support Omega; Entity EdgesOf[All]; }
+    }
+
+    Constraint{
+      { NameOfCoef ee; EntityType EdgesOf; NameOfConstraint Dirichlet; }
+    }
+  }
+}
+
+Formulation{
+  { Name Maxwell; Type FemEquation;
+    Quantity{
+      { Name e; Type Local; NameOfSpace HCurl; }
+    }
+    Equation{
+      Galerkin{ [         Dof{d e} , {d e}];
+        In Omega; Integration IP2; Jacobian JVol; }
+      Galerkin{ [-k[]^2 * Dof{  e} , {  e}];
+        In Omega; Integration IP2; Jacobian JVol; }
+    }
+  }
+}
+
+Resolution{
+  { Name Maxwell;
+    System{
+      { Name A; NameOfFormulation Maxwell; Type Complex; }
+    }
+    Operation{
+      Generate[A];
+
+      If(imposeRHS)
+        CopyRightHandSide[b(), A];
+      EndIf
+
+      Solve[A];
+      CopySolution[A, x()];
+    }
+  }
+}
+
+PostProcessing{
+  { Name Maxwell; NameOfFormulation Maxwell;
+    Quantity{
+      { Name e; Value{ Local{ [{e}]; In Omega; Jacobian JVol; } } }
+    }
+  }
+}
+
+PostOperation{
+  { Name Maxwell; NameOfPostProcessing Maxwell;
+    Operation{
+      Print[e, OnElementsOf Omega, File "maxwell.pos"];
+    }
+  }
+}

ref.pro

+Group{
+  Omega = Region[1]; // Domain
+  Bndry = Region[2]; // Boundary
+}
+
+Function{
+  // Physical data //
+  C0 = 299792458; // [m/s]
+
+  // Eigenvalue solver data //
+  DefineConstant[nEig   = 10,    // [-]
+                 target = 8e17]; // [(rad/s)^2]
+}
+
+Jacobian{
+  { Name JVol;
+    Case{
+      { Region All; Jacobian Vol; }
+    }
+  }
+}
+
+Integration{
+  { Name IP2;
+    Case{
+      { Type Gauss;
+        Case{
+          { GeoElement Line;        NumberOfPoints 2; }
+          { GeoElement Triangle;    NumberOfPoints 3; }
+          { GeoElement Tetrahedron; NumberOfPoints 4; }
+        }
+      }
+    }
+  }
+}
+
+Constraint{
+  { Name Dirichlet; Type Assign;
+    Case{
+      { Region Bndry; Value 0; }
+    }
+  }
+}
+
+FunctionSpace{
+  { Name HCurl; Type Form1;
+    BasisFunction{ // Nedelec
+      { Name se; NameOfCoef ee; Function BF_Edge;
+        Support Omega; Entity EdgesOf[All]; }
+    }
+
+    Constraint{
+      { NameOfCoef ee; EntityType EdgesOf; NameOfConstraint Dirichlet; }
+    }
+  }
+}
+
+Formulation{
+  { Name Eig; Type FemEquation;
+    Quantity{
+      { Name e; Type Local; NameOfSpace HCurl; }
+    }
+    Equation{
+      Galerkin{        [         Dof{d e} , {d e}];
+        In Omega; Integration IP2; Jacobian JVol; }
+      Galerkin{ DtDtDof[1/C0^2 * Dof{  e} , {  e}];
+        In Omega; Integration IP2; Jacobian JVol; }
+    }
+  }
+}
+
+Resolution{
+  { Name Eig;
+    System{
+      { Name A; NameOfFormulation Eig; Type Complex; }
+    }
+    Operation{
+      GenerateSeparate[A];
+      EigenSolve[A, nEig, target, 0];
+    }
+  }
+}
+
+PostProcessing{
+  { Name Eig; NameOfFormulation Eig;
+    Quantity{
+      { Name e; Value{ Local{ [{e}]; In Omega; Jacobian JVol; } } }
+    }
+  }
+}
+
+PostOperation{
+  { Name Eig; NameOfPostProcessing Eig;
+    Operation{
+      Print[e, OnElementsOf Omega, File "eig.pos", EigenvalueLegend];
+    }
+  }
+}

solver.py

+import getdp as GetDP
+import numpy as np
+
+class GetDPWave:
+    """A GetDP solver using complex arithmetic for time-harmonic wave problems
+
+    Only one instance of this class makes sens
+
+    The .pro file should use the following variables:
+    angularFreqRe -- the real part of the angular frequency to use
+    angularFreqIm -- the imaginary part of the angular frequency to use
+    b -- the right hand side (RHS) to use
+    x -- the solution vector computed by GetDP
+    imposeRHS -- should the imposed RHS be taken into account?
+    """
+
+    def __init__(self, pro, mesh, resolution):
+        """Instanciates a new SolverGetDP with a full '-solve'
+
+        Keyword arguments:
+        pro -- the .pro file to use
+        mesh -- the .msh file to use
+        resolution -- the resolution (from the .pro file) to use
+        """
+        # Save
+        self.__pro        = pro
+        self.__mesh       = mesh
+        self.__resolution = resolution
+
+        # Generate DoFs and initialise RHS and solution vectors in GetDP
+        GetDP.GetDP(["getdp",  self.__pro,
+                     "-msh",   self.__mesh,
+                     "-solve", self.__resolution,
+                     "-v", "2"])
+
+    def solution(self):
+        """Returns the solution"""
+        return self.__toNumpy(GetDP.GetDPGetNumber("x"))
+
+    def size(self):
+        """Returns the number of degrees of freedom"""
+        return self.solution().shape[0]
+
+    def solve(self, b, w):
+        """Solves with b as RHS and w as complex angular frequency"""
+        GetDP.GetDPSetNumber("angularFreqRe", np.real(w).tolist())
+        GetDP.GetDPSetNumber("angularFreqIm", np.imag(w).tolist())
+        GetDP.GetDPSetNumber("b", self.__toGetDP(b))
+        GetDP.GetDPSetNumber("imposeRHS", 1)
+        GetDP.GetDP(["getdp", self.__pro,
+                     "-msh",  self.__mesh,
+                     "-cal",  self.__resolution,
+                     "-v", "2"])
+
+    @staticmethod
+    def __toNumpy(vGetDP):
+        """Takes a GetDP list and returns a numpy array"""
+        size   = vGetDP.size() / 2
+        vNumpy = np.empty((size, 1), dtype=complex)
+
+        for i in range(size):
+            vNumpy[i] = complex(vGetDP[i*2], vGetDP[i*2 + 1])
+
+        return vNumpy
+
+    @staticmethod
+    def __toGetDP(vNumpy):
+        """Takes a numpy array and returns a GetDP list"""
+        size   = vNumpy.shape[0]
+        vGetDP = list()
+
+        for i in range(size):
+            vGetDP.append(float(np.real(vNumpy[i])))
+            vGetDP.append(float(np.imag(vNumpy[i])))
+
+        return vGetDP

square.geo

+DefineConstant[l   = {01, Name "Input/00Geometry/00Length [m]"},
+               msh = {10, Name "Input/01Mesh/00Density [point per length]"}];
+
+Point(1) = {+l/2, -l/2, 0, 1};
+Point(2) = {+l/2, +l/2, 0, 1};
+Point(3) = {-l/2, +l/2, 0, 1};
+Point(4) = {-l/2, -l/2, 0, 1};
+
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+
+Line     Loop(1) = {1, 2, 3, 4};
+Plane Surface(1) = {1};
+
+Transfinite    Line {1, 2, 3, 4} = (msh * l) Using Progression 1;
+Transfinite Surface {1};
+
+Physical Surface(1) = {1};
+Physical    Line(2) = {1, 2, 3, 4};