updates

base/PyNucleus_base/performanceLogger.pxd

@@ -24,8 +24,8 @@ cdef class Timer(FakeTimer):
         double startTime_unsynced
         double elapsed
         double elapsed_unsynced
-        double startMem
-        double endMem
+        readonly double startMem
+        readonly double endMem
         str key
         FakePLogger parent
         BOOL_t manualDataEntry

base/PyNucleus_base/performanceLogger.pyx

@@ -62,6 +62,8 @@ cdef class Timer(FakeTimer):
         else:
             self.memoryRegionsAreEnabled = memRegionsAreEnabled
         self.memoryProfiling = self.parent.memoryProfiling
+        self.startMem = 0.0
+        self.endMem = 0.0
         if self.sync:
             assert self.comm is not None
 
@@ -102,8 +104,15 @@ cdef class Timer(FakeTimer):
     def getIntervalUnsynced(self):
         return self.elapsed_unsynced
 
+    def getMemoryDiff(self):
+        if self.memoryProfiling:
+            return self.endMem-self.startMem
+        else:
+            return 0.0
+
     interval = property(fget=getInterval)
     interval_unsynced = property(fget=getIntervalUnsynced)
+    diffMemory = property(fget=getMemoryDiff)
 
 
 cdef class LoggingTimer(Timer):

base/PyNucleus_base/utilsFem.py

@@ -411,25 +411,22 @@ def updateFromDefaults(params, defaults):
             updateFromDefaults(params[key], defaults[key])
 
 
-KEY_VAL_FORMAT = '{:<54}{}'
-
-
-def getMPIinfo():
+def getMPIinfo(grp):
     from sys import modules
     if 'mpi4py.MPI' in modules:
         import mpi4py
         mpi4py.initialize = False
         from mpi4py import MPI
         if not MPI.Is_initialized():
-            return ''
+            return
         t = {MPI.THREAD_SINGLE: 'single',
              MPI.THREAD_FUNNELED: 'funneled',
              MPI.THREAD_SERIALIZED: 'serialized',
              MPI.THREAD_MULTIPLE: 'multiple'}
         hosts = MPI.COMM_WORLD.gather(MPI.Get_processor_name())
         if MPI.COMM_WORLD.rank == 0:
             hosts = ','.join(set(hosts))
-        s = ['{}'.format(MPI.Get_library_version()[:-1])]
+        grp.add('MPI library', '{}'.format(MPI.Get_library_version()[:-1]))
         for label, value in [('MPI standard supported:', MPI.Get_version()),
                              ('Vendor:', MPI.get_vendor()),
                              ('Level of thread support:', t[MPI.Query_thread()]),
@@ -438,13 +435,10 @@ def getMPIinfo():
                              ('Thread level requested:', mpi4py.rc.thread_level),
                              ('Hosts:', hosts),
                              ('Communicator size:', MPI.COMM_WORLD.size)]:
-            s.append(KEY_VAL_FORMAT.format(label, value))
-        return '\n'.join(s)
-    else:
-        return ''
+            grp.add(label, value)
 
 
-def getEnvVariables(envVars=[('OMP_NUM_THREADS', True)]):
+def getEnvVariables(grp, envVars=[('OMP_NUM_THREADS', True)]):
     from os import environ
     s = []
     for var, printNotSet in envVars:
@@ -454,23 +448,22 @@ def getEnvVariables(envVars=[('OMP_NUM_THREADS', True)]):
             varVal = 'not set'
         else:
             continue
-        s.append(KEY_VAL_FORMAT.format(var+':', varVal))
+        grp.add(var, varVal)
     return '\n'.join(s)
 
 
-def getSystemInfo(argv=None, envVars=[('OMP_NUM_THREADS', True)]):
+def getSystemInfo(grp, argv=None, envVars=[('OMP_NUM_THREADS', True)]):
     from sys import executable
-    s = '\n'
     if argv is not None:
-        s += KEY_VAL_FORMAT.format('Running:', executable + ' ' + ' '.join(argv)) + '\n'
+        grp.add('Running', executable + ' ' + ' '.join(argv))
     else:
-        s += KEY_VAL_FORMAT.format('Running:', executable) + '\n'
+        grp.add('Running', executable)
     import mpi4py
     mpi4py.initialize = False
     from mpi4py import MPI
     if MPI.Is_initialized():
-        s += getMPIinfo()+'\n'
-    s += getEnvVariables(envVars)+'\n'
+        getMPIinfo(grp)
+    getEnvVariables(grp, envVars)
     import pkg_resources
     from PyNucleus import subpackages
     versions = {}
@@ -481,17 +474,21 @@ def getSystemInfo(argv=None, envVars=[('OMP_NUM_THREADS', True)]):
         except KeyError:
             versions[version] = [pkg]
     for version in versions:
-        s += KEY_VAL_FORMAT.format(','.join(versions[version])+':', version)+'\n'
+        grp.add(','.join(versions[version]), version)
+
+    import importlib
+
     versions = {}
     for pkg in sorted(subpackages.keys()):
         version = pkg_resources.get_distribution('PyNucleus_'+pkg).version
+        module = importlib.import_module('PyNucleus_'+pkg+'.config')
+        sha = module.gitSHA
         try:
-            versions[version].append(pkg)
+            versions[(version, sha)].append(pkg)
         except KeyError:
-            versions[version] = [pkg]
+            versions[(version, sha)] = [pkg]
     for version in versions:
-        s += KEY_VAL_FORMAT.format('PyNucleus_'+(','.join(versions[version]))+':', version)+'\n'
-    return s
+        grp.add('PyNucleus_'+(','.join(versions[version])), version)
 
 
 class MPIFileHandler(logging.Handler):
@@ -504,6 +501,7 @@ def __init__(self, filename, comm, mode=MPI.MODE_WRONLY | MPI.MODE_CREATE):
         # keep the absolute path, otherwise derived classes which use this
         # may come a cropper when the current directory changes
         self.baseFilename = os.path.abspath(filename)
+        assert len(self.baseFilename) <= 245, 'The length of the log file path \"{}\" is too long and will probably crash MPI. Try running with \"--disableFileLog\"'.format(self.baseFilename)
         if Path(self.baseFilename).exists() and comm.rank == 0:
             from os import remove
             remove(self.baseFilename)
@@ -1095,11 +1093,12 @@ def process(self, override={}):
         if params['displayConfig']:
             from pprint import pformat
             self.logger.info('\n'+pformat(params))
+
+        from sys import argv
+        sysInfo = self.addOutputGroup('sysInfo')
+        getSystemInfo(argv=argv, grp=sysInfo)
         if not params['disableHeader']:
-            from sys import argv
-            sysInfo = getSystemInfo(argv)
-            if self.isMaster:
-                self.logger.info(sysInfo)
+            self.logger.info('\n'+str(sysInfo))
         if params['logDependencies']:
             dependencyLogger.setLevel(logging.DEBUG)
         else:

docs/example1.rst

@@ -3,7 +3,7 @@ Example 1 - A simple PDE problem
 ================================
 
 In this first example, we will construct a finite element discretization of a classical PDE problem and solve it.
-The full code of this example can be found in ``drivers/example1.py``.
+The full code of this example can be found in ``examples/example1.py``.
 
 Factories
 ---------

docs/example2.rst

@@ -3,7 +3,7 @@ Example 2 - Nonlocal problems
 =============================
 
 I this second example, we will assemble and solve several nonlocal equations.
-The full code of this example can be found in ``drivers/example2.py``.
+The full code of this example can be found in ``examples/example2.py``.
 
 PyNucleus can assemble operators of the form
 

drivers/example2.py → examples/example2.py

@@ -94,8 +94,8 @@
 dmConst = dofmapFactory('P1', meshConst, nIConst['domain'])
 dmConstInteraction = dmConst.getComplementDoFMap()
 
-A_const = dmConst.assembleNonlocal(kernelConst, matrixFormat='sparsified')
-B_const = dmConst.assembleNonlocal(kernelConst, dm2=dmConstInteraction, matrixFormat='sparsified')
+A_const = dmConst.assembleNonlocal(kernelConst, matrixFormat='sparse')
+B_const = dmConst.assembleNonlocal(kernelConst, dm2=dmConstInteraction, matrixFormat='sparse')
 
 g = functionFactory('Lambda', lambda x: -(x[0]**2 + x[1]**2)/4)
 g_interp = dmConstInteraction.interpolate(g)
@@ -109,7 +109,7 @@
 
 solver(b, u)
 
-u_global = dmConst.augmentWithBoundaryData(u, g_interp)
+u_global = u.augmentWithBoundaryData(g_interp)
 
 plt.figure().gca().set_title('Numerical solution, constant kernel')
 u_global.plot()
@@ -121,3 +121,4 @@
 
 ######################################################################
 plt.show()
+

nl/PyNucleus_nl/clusterMethodCy.pxd

@@ -13,6 +13,7 @@ from PyNucleus_fem.quadrature cimport (simplexDuffyTransformation,
 from PyNucleus_fem.functions cimport function
 from PyNucleus_base.myTypes cimport INDEX_t, REAL_t, BOOL_t
 from PyNucleus_base.linear_operators cimport (LinearOperator,
+                                              VectorLinearOperator,
                                               Dense_LinearOperator,
                                               CSR_LinearOperator,
                                               SSS_LinearOperator)
@@ -76,6 +77,7 @@ cdef class tree_node:
         public REAL_t[:, ::1] transferOperator
         public REAL_t[:, :, ::1] value
         public REAL_t[::1] coefficientsUp, coefficientsDown
+        public REAL_t[:, ::1] coefficientsDownVec
         public BOOL_t mixed_node
         public BOOL_t canBeAssembled
         public INDEX_t levelNo
@@ -89,15 +91,17 @@ cdef class tree_node:
     cdef BOOL_t get_is_leaf(self)
     cdef INDEX_t _getLevels(self)
     cdef INDEX_t _getParentLevels(self)
-    cdef void prepareTransferOperators(self, INDEX_t m, transferMatrixBuilder tMB=*)
+    cdef void prepareTransferOperators(self, INDEX_t m, INDEX_t valueSize, transferMatrixBuilder tMB=*)
     cdef void upwardPass(self, REAL_t[::1] x, INDEX_t componentNo=*, BOOL_t skip_leaves=*, BOOL_t local=*)
-    cdef void resetCoefficientsDown(self)
+    cdef void resetCoefficientsDown(self, BOOL_t vecValued=*)
     cdef void resetCoefficientsUp(self)
     cdef void downwardPass(self, REAL_t[::1] y, INDEX_t componentNo=*, BOOL_t local=*)
+    cdef void downwardPassVec(self, REAL_t[:, ::1] y, INDEX_t componentNo=*, BOOL_t local=*)
     cpdef INDEX_t findCell(self, meshBase mesh, REAL_t[::1] vertex, REAL_t[:, ::1] simplex, REAL_t[::1] bary)
     cpdef set findCells(self, meshBase mesh, REAL_t[::1] vertex, REAL_t r, REAL_t[:, ::1] simplex)
     cdef BOOL_t trim(self, bitArray keep)
     cdef void upwardPassMatrix(self, dict coefficientsUp)
+    cpdef void relabelDoFs(self, INDEX_t[::1] translate)
 
 
 cdef class productIterator:
@@ -113,7 +117,9 @@ cdef class farFieldClusterPair:
     cdef:
         public tree_node n1, n2
         public REAL_t[:, ::1] kernelInterpolant
+        public REAL_t[:, :, ::1] kernelInterpolantVec
     cpdef void apply(self, REAL_t[::1] x, REAL_t[::1] y)
+    cpdef void applyVec(farFieldClusterPair self, REAL_t[::1] x, REAL_t[:, ::1] y)
 
 
 cdef class H2Matrix(LinearOperator):
@@ -135,6 +141,25 @@ cdef class H2Matrix(LinearOperator):
                                tree_node left) except -1
 
 
+cdef class VectorH2Matrix(VectorLinearOperator):
+    cdef:
+        public VectorLinearOperator Anear
+        public dict Pfar
+        public tree_node tree
+        public FakePLogger PLogger
+        public LinearOperator basis
+        public BOOL_t skip_leaves_upward
+        public REAL_t[::1] leafCoefficientsUp
+    # cdef INDEX_t matvec(self,
+    #                     REAL_t[::1] x,
+    #                     REAL_t[:, ::1] y) except -1
+    # cdef INDEX_t matvec_submat(self,
+    #                            REAL_t[::1] x,
+    #                            REAL_t[:, ::1] y,
+    #                            list right_list,
+    #                            tree_node left) except -1
+
+
 cdef class DistributedH2Matrix_globalData(LinearOperator):
     cdef:
         public LinearOperator localMat