Reuse the LU decomposition the GetDP way: default behaviour

Merge commit '9e5bef6d'

README.md

@@ -28,11 +28,13 @@ cim.py and [GetDP](http://getdp.info) are exchanging the following parameters:
 
     angularFreqRe  the real part of the frequency
     angularFreqIm  the imaginary part of the frequency
-    x()            the solution vector of the linear system
-    b()            the right hand side (RHS) of the linear system
-    imposeRHS      flag signalling if an other RHS should be imposed
-    doPostpro      flag signalling if only a post-processing should be done
-    doApply        flag signalling if x() should be applied to the linear system
+    nRHS           the number of right hand side that should be considered (default: 1)
+    b_i()/b~{i}()  the ith right hand side (RHS) of the linear system (first index is 0)
+    x_i()/x~{i}()  the ith solution vector of the linear system (first index is 0)
+    doInit         flag signalling if an initialisation should be done (default: 0)
+    doSolve        flag signalling if a solve should be done for all RHS (default: 0)
+    doPostpro      flag signalling if a post-processing should be done (default: 0)
+    doApply        flag signalling if x() should be applied to the linear system (default: 0)
     fileName       name of the post-processing file
 
 More information can be found in [bin/solver.py](bin/solver.py). Furthermore,

bin/beyn.py

@@ -101,7 +101,7 @@ def simple(operator, origin, radius,
     norm   = np.empty((nFound))
     for i in range(nFound):
         operator.apply(Q[:, i], myLambda[i]) # Apply eigenpair
-        r = operator.solution()              # Get residual
+        r = operator.solution(0)             # Get residual
         norm[i] = np.linalg.norm(r)          # Save the norm
 
     # Done
@@ -174,18 +174,15 @@ def multiSolve(solver, B, w):
     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)
+    # Solve
+    solver.solve(B, w)
 
-    # Loop and solve
+    # Get solutions
+    nSolve = B.shape[1]
+    size   = solver.size()
+    X      = np.matlib.empty((size, nSolve), dtype=complex)
     for i in range(nSolve):
-        b = B[:, i]
-        solver.solve(b, w)
-        X[:, i] = solver.solution()
+        X[:, i] = solver.solution(i)
 
     # Done
     return X

bin/solver.py

@@ -18,11 +18,13 @@ class GetDPWave:
     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?
-    doPostpro -- should only a view be created for x?
-    doApply -- should only the application of x be computed?
+    nRHS -- the number of right hand side that should be considered (default: 1)
+    b_i / b~{i} -- the ith right hand side (RHS) to use (first index is 0)
+    x_i / x~{i} -- the ith solution vector computed by GetDP (first index is 0)
+    doInit -- should some initialization be done (default: 0)?
+    doSolve -- should Ax_i = b_i be solved for all i (default: 0)?
+    doPostpro -- should only a view be created for x (default: 0)?
+    doApply -- should only the application of x be computed (default: 0)?
     fileName -- post-processing file name
     """
 
@@ -44,19 +46,20 @@ class GetDPWave:
         self.__resolution = resolution
         self.__optional   = optional
 
-        # Generate DoFs and initialise RHS and solution vectors in GetDP
+        # Generate DoFs and assemble a first system
+        GetDP.GetDPSetNumber("doInit", 1);
         GetDP.GetDP(["getdp",     self.__pro,
                      "-msh",      self.__mesh,
                      "-solve",    self.__resolution,
                      "-v", "2"] + self.__optional)
 
-    def solution(self):
-        """Returns the solution"""
-        return self.__toNumpy(GetDP.GetDPGetNumber("x"))
+    def solution(self, i):
+        """Returns the solution for the ith RHS (first index is 0)"""
+        return self.__toNumpy(GetDP.GetDPGetNumber("x_" + str(i)))
 
     def size(self):
         """Returns the number of degrees of freedom"""
-        return self.solution().shape[0]
+        return self.solution(0).shape[0]
 
     def apply(self, x, w):
         """Applies x to the operator with a pulsation of w
@@ -65,23 +68,34 @@ class GetDPWave:
         """
         GetDP.GetDPSetNumber("angularFreqRe", np.real(w).tolist())
         GetDP.GetDPSetNumber("angularFreqIm", np.imag(w).tolist())
-        GetDP.GetDPSetNumber("x", self.__toGetDP(x))
+        GetDP.GetDPSetNumber("x_0", self.__toGetDP(x))
         GetDP.GetDPSetNumber("doApply", 1)
         GetDP.GetDP(["getdp",  self.__pro,
                      "-msh",   self.__mesh,
                      "-cal"] + self.__optional)
-        GetDP.GetDPSetNumber("doApply", 0)
 
     def solve(self, b, w):
-        """Solves with b as RHS and w as complex angular frequency"""
+        """Solves with b as RHS and w as complex angular frequency
+
+        b is a matrix and each column is a different RHS
+        """
+        # Number of RHS
+        nRHS = b.shape[1]
+
+        # Set variables
+        GetDP.GetDPSetNumber("nRHS", nRHS)
         GetDP.GetDPSetNumber("angularFreqRe", np.real(w).tolist())
         GetDP.GetDPSetNumber("angularFreqIm", np.imag(w).tolist())
-        GetDP.GetDPSetNumber("b", self.__toGetDP(b))
-        GetDP.GetDPSetNumber("imposeRHS", 1)
+
+        # Set RHS
+        for i in range(nRHS):
+            GetDP.GetDPSetNumber("b_" + str(i), self.__toGetDP(b[:, i]))
+
+        # Solve
+        GetDP.GetDPSetNumber("doSolve", 1)
         GetDP.GetDP(["getdp",  self.__pro,
                      "-msh",   self.__mesh,
                      "-cal"] + self.__optional)
-        GetDP.GetDPSetNumber("imposeRHS", 0)
 
     def view(self, x, fileName):
         """Generates a post-processing view
@@ -92,13 +106,12 @@ class GetDPWave:
 
         This method generates a linear system
         """
-        GetDP.GetDPSetNumber("x", self.__toGetDP(x))
+        GetDP.GetDPSetNumber("x_0", self.__toGetDP(x))
         GetDP.GetDPSetNumber("doPostpro", 1)
         GetDP.GetDPSetString("fileName", fileName)
         GetDP.GetDP(["getdp",  self.__pro,
                      "-msh",   self.__mesh,
                      "-cal"] + self.__optional)
-        GetDP.GetDPSetNumber("doPostpro", 0)
 
     @staticmethod
     def __toNumpy(vGetDP):

example/maxwell_linear/maxwell.pro

@@ -16,11 +16,15 @@ Function{
          angularFreqRe, angularFreqIm];
 
   // Algebraic data //
-  DefineConstant[x() = {},  // Solution
-                 b() = {}]; // Right hand side
+  DefineConstant[nRHS = 1];       // Number of RHS for this run
+  For i In {0:nRHS-1}
+    DefineConstant[x~{i}() = {},  // Solution
+                   b~{i}() = {}]; // Right hand side
+  EndFor
 
   // Control data //
-  DefineConstant[imposeRHS = 0,          // Should I use an imposed RHS?
+  DefineConstant[doInit    = 0,          // Should I initialize some stuff?
+                 doSolve   = 0,          // Should I solve Ax = b?
                  doPostpro = 0,          // Should I only create a view for x()?
                  doApply   = 0,          // Should I only apply x(): x <- Ax?
                  fileName  = "eig.pos"]; // Postpro file name
@@ -91,21 +95,32 @@ Resolution{
     Operation{
       Generate[A];
 
-      If(imposeRHS)
-        CopyRightHandSide[b(), A];
+      If(doInit)
+        CopySolution[A, x~{0}()];
       EndIf
 
-      If(!doPostpro && !doApply)
+      If(doSolve)
+        // Full solve for first RHS
+        CopyRightHandSide[b~{0}(), A];
         Solve[A];
-        CopySolution[A, x()];
+        CopySolution[A, x~{0}()];
+
+        // SolveAgain for other RHSs
+        For i In {1:nRHS-1}
+          CopyRightHandSide[b~{i}(), A];
+          SolveAgain[A];
+          CopySolution[A, x~{i}()];
+        EndFor
       EndIf
+
       If(doApply)
-        CopySolution[x(), A];
+        CopySolution[x~{0}(), A];
         Apply[A];
-        CopySolution[A, x()];
+        CopySolution[A, x~{0}()];
       EndIf
+
       If(doPostpro)
-        CopySolution[x(), A];
+        CopySolution[x~{0}(), A];
         PostOperation[Maxwell];
       EndIf
     }

example/maxwell_sibc/cimParameters.pro

@@ -11,11 +11,15 @@ Function{
                     angularFreqIm];  // [rad/m]
 
   // Algebraic data //
-  DefineConstant[x() = {},  // Solution
-                 b() = {}]; // Right hand side
+  DefineConstant[nRHS = 1];       // Number of RHS for this run
+  For i In {0:nRHS-1}
+    DefineConstant[x~{i}() = {},  // Solution
+                   b~{i}() = {}]; // Right hand side
+  EndFor
 
   // Control data //
-  DefineConstant[imposeRHS = 0,          // Should I use an imposed RHS?
+  DefineConstant[doInit    = 0,          // Should I initialize some stuff?
+                 doSolve   = 0,          // Should I solve Ax = b?
                  doPostpro = 0,          // Should I only create a view for x()?
                  doApply   = 0,          // Should I only apply x(): x <- Ax?
                  fileName  = "eig.pos"]; // Postpro file name

example/maxwell_sibc/cimResolution.pro

@@ -9,23 +9,32 @@ Resolution{
     Operation{
       Generate[A];
 
-      If(imposeRHS)
-        CopyRightHandSide[b(), A];
+      If(doInit)
+        CopySolution[A, x~{0}()];
       EndIf
 
-      If(!doPostpro && !doApply)
+      If(doSolve)
+        // Full solve for first RHS
+        CopyRightHandSide[b~{0}(), A];
         Solve[A];
-        CopySolution[A, x()];
+        CopySolution[A, x~{0}()];
+
+        // SolveAgain for other RHSs
+        For i In {1:nRHS-1}
+          CopyRightHandSide[b~{i}(), A];
+          SolveAgain[A];
+          CopySolution[A, x~{i}()];
+        EndFor
       EndIf
 
       If(doApply)
-        CopySolution[x(), A];
+        CopySolution[x~{0}(), A];
         Apply[A];
-        CopySolution[A, x()];
+        CopySolution[A, x~{0}()];
       EndIf
 
       If(doPostpro)
-        CopySolution[x(), A];
+        CopySolution[x~{0}(), A];
         PostOperation[Post];
       EndIf
     }

example/maxwell_sibc_divfree/Makefile

-BIN=../../bin
-
-sphere: init cim
-
-init:
-	gmsh sphere.geo -3
-
-cim:
-	python ${BIN}/cim.py sphere.pro sphere.msh Solve 8.22543e9+2.46361e1j 1e8 -nodes 10 -lStart 4
-
-clean:
-	rm -f *.pos
-	rm -f *.msh
-	rm -f *.pre
-	rm -f *.res

example/maxwell_sibc_divfree/README.md

-Same example as [maxwell\_sibc](../maxwell_sibc), but imposes a divergence free
-source for Beyn. This is achieved by taking the curl of the random Beyn source.
-
-This example requires a [GetDP](http://getdp.info) version compiled with [Gmsh](
-http://gmsh.info) to handle Fields.

example/maxwell_sibc_divfree/cimParameters.pro

-/////////////////////////////////////////////////////////////////
-// Parameters for contour integral method.                     //
-// The master file should inclue "cimResolution.pro" later on. //
-/////////////////////////////////////////////////////////////////
-Function{
-  // Physical data //
-  DefineConstant[angularFreqRe = 1,  // [rad/s]
-                 angularFreqIm = 0]; // [rad/s]
-
-  Puls[]  = Complex[angularFreqRe,
-                    angularFreqIm];  // [rad/m]
-
-  // Algebraic data //
-  DefineConstant[x() = {},  // Solution
-                 b() = {}]; // Right hand side
-
-  // Control data //
-  DefineConstant[imposeRHS = 0,          // Should I use an imposed RHS?
-                 doPostpro = 0,          // Should I only create a view for x()?
-                 doApply   = 0,          // Should I only apply x(): x <- Ax?
-                 fileName  = "eig.pos"]; // Postpro file name
-}

example/maxwell_sibc_divfree/cimResolution.pro

-/////////////////////////////////////////////////////////////////////
-// Resolution for contour integral method.                         //
-// Transforms Beyn random source into a divergence free one.       //
-// Warning, the master file should include "cimParameters.pro",    //
-// and should also define "FormulationCIM" the formulation to use. //
-/////////////////////////////////////////////////////////////////////
-Resolution{
-  { Name Solve;
-    System{
-      { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; }
-      { Name B; NameOfFormulation Source;         Type ComplexValue; }
-    }
-    Operation{
-      // Divergence free source
-      If(imposeRHS)
-        Generate[B];
-        CopyRightHandSide[b(), B];
-        Solve[B];
-        PostOperation[Source];
-      EndIf
-
-      // Beyn (source given by dSrc)
-      Generate[A];
-
-      If(!doPostpro && !doApply)
-        Solve[A];
-        CopySolution[A, x()];
-      EndIf
-
-      If(doApply)
-        CopySolution[x(), A];
-        Apply[A];
-        CopySolution[A, x()];
-      EndIf
-
-      If(doPostpro)
-        CopySolution[x(), A];
-        PostOperation[Post];
-      EndIf
-    }
-  }
-}

example/maxwell_sibc_divfree/sphere.dat

-// Geometry //
-DefineConstant[R  = { 100e-3, Name "Input/00Geometry/00Radius [m]" }];
-DefineConstant[K1 = { 1/R,    Name "Input/00Geometry/01Curvature 1 [m^-1]",
-                              ReadOnly 1}];
-DefineConstant[K2 = { 1/R,    Name "Input/00Geometry/02Curvature 2 [m^-1]",
-                              ReadOnly 1}];
-
-// Mesh //
-DefineConstant[md = { 10,     Name "Input/01Mesh/00Density [Radius^-1]" }];
-cl = R/md;
-
-// Material //
-DefineConstant[Sigma = { 1e15,
-                         Name "Input/02Material/00Conductivity [S m^-1]" }];

example/maxwell_sibc_divfree/sphere.geo

-Include "sphere.dat";
-
-Point(0) = { 0,  0,  0, cl};
-
-Point(1) = {+R,  0,  0, cl};
-Point(2) = { 0, +R,  0, cl};
-Point(3) = {-R,  0,  0, cl};
-Point(4) = { 0, -R,  0, cl};
-Point(5) = { 0,  0, +R, cl};
-Point(6) = { 0,  0, -R, cl};
-
-Circle(01) = {1, 0, 2};
-Circle(02) = {2, 0, 3};
-Circle(03) = {3, 0, 4};
-Circle(04) = {4, 0, 1};
-
-Circle(05) = {5, 0, 1};
-Circle(06) = {1, 0, 6};
-Circle(07) = {6, 0, 3};
-Circle(08) = {3, 0, 5};
-
-Circle(09) = {2, 0, 6};
-Circle(10) = {6, 0, 4};
-Circle(11) = {4, 0, 5};
-Circle(12) = {5, 0, 2};
-
-Line Loop(1) = {1, -12, 5};
-Line Loop(2) = {1, 9, -6};
-Line Loop(3) = {9, 7, -2};
-Line Loop(4) = {2, 8, 12};
-Line Loop(5) = {11, 5, -4};
-Line Loop(6) = {4, 6, 10};
-Line Loop(7) = {10, -3, -7};
-Line Loop(8) = {3, 11, -8};
-
-Ruled Surface(1) = {1};
-Ruled Surface(2) = {2};
-Ruled Surface(3) = {3};
-Ruled Surface(4) = {4};
-Ruled Surface(5) = {5};
-Ruled Surface(6) = {6};
-Ruled Surface(7) = {7};
-Ruled Surface(8) = {8};
-
-Surface Loop(1) = {3, 2, 1, 4, 8, 7, 6, 5};
-Volume(1) = {1};
-
-Physical  Volume(1) = {1};
-Physical Surface(2) = {1, 2, 3, 4, 5, 6, 7, 8};
-Physical   Point(3) = {0};

example/maxwell_sibc_divfree/sphere.pro

-Include "sphere.dat";
-
-Group{
-  Omega = Region[1];
-  Bndry = Region[2];
-}
-
-Include "./cimParameters.pro"; // Functions for cim.py
-
-Function{
-  // Imaginary Unit //
-  J[] = Complex[0, 1];
-
-  // Physical Constants //
-  C0   = 299792458;      // [m/s]
-  Mu0  = 4*Pi*1e-7;      // [H/m]
-  Eps0 = 1/(Mu0 * C0^2); // [F/m]
-
-  // Conductivity (coming from sphere.dat) //
-  Sigma[] = Sigma; // [S/m]
-
-  // Leontovich SIBC (H-Formulation) //
-  SL[] = Sqrt[(J[] * Re[Puls[]] * Mu0) / Sigma[]]; // H-Formulation
-  SIbc[] = -1.0 / SL[];                            // E-Formulation
-
-  // Curl of Beyn source //
-  dSrc[] = ComplexVectorField[XYZ[]]{1};
-}
-
-Integration{
-  { Name I1;
-    Case{
-      { Type Gauss;
-        Case{
-          { GeoElement Tetrahedron; NumberOfPoints 4; } // O1 * O1
-          { GeoElement Triangle;    NumberOfPoints 3; } // O1 * O1
-          { GeoElement Line;        NumberOfPoints 2; } // O1 * O1
-        }
-      }
-    }
-  }
-}
-
-Jacobian{
-  { Name Jac;
-    Case{
-      { Region Omega; Jacobian Vol; }
-      { Region Bndry; Jacobian Sur; }
-    }
-  }
-}
-
-FunctionSpace{
-  { Name HCurl; Type Form1;
-    BasisFunction{
-      { Name sx; NameOfCoef ux; Function BF_Edge;
-        Support Region[{Omega, Bndry}]; Entity EdgesOf[All]; }
-    }
-  }
-}
-
-Formulation{
-  { Name Source; Type FemEquation;
-    Quantity{
-      { Name s; Type Local; NameOfSpace HCurl; }
-    }
-    Equation{
-      Galerkin{ [Dof{s}, {s}];
-        In Omega; Jacobian Jac; Integration I1; }
-    }
-  }
-  { Name FormulationCIM; Type FemEquation;
-    Quantity{
-      { Name e; Type Local; NameOfSpace HCurl; }
-    }
-    Equation{
-      // Maxwell
-      Galerkin{ [           1/Mu0 * Dof{d e}, {d e}];
-        In Omega; Jacobian Jac; Integration I1; }
-      Galerkin{ [-Puls[]^2 * Eps0 * Dof{  e}, {  e}];
-        In Omega; Jacobian Jac; Integration I1; }
-
-      // IBC
-      Galerkin{ [-J[] * Puls[] * SIbc[] * Dof{e}, {e}];
-        In Bndry; Jacobian Jac; Integration I1; }
-
-      // Source
-      If(imposeRHS)
-        Galerkin{ [dSrc[], {e}];
-          In Omega; Jacobian Jac; Integration I1; }
-      EndIf
-    }
-  }
-}
-
-Include "./cimResolution.pro"; // Resolution for cim.py
-
-PostProcessing{
-  { Name Source; NameOfFormulation Source;
-    Quantity{
-      { Name dSrc; Value{ Term{ [{d s}]; In Omega; Jacobian Jac; } } }
-    }
-  }
-  { Name Post;   NameOfFormulation FormulationCIM;
-    Quantity{
-      { Name E;    Value{ Term{ [{  e}]; In Omega; Jacobian Jac; } } }
-    }
-  }
-}
-
-PostOperation{
-  { Name Source; NameOfPostProcessing Source;
-    Operation{
-      Print[dSrc, OnElementsOf Omega, StoreInField 1];
-    }
-  }
-  { Name Post;   NameOfPostProcessing Post;
-    Operation{
-      Print[E,    OnElementsOf Omega, File fileName];
-    }
-  }
-}