From d255620ca24d9e9eace3faea50232afa19e63644 Mon Sep 17 00:00:00 2001
From: Carl-Martin Pfeiler <carl-martin.pfeiler@asc.tuwien.ac.at>
Date: Mon, 19 Aug 2019 15:48:04 +0200
Subject: [PATCH] Jacobi Preconditioner for MPS

---
 integrators/_details/interfaces/__init__.py   |  4 +-
 .../_details/interfaces/mpsPreconditioner.py  | 72 +++++++++++++++++++
 integrators/mps.py                            |  9 ++-
 3 files changed, 82 insertions(+), 3 deletions(-)
 create mode 100644 integrators/_details/interfaces/mpsPreconditioner.py

diff --git a/integrators/_details/interfaces/__init__.py b/integrators/_details/interfaces/__init__.py
index cb04408..fea2ff1 100644
--- a/integrators/_details/interfaces/__init__.py
+++ b/integrators/_details/interfaces/__init__.py
@@ -8,6 +8,8 @@ from commics.integrators._details.interfaces.parameterInterface import Parameter
 
 from commics.integrators._details.interfaces.tpsPreconditioner import TPSPreconditioner
 from commics.integrators._details.interfaces.tpsPreconditionerJacobiStationary import TPSPreconditionerJacobiStationary
+from commics.integrators._details.interfaces.mpsPreconditioner import MPSPreconditioner
 
 __all__ = ["CommunicationInterface", "ComputeInterface", "ContributionInterface", "WriteInterface", "ParameterInterface" \
-  , "TPSPreconditioner", "TPSPreconditionerJacobiStationary"]
+  , "TPSPreconditioner", "TPSPreconditionerJacobiStationary" \
+  , "MPSPreconditioner"]
diff --git a/integrators/_details/interfaces/mpsPreconditioner.py b/integrators/_details/interfaces/mpsPreconditioner.py
new file mode 100644
index 0000000..4382f38
--- /dev/null
+++ b/integrators/_details/interfaces/mpsPreconditioner.py
@@ -0,0 +1,72 @@
+from . import _preconditionerInterface
+
+################################################################################
+
+
+class MPSPreconditioner(_preconditionerInterface._PreconditionerInterface):
+
+
+#------------------------------------------------------------------------------#
+
+  
+  def _SetUpPreconditioner(self):
+
+    if(self.parameters.preconditionerType == "jacobi"):
+      self._SetUpPreconditioner_Jacobi()
+    elif(self.parameters.preconditionerType == None):
+      self._SetUpPreconditioner_None()
+    else:
+      raise ValueError("preconditionerType = " \
+      , self.parameters.preconditionerType, "not known")
+
+
+#------------------------------------------------------------------------------#
+
+
+  def _UpdatePreconditioner(self):
+    if(self.parameters.preconditionerType == "jacobi"):
+      self._UpdatePreconditioner_Jacobi()
+    elif(self.parameters.preconditionerType == None):
+      self._UpdatePreconditioner_None()
+    else:
+      raise ValueError("preconditionerType = ", self.parameters.preconditionerType \
+        , "not known")
+
+
+#------------------------------------------------------------------------------#
+
+  
+  def _SetUpPreconditioner_Jacobi(self):
+    import scipy.sparse.linalg
+    pjac = 1.0 / self._LHS_stat.diagonal()
+
+    jacobi3D = lambda x: pjac * x
+    self._preconditioner = scipy.sparse.linalg.LinearOperator( \
+      (pjac.size, pjac.size), matvec=jacobi3D)
+
+
+#------------------------------------------------------------------------------#
+
+  
+  def _SetUpPreconditioner_None(self):
+    self._preconditioner = None
+
+
+#------------------------------------------------------------------------------#
+
+
+  def _UpdatePreconditioner_Jacobi(self):
+    pass
+
+
+#------------------------------------------------------------------------------#
+
+
+  def _UpdatePreconditioner_None(self):
+    pass
+
+
+#------------------------------------------------------------------------------#
+
+
+################################################################################
diff --git a/integrators/mps.py b/integrators/mps.py
index b3a62cc..98b0701 100644
--- a/integrators/mps.py
+++ b/integrators/mps.py
@@ -1,11 +1,13 @@
 ### MPS ###
 from . import _integrator
 from . import _forms
+from . import _details
 
 ################################################################################
 
 
 class MPS(_integrator._Integrator, \
+  _details.interfaces.mpsPreconditioner.MPSPreconditioner, \
   _forms.blf_common.BilinearForms_Common, \
   _forms.lf_common.LinearForms_Common, \
   _forms.lf_eta.LinearForms_Eta):
@@ -39,6 +41,8 @@ class MPS(_integrator._Integrator, \
     _integrator._Integrator.GetReady(self)
     self._Mag.Project()
 
+    self._SetUpPreconditioner()
+
 
 #------------------------------------------------------------------------------#
 
@@ -254,7 +258,8 @@ class MPS(_integrator._Integrator, \
 #------------------------------------------------------------------------------#
 
 
-  def _StepPrepareWithoutTaskManager(self): pass
+  def _StepPrepareWithoutTaskManager(self):
+    self._UpdatePreconditioner()
 
 
 #------------------------------------------------------------------------------#
@@ -337,7 +342,7 @@ class MPS(_integrator._Integrator, \
 
       self._Mag_eta.m.vec.FV().NumPy()[:], succ \
         = scipy.sparse.linalg.gmres(LHS, RHS , x0=None \
-        , tol=self._solvetol, maxiter=4000, M=None, callback=self._callback)
+        , tol=self._solvetol, maxiter=4000, M=self._preconditioner, callback=self._callback)
       
       diff_gf.vec.FV().NumPy()[:] = self._Mag_eta.m.vec.FV().NumPy()-oldIterate
       itError = sqrt( Integrate(diff_cf*diff_cf, self._GetMesh()) )
-- 
GitLab