first SDC run works

Author: pancetta

pySDC/implementations/datatype_classes/petsc_vec.py

@@ -0,0 +1,106 @@
+from petsc4py import PETSc
+
+from pySDC.core.Errors import DataError
+
+
+class petsc_vec(PETSc.Vec):
+    __array_priority__ = 1000  # otherwise rmul with float64 does not work (don't ask, won't tell)
+
+    def __new__(cls, init=None, val=0.0):
+        if isinstance(init, petsc_vec) or isinstance(init, PETSc.Vec):
+            obj = PETSc.Vec().__new__(cls)
+            init.copy(obj)
+        elif isinstance(init, PETSc.DMDA):
+            tmp = init.createGlobalVector()
+            obj = petsc_vec(tmp)
+            objarr = init.getVecArray(obj)
+            objarr[:] = val
+        else:
+            obj = PETSc.Vec().__new__(cls)
+        return obj
+
+    def __abs__(self):
+        """
+        Overloading the abs operator
+
+        Returns:
+            float: absolute maximum of all vec values
+        """
+        # take absolute values of the mesh values (INF = 3)
+        return self.norm(3)
+
+    def isend(self, dest=None, tag=None, comm=None):
+        """
+        Routine for sending data forward in time (non-blocking)
+
+        Args:
+            dest (int): target rank
+            tag (int): communication tag
+            comm: communicator
+
+        Returns:
+            request handle
+        """
+        return comm.Issend(self[:], dest=dest, tag=tag)
+
+    def irecv(self, source=None, tag=None, comm=None):
+        """
+        Routine for receiving in time
+
+        Args:
+            source (int): source rank
+            tag (int): communication tag
+            comm: communicator
+
+        Returns:
+            None
+        """
+        return comm.Irecv(self[:], source=source, tag=tag)
+
+    def bcast(self, root=None, comm=None):
+        """
+        Routine for broadcasting values
+
+        Args:
+            root (int): process with value to broadcast
+            comm: communicator
+
+        Returns:
+            broadcasted values
+        """
+        comm.Bcast(self[:], root=root)
+        return self
+
+
+class petsc_vec_imex(object):
+    """
+    RHS data type for Vec with implicit and explicit components
+
+    This data type can be used to have RHS with 2 components (here implicit and explicit)
+
+    Attributes:
+        impl (petsc_vec): implicit part
+        expl (petsc_vec): explicit part
+    """
+
+    def __init__(self, init, val=0.0):
+        """
+        Initialization routine
+
+        Args:
+            init: can either be a tuple (one int per dimension) or a number (if only one dimension is requested)
+                  or another imex_mesh object
+            val (float): an initial number (default: 0.0)
+        Raises:
+            DataError: if init is none of the types above
+        """
+
+        if isinstance(init, type(self)):
+            self.impl = petsc_vec(init.impl)
+            self.expl = petsc_vec(init.expl)
+        elif isinstance(init, PETSc.DMDA):
+            self.impl = petsc_vec(init, val=val)
+            self.expl = petsc_vec(init, val=val)
+        # something is wrong, if none of the ones above hit
+        else:
+            raise DataError('something went wrong during %s initialization' % type(self))

pySDC/implementations/problem_classes/HeatEquation_2D_PETSc_forced.py

@@ -3,7 +3,7 @@
 
 from pySDC.core.Errors import ParameterError, ProblemError
 from pySDC.core.Problem import ptype
-from pySDC.implementations.datatype_classes.petsc_dmda_grid import petsc_data, rhs_imex_petsc_data
+from pySDC.implementations.datatype_classes.petsc_vec import petsc_vec, petsc_vec_imex
 
 
 # noinspection PyUnusedLocal
@@ -20,7 +20,7 @@ class heat2d_petsc_forced(ptype):
         ksp: PETSc linear solver object
     """
 
-    def __init__(self, problem_params, dtype_u=petsc_data, dtype_f=rhs_imex_petsc_data):
+    def __init__(self, problem_params, dtype_u=petsc_vec, dtype_f=petsc_vec_imex):
         """
         Initialization routine
 
@@ -165,10 +165,10 @@ def eval_f(self, u, t):
 
         f = self.dtype_f(self.init)
         # evaluate Au for implicit part
-        self.A.mult(u.values, f.impl.values)
+        self.A.mult(u, f.impl)
 
         # evaluate forcing term for explicit part
-        fa = self.init.getVecArray(f.expl.values)
+        fa = self.init.getVecArray(f.expl)
         xv, yv = np.meshgrid(range(self.xs, self.xe), range(self.ys, self.ye), indexing='ij')
         fa[self.xs:self.xe, self.ys:self.ye] = -np.sin(np.pi * self.params.freq * xv * self.dx) * \
             np.sin(np.pi * self.params.freq * yv * self.dy) * \
@@ -192,10 +192,9 @@ def solve_system(self, rhs, factor, u0, t):
 
         me = self.dtype_u(u0)
         self.ksp.setOperators(self.Id - factor * self.A)
-        self.ksp.solve(rhs.values, me.values)
+        self.ksp.solve(rhs, me)
         self.ksp_ncalls += 1
         self.ksp_itercount += int(self.ksp.getIterationNumber())
-
         return me
 
     def u_exact(self, t):
@@ -210,7 +209,7 @@ def u_exact(self, t):
         """
 
         me = self.dtype_u(self.init)
-        xa = self.init.getVecArray(me.values)
+        xa = self.init.getVecArray(me)
         for i in range(self.xs, self.xe):
             for j in range(self.ys, self.ye):
                 xa[i, j] = np.sin(np.pi * self.params.freq * i * self.dx) * \

pySDC/playgrounds/other/petsc_datatype.py

@@ -1,25 +1,39 @@
 from petsc4py import PETSc
-
+import numpy as np
 
 class mymesh(PETSc.Vec):
-
+    __array_priority__ = 1000
     def __new__(cls, init=None, val=0.0):
-        if isinstance(init, mymesh):
+        if isinstance(init, mymesh) or isinstance(init, PETSc.Vec):
             obj = PETSc.Vec().__new__(cls)
             init.copy(obj)
-            # obj.
-            # obj[:] = init[:]
-            # obj.part1 = obj[0, :]
-            # obj.part2 = obj[1, :]
-        elif isinstance(init, tuple):
-            obj = PETSc.Vec().__new__(cls)
-            obj.createMPI(size=init[0], comm=init[1])
-            objarr = obj.getArray()
+        elif isinstance(init, PETSc.DMDA):
+            tmp = init.createGlobalVector()
+            obj = mymesh(tmp)
+            objarr = init.getVecArray(obj)
             objarr[:] = val
         else:
             obj = PETSc.Vec().__new__(cls)
         return obj
 
+    # def __rmul__(self, other):
+    #     print('here')
+    #     tmp = self.getArray()
+    #     tmp[:] *= other
+    #     return self
+
+    def __abs__(self):
+        """
+        Overloading the abs operator
+
+        Returns:
+            float: absolute maximum of all mesh values
+        """
+        # take absolute values of the mesh values
+        # print(self.norm(3))
+        return self.norm(3)
+
+
     # def __array_finalize__(self, obj):
     #     """
     #     Finalizing the datatype. Without this, new datatypes do not 'inherit' the communicator.
@@ -29,39 +43,49 @@ def __new__(cls, init=None, val=0.0):
     #     self.part1 = getattr(obj, 'part1', None)
     #     self.part2 = getattr(obj, 'part2', None)
 
-u = mymesh((10, PETSc.COMM_WORLD), val=9)
-uarr = u.getArray()
-# uarr[:] = 1
-print(uarr, u.getOwnershipRange())
+
+# u = mymesh((16, PETSc.COMM_WORLD), val=9)
+# uarr = u.getArray()
+# # uarr[:] = np.random.rand(4,4)
+# print(uarr, u.getOwnershipRanges())
+# v = mymesh(u)
+# varr = v.getArray()
+# print(varr, v.getOwnershipRange() )
+# if v.comm.getRank() == 0:
+#     uarr[0] = 7
+#
+# print(uarr)
+# print(varr)
+#
+# w = u + v
+# warr = w.getArray()
+# print(warr, w.getOwnershipRange())
+
+da = PETSc.DMDA().create([4, 4], stencil_width=1)
+
+u = mymesh(da, val=9.0)
+uarr = da.getVecArray(u)
+
 v = mymesh(u)
-varr = v.getArray()
-print(varr, v.getOwnershipRange() )
+varr = da.getVecArray(v)
 if v.comm.getRank() == 0:
-    uarr[0] = 7
+    uarr[0, 0] = 7
+print(uarr[:], uarr.shape, type(u))
+print(varr[:], v.getOwnershipRange(), type(v))
 
-print(uarr)
-print(varr)
+w = np.float(1.0) * (u + v)
+warr = da.getVecArray(w)
+print(warr[:], da.getRanges(), type(w))
 
-w = u + v
-warr = w.getArray()
-print(warr, w.getOwnershipRange())
-exit()
+print(w.norm(3))
+print(abs(w))
 
+# v = PETSc.Vec().createMPI(16)
+# varr = v.getArray()
+# varr[:] = 9.0
+# a = np.float64(1.0) * v
+# print(type(a))
 
-m = mymesh((10, np.dtype('float64')))
-m.part1[:] = 1
-m.part2[:] = 2
-print(m.part1, m.part2)
-print(m)
-print()
 
-n = mymesh(m)
-n.part1[0] = 10
-print(n.part1, n.part2)
-print(n)
-print()
+exit()
 
-print(o.part1, o.part2)
-print(o)
-print()
-# print(n + m)

pySDC/tutorial/step_7/C_pySDC_with_PETSc.py

@@ -56,7 +56,7 @@ def main():
     problem_params['nu'] = 1.0  # diffusion coefficient
     problem_params['freq'] = 2  # frequency for the test value
     problem_params['cnvars'] = [(65, 65)]  # number of degrees of freedom for the coarsest level
-    problem_params['refine'] = [1, 0]  # number of refinements
+    problem_params['refine'] = [1]  # number of refinements
     problem_params['comm'] = space_comm  # pass space-communicator to problem class
     problem_params['sol_tol'] = 1E-12  # set tolerance to PETSc' linear solver