added ml, splu impl for pract. ML as default

integrators/_details/interfaces/tpsPreconditioner.py

@@ -108,6 +108,7 @@ class TPSPreconditioner(_preconditionerInterface._PreconditionerInterface):
     import scipy.sparse.linalg
     import scipy.sparse as sp
     import numpy as np
+    import pyamg
 
     if(self.parameters.alphaPreconditioner == None):
       A = self._A_stat[:self._A_stat.shape[0]//3, :self._A_stat.shape[0]//3]
@@ -115,16 +116,34 @@ class TPSPreconditioner(_preconditionerInterface._PreconditionerInterface):
       A = self._A_stat_precond[:self._A_stat_precond.shape[0]//3 \
         , :self._A_stat_precond.shape[0]//3]
 
-    
     self._Compute_Q()
-    metaPre = 1.0 / A.diagonal()
-    metaP = lambda x: metaPre*x
-    P = scipy.sparse.linalg.LinearOperator(A.shape, matvec=metaP)
-    self._prePractical = lambda y: np.concatenate(( \
-      scipy.sparse.linalg.cg(A, y[:y.size//3], tol=self._solvetol/10, M=P)[0], \
-      scipy.sparse.linalg.cg(A, y[y.size//3:2*y.size//3], tol=self._solvetol/10, M=P)[0], \
-      scipy.sparse.linalg.cg(A, y[2*y.size//3:], tol=self._solvetol/10, M=P)[0]))
 
+    # trying ML, SPLU, CG preconditioner, ML as default
+    # (CG not working well for h << lex)
+    # it seems ML faster, but SPLU less iterations
+    # way faster than true inverse, especially for large NV
+    opt = "ml" # "ml", "splu"
+    if opt == "ml":
+      ml = pyamg.ruge_stuben_solver(A)
+      pre = ml.aspreconditioner()
+      
+      self._prePractical = lambda y: np.concatenate(( \
+        (pre * y[:y.size//3], pre * y[y.size//3:2*y.size//3], pre * y[2*y.size//3:])))
+    elif opt == "splu":
+      cgs = pyamg.coarse_grid_solver("splu")
+      self._prePractical = lambda y: np.concatenate(( \
+        (cgs(A, y[:y.size//3]), \
+        cgs(A, y[y.size//3:2*y.size//3]), \
+        cgs(A, y[2*y.size//3:]))))
+    elif opt == "cg":
+      metaPre = 1.0 / A.diagonal()
+      metaP = lambda x: metaPre*x
+      P = scipy.sparse.linalg.LinearOperator(A.shape, matvec=metaP)
+      self._prePractical = lambda y: np.concatenate(( \
+        scipy.sparse.linalg.cg(A, y[:y.size//3], tol=self._solvetol/10, M=P)[0], \
+        scipy.sparse.linalg.cg(A, y[y.size//3:2*y.size//3], tol=self._solvetol/10, M=P)[0], \
+        scipy.sparse.linalg.cg(A, y[2*y.size//3:], tol=self._solvetol/10, M=P)[0]))
+  
     practical2D = lambda x: \
       self._Q.dot(self._prePractical(self._Q.transpose().dot(x)))
 
@@ -147,9 +166,11 @@ class TPSPreconditioner(_preconditionerInterface._PreconditionerInterface):
       A = self._A_stat_precond[:self._A_stat_precond.shape[0]//3 \
         , :self._A_stat_precond.shape[0]//3]
 
-    # trying AMG, SPLU, CG preconditioner, CG as default. 
+    # trying ML, SPLU, CG preconditioner, ML as default
+    # (CG not working well for h << lex)
+    # it seems ML faster, but SPLU less iterations
     # way faster than true inverse, especially for large NV
-    opt = "cg" # "ml", "splu"
+    opt = "ml" # "ml", "splu"
     if opt == "ml":
       ml = pyamg.ruge_stuben_solver(A)
       pre = ml.aspreconditioner()