diff --git a/integrators/_details/interfaces/tpsPreconditioner.py b/integrators/_details/interfaces/tpsPreconditioner.py
index bb12755cb5ebe88798d3d7ab9d119dbc9b8c5f5d..c8df2a7b7d21cfee2deaf1dd27fcf32bc21b0e35 100644
--- a/integrators/_details/interfaces/tpsPreconditioner.py
+++ b/integrators/_details/interfaces/tpsPreconditioner.py
@@ -118,23 +118,10 @@ class TPSPreconditioner(_preconditionerInterface._PreconditionerInterface):
 
     self._Compute_Q()
 
-    # 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()
-      pre1D = lambda z: pre * z
-    elif opt == "splu":
-      cgs = pyamg.coarse_grid_solver("splu")
-      pre1D = lambda z: cgs(A, z)
-    elif opt == "cg":
-      metaPre = 1.0 / A.diagonal()
-      metaP = lambda x: metaPre*x
-      P = scipy.sparse.linalg.LinearOperator(A.shape, matvec=metaP)
-      pre1D = lambda z: scipy.sparse.linalg.cg(A, z, tol=self._solvetol/10, M=P)[0]
+    # realize "inverse" of stationary part via multilevel solver
+    ml = pyamg.ruge_stuben_solver(A)
+    pre = ml.aspreconditioner()
+    pre1D = lambda z: pre * z
   
     self._prePractical = lambda y: np.concatenate((
       pre1D(y[:y.size//3]), pre1D(y[y.size//3:2*y.size//3]), pre1D(y[2*y.size//3:])))
@@ -161,23 +148,10 @@ class TPSPreconditioner(_preconditionerInterface._PreconditionerInterface):
       A = self._A_stat_precond[:self._A_stat_precond.shape[0]//3 \
         , :self._A_stat_precond.shape[0]//3]
 
-    # 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()
-      pre1D = lambda z: pre * z
-    elif opt == "splu":
-      cgs = pyamg.coarse_grid_solver("splu")
-      pre1D = lambda z: cgs(A, z)
-    elif opt == "cg":
-      metaPre = 1.0 / A.diagonal()
-      MetaPre = lambda x: metaPre * x
-      P = scipy.sparse.linalg.LinearOperator(A.shape, matvec=MetaPre)
-      pre1D = lambda z: scipy.sparse.linalg.cg(A, z, tol=self._solvetol/10, M=P)[0]
+    # realize "inverse" of stationary part via multilevel solver
+    ml = pyamg.ruge_stuben_solver(A)
+    pre = ml.aspreconditioner()
+    pre1D = lambda z: pre * z
 
     stationary2D = lambda x: \
       np.vstack((pre1D(x[0::2]), pre1D(x[1::2]))).reshape((-1,),order='F')