first step with real-part operation correction with Larange interpolation: proof of concept

bin/args.py

@@ -10,7 +10,7 @@ See the LICENSE.txt and README.md for more license and copyright information.
 import argparse
 
 def parse():
-    """Parses the python script arguments
+    """Parses python script arguments
 
     Returns an argparse
     """
@@ -69,8 +69,7 @@ desc = ("a non-linear eigenvalue solver " +
 
 
 class MyFormatter(argparse.HelpFormatter):
-    """ Custom Formatter
-    """
+    """ Custom Formatter"""
 
     def __init__(self, prog):
         super(MyFormatter, self).__init__(prog,

bin/beyn.py

@@ -44,9 +44,6 @@ class Beyn:
         self.__rankTol  = rankTol
         self.__verbose  = verbose
 
-        # Display the parameters
-        if(verbose): self.__display()
-
     def simple(self):
         """Solves the eigenvalue problem using Beyn's algorithm (simple version)
 
@@ -131,6 +128,32 @@ class Beyn:
         if(verbose): print "Done!"
         return myLambda, Q, norm
 
+    def quite(self, beQuite):
+        """Sets this Beyn solver verbosity
+
+        If beQuite is True, this solver becomes quite; it is verbose otherwise
+        """
+        self.__verbose = not beQuite
+
+    def display(self):
+        """Prints this Beyn solver parameters"""
+        print "Beyn's contour integral method (simple)"
+        print "---------------------------------------"
+        print " # Nodes used for the trapezoidal rule:"+ " " +str(self.__nNode)
+        print " # Maximum number of iterations:       "+ " " +str(self.__maxIt)
+        print " # Initial size of col(A0):            "+ " " +str(self.__lStart)
+        print " # Step size for increasing col(A0):   "+ " " +str(self.__lStep)
+        print " # Rank test relative tolerance:       ",
+        print format(self.__rankTol, ".2e")
+        print "---------------------------------------"
+        print " # Cirular path origin:                ",
+        print "(" + format(np.real(self.__origin), "+.2e") + ")",
+        print "+",
+        print "(" + format(np.imag(self.__origin), "+.2e") + ")j"
+        print " # Cirular path radius:                ",
+        print format(self.__radius, "+.2e")
+        print "---------------------------------------"
+
     def __rank(self, A):
         """Returns the rank of a given matrix and its singular values"""
         S   = np.linalg.svd(A, full_matrices=False, compute_uv=0)
@@ -195,22 +218,3 @@ class Beyn:
     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(self):
-        """Prints this Beyn parameters"""
-        print "Beyn's contour integral method (simple)"
-        print "---------------------------------------"
-        print " # Nodes used for the trapezoidal rule:"+ " " +str(self.__nNode)
-        print " # Maximum number of iterations:       "+ " " +str(self.__maxIt)
-        print " # Initial size of col(A0):            "+ " " +str(self.__lStart)
-        print " # Step size for increasing col(A0):   "+ " " +str(self.__lStep)
-        print " # Rank test relative tolerance:       ",
-        print format(self.__rankTol, ".2e")
-        print "---------------------------------------"
-        print " # Cirular path origin:                ",
-        print "(" + format(np.real(self.__origin), "+.2e") + ")",
-        print "+",
-        print "(" + format(np.imag(self.__origin), "+.2e") + ")j"
-        print " # Cirular path radius:                ",
-        print format(self.__radius, "+.2e")
-        print "---------------------------------------"

bin/cim.py

@@ -8,10 +8,11 @@ Technische Universitaet Darmstadt.
 See the LICENSE.txt and README.md for more license and copyright information.
 """
 
-import solver as Solver
-import beyn   as Beyn
-import numpy  as np
+from solver       import GetDPWave
+from beyn         import Beyn
+from fixed_point  import FixedPoint
 import args       as args
+import numpy      as np
 import time       as timer
 
 # Start timer
@@ -28,18 +29,25 @@ for i in arg.setnumber:
         setnumber.append(j)
 
 # Initialise GetDP
-operator = Solver.GetDPWave(arg.pro, arg.mesh, arg.resolution, setnumber)
+operator = GetDPWave(arg.pro, arg.mesh, arg.resolution, setnumber)
 
 # Initialise Beyn's solver
-beyn = Beyn.Beyn(operator, arg.origin, arg.radius,
+beyn = Beyn(operator, arg.origin, arg.radius,
             arg.nodes, arg.maxIt, arg.lStart, arg.lStep, arg.rankTol,
             not arg.quiet)
 
+# Display parameters
+if(not arg.quiet):
+    beyn.display()
+
 # Compute
 if(not arg.non_holomorphic):
     l, V, r = beyn.simple()
 else:
-    raise ValueError('Non-holomorphic case not implemented... yet')
+    fixed = FixedPoint(beyn) # This variable ('fixed') is meant to be called
+    fixed.solve()            # elsewhere (i.e. Python[] call in GetDP) thanks to
+                             # the python interpreter.
+                             # Please don't change its name!
 
 # Display the computed eigenvalues
 if(not arg.quiet):

bin/fixed_point.py

+"""
+cim.py, a non-linear eigenvalue solver.
+Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
+Institut fuer Theorie Elektromagnetischer Felder (TEMF),
+Technische Universitaet Darmstadt.
+
+See the LICENSE.txt and README.md for more license and copyright information.
+"""
+
+import beyn as Beyn
+
+"""
+At some point I will need:
+*/ a cluster identification
+*/ an error with respect to the *NON*-holomorphic problem
+*/ a way to switch between holomorphic and non-holomorphic runs
+*/ a corrector?
+"""
+
+class FixedPoint:
+    """A fixed-point solver for handling non-holomorphic operator
+
+    The fixed-point approach follows the strategy bellow:
+    TODO
+
+    if cluster = [], use corrector initial guess
+    if needTrueResidual = True, corrector becomes non-holomorphic
+    """
+
+    def __init__(self, solver):
+        """Instantiates a new FixedPoint non-holomorphic solver
+
+        Keyword arguments:
+        solver -- the solver used to compute holomorphic solutions
+        """
+        # Save
+        self.__solver = solver
+
+        # State
+        self.__cluster     = []
+        self.__holomorphic = True
+
+    def solve(self):
+        """Solves the eigenvalue problem by a fixed-point iteration
+
+        It returns the computed eigenvalues, eigenvectors, TODO
+        """
+        # (Re)initialize
+        self.__cluster     = []
+        self.__holomorphic = True
+
+        # First run with corrector initial guess
+        l, V, r = self.__solver.simple()
+
+        print l
+
+    def cluster(self):
+        """Returns the eigenvalues (by cluster) found by the last iteration"""
+        return self.__cluster
+
+    def needNonHolomorphic(self):
+        """This method returns True if this FixedPoint needs the true,
+        non-holomorphic, operator and not its holomorphic approximation"""
+        return self.__trueRes

bin/interpolator.py

+"""
+cim.py, a non-linear eigenvalue solver.
+Copyright (C) 2017 N. Marsic, F. Wolf, S. Schoeps and H. De Gersem,
+Institut fuer Theorie Elektromagnetischer Felder (TEMF),
+Technische Universitaet Darmstadt.
+
+See the LICENSE.txt and README.md for more license and copyright information.
+"""
+
+from numpy import poly1d
+from scipy import interpolate
+from abc   import ABCMeta, abstractmethod
+
+class Interpolator:
+    """A generic class for interpolation"""
+    @abstractmethod
+    def __init__(self):
+        pass
+
+    def at(self, x):
+        """Returns the value of this Interpolation at a given point"""
+        return self._p(x) # We use poly1d from scipy
+
+    __metaclass__ = ABCMeta
+
+
+class Lagrange(Interpolator):
+    """A class for Lagrange interpolation"""
+    def __init__(self, xdata, ydata):
+        """Instantiates a new Interpolator, such that:
+        p(xdata) = ydata for all xdata, where p is build with a Lagrange basis
+        """
+        self._p = interpolate.lagrange(xdata, ydata)

bin/solver.py

@@ -62,7 +62,7 @@ class GetDPWave:
         return self.solution(0).shape[0]
 
     def apply(self, x, w):
-        """Applies x to the operator with a pulsation of w
+        """Applies x to the operator with an angular frequency of w
 
         This method updates self.solution()
         """

example/maxwell_sibc_lagrange/Makefile

@@ -6,12 +6,12 @@ init:
 	gmsh sphere.geo -3
 
 cim:
-	python ${BIN}/cim.py sphere.pro sphere.msh Solve 1e10+1e9j 4e9 -nodes 20 -lStart 20 -setnumber isSC 0
+#	python ${BIN}/cim.py sphere.pro sphere.msh Solve 1e10+1e9j 4e9 -nodes 30 -lStart 20 -lStep 30 -setnumber isSC 0 -non_holomorphic
 
-#	python ${BIN}/cim.py sphere.pro sphere.msh Solve 7.4e9+7.3e8j 1e9 -nodes 10 -lStart 4 -setnumber isSC 0
+	python ${BIN}/cim.py sphere.pro sphere.msh Solve 7.4e9+7.3e8j 1e9 -nodes 10 -lStart 4 -setnumber isSC 0 -non_holomorphic
 #7.474355560164e+09 + 7.742536715699e+08j
 
-#	python ${BIN}/cim.py sphere.pro sphere.msh Solve 1.0e10+1.1e9j 1e9 -nodes 10 -lStart 6 -setnumber isSC 0
+#	python ${BIN}/cim.py sphere.pro sphere.msh Solve 1.0e10+1.1e9j 1e9 -nodes 10 -lStart 6 -setnumber isSC 0 -non_holomorphic
 #1.050336615299e+10 + 1.186115651023e+09j
 
 clean:

example/maxwell_sibc_lagrange/cimResolution.pro

@@ -7,6 +7,10 @@ Resolution{
   { Name Solve;
     System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
     Operation{
+      If(!doInit)
+        Evaluate[Python[]{"real_corrector_init.py"}];
+      EndIf
+
       Generate[A];
 
       If(doInit)

example/maxwell_sibc_lagrange/real_corrector_handl.py

+"""Real-part operator corrector
+
+Approximates the real-part operator by an holomorphic function, constructed with
+a Lagrange polynomial defined at some correction points. Please make sure that
+GetDP calls real_corrector_init.py before calling this file.
+
+The naming convention below must be followed:
+-- a FixedPoint object named 'fixed' is available and prepared by cim.py
+-- a Lagrange Interpolator object named 'lagrange' coming from initialization
+-- the variable 'output' is read back by GetDP as the Pyhton[] call output
+"""
+
+# Input angular frequency
+omega = input[0]
+
+# Compute holomorphic approximation of the real part operator
+if(not fixed.cluster()):
+    output = omega                            # If no cluster: Re(omega) ~ omega
+else:
+    output = (omega + lagrange.at(omega)) / 2 # Else: correct

example/maxwell_sibc_lagrange/real_corrector_init.py

+"""Real-part operator corrector: initialization
+
+See: real_corrector_handl.py for more information
+"""
+
+from interpolator import Lagrange
+
+# Previous eigenvalue clusters
+cluster = fixed.cluster()
+
+# Initialize interpolator
+if(cluster):
+    lagrange = Lagrange(cluster, np.conj(cluster).tolist())

example/maxwell_sibc_lagrange/sphere.pro

@@ -16,27 +16,24 @@ Function{
   Mu0  = 4*Pi*1e-7;      // [H/m]
   Eps0 = 1/(Mu0 * C0^2); // [F/m]
 
+  // Correction for real-part operator //
+  If(doInit)
+    PyPulsRe[] = 1;
+  EndIf
+  If(!doInit)
+    PyPulsRe[] = Python[Puls[]]{"real_corrector_handl.py"};
+  EndIf
+
   // Conductivity (coming from sphere.dat) //
   If(isSC == 1)                                             // Superconductor
-    Sigma[] = Sigma - J[] / (Mu0 * Lambda^2 * Re[Puls[]]); // [S/m]
+    Sigma[] = Sigma - J[] / (Mu0 * Lambda^2 * PyPulsRe[]]); // [S/m]
   EndIf
   If(isSC == 0)                                             // Normal conductor
     Sigma[] = Sigma;                                        // [S/m]
   EndIf
 
   // Leontovich SIBC (H-Formulation) //
-  PulsG1[]   = 7.474355560164e+09 + 7.742536715699e+08*J[];
-  PulsG2[]   = 1.050336615299e+10 + 1.186115651023e+09*J[];
-
-  PulsBar1[] = Re[PulsG1[]] - J[]*Im[PulsG1[]];
-  PulsBar2[] = Re[PulsG2[]] - J[]*Im[PulsG2[]];
-
-  L1[] = (Puls[] - PulsG2[]) / (PulsG1[] - PulsG2[]);
-  L2[] = (Puls[] - PulsG1[]) / (PulsG2[] - PulsG1[]);
-
-  //PulsCor[] = PulsBar1[]*Puls[]/PulsG1[]; //Re[Puls[]] + J[]*Im[Puls[]];
-  PulsCor[] = PulsBar1[]*L1[] + PulsBar2[]*L2[];
-  SL[] = Sqrt[(J[] * (Puls[]+PulsCor[])/2 * Mu0) / Sigma[]]; // H-Formulation
+  SL[] = Sqrt[(J[] * PyPulsRe[] * Mu0) / Sigma[]]; // H-Formulation
   SIbc[] = -1.0 / SL[];                            // E-Formulation
 }