Updates

Dockerfile

@@ -22,7 +22,7 @@ RUN sed -i 's/Components: main/Components: main contrib non-free/' /etc/apt/sour
         libarpack2-dev \
         mpi-default-bin mpi-default-dev \
         python3 python3-dev python-is-python3 python3-pip \
-        python3-numpy python3-scipy python3-matplotlib python3-mpi4py cython3 python3-yaml python3-h5py python3-tk jupyter-notebook \
+        python3-numpy python3-scipy python3-matplotlib python3-mpi4py cython3 python3-yaml python3-h5py python3-tk jupyter-notebook python3-meshio python3-gmsh \
   --no-install-recommends \
   && rm -rf /var/lib/apt/lists/*
 
@@ -31,7 +31,8 @@ RUN sed -i 's/Components: main/Components: main contrib non-free/' /etc/apt/sour
 ENV OMPI_MCA_hwloc_base_binding_policy=hwthread \
     MPIEXEC_FLAGS=--allow-run-as-root \
     OMPI_ALLOW_RUN_AS_ROOT=1 \
-    OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1
+    OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1 \
+    PIP_BREAK_SYSTEM_PACKAGES=1
 
 COPY . /pynucleus
 

base/PyNucleus_base/LinearOperator_{SCALAR}.pxi

@@ -41,11 +41,15 @@ cdef class {SCALAR_label}LinearOperator:
                  BOOL_t no_overwrite=False,
                  BOOL_t trans=False):
         if not trans:
+            assert self.num_columns == x.shape[0], (self.num_columns, x.shape[0])
+            assert self.num_rows == y.shape[0], (self.num_rows, y.shape[0])
             if no_overwrite:
                 self.matvec_no_overwrite(x, y)
             else:
                 self.matvec(x, y)
         else:
+            assert self.num_rows == x.shape[0], (self.num_rows, x.shape[0])
+            assert self.num_columns == y.shape[0], (self.num_columns, y.shape[0])
             if no_overwrite:
                 self.matvecTrans_no_overwrite(x, y)
             else:
@@ -93,6 +97,9 @@ cdef class {SCALAR_label}LinearOperator:
     def __sub__(self, x):
         return self + (-1.*x)
 
+    def __matmul__(self, {SCALAR_label}LinearOperator x):
+        return {SCALAR_label}Product_Linear_Operator(self, x)
+
     def __mul__(self, x):
         cdef:
             np.ndarray[{SCALAR}_t, ndim=1] y
@@ -110,7 +117,7 @@ cdef class {SCALAR_label}LinearOperator:
                 return self.dotMV(x)
             elif isinstance(self, {SCALAR_label}TimeStepperLinearOperator) and isinstance(x, (float, int, {SCALAR})):
                 tsOp = self
-                return {SCALAR_label}TimeStepperLinearOperator(self, tsOp.M, tsOp.S, tsOp.facS*x, tsOp.facM*x)
+                return {SCALAR_label}TimeStepperLinearOperator(tsOp.M, tsOp.S, tsOp.facS*x, tsOp.facM*x)
             elif isinstance(self, {SCALAR_label}LinearOperator) and isinstance(x, (float, int, {SCALAR})):
                 return {SCALAR_label}Multiply_Linear_Operator(self, x)
             elif isinstance(x, {SCALAR_label}LinearOperator) and isinstance(self, (float, int, {SCALAR})):
@@ -128,8 +135,14 @@ cdef class {SCALAR_label}LinearOperator:
                 raise NotImplementedError('Cannot multiply {} with {}:\n{}'.format(self, x, e))
 
     def __rmul__(self, x):
+        cdef:
+            {SCALAR_label}TimeStepperLinearOperator tsOp
         if isinstance(x, (float, int, {SCALAR})):
-            return {SCALAR_label}Multiply_Linear_Operator(self, x)
+            if isinstance(self, {SCALAR_label}TimeStepperLinearOperator):
+                tsOp = self
+                return {SCALAR_label}TimeStepperLinearOperator(tsOp.M, tsOp.S, tsOp.facS*x, tsOp.facM*x)
+            else:
+                return {SCALAR_label}Multiply_Linear_Operator(self, x)
         else:
             raise NotImplementedError('Cannot multiply with {}'.format(x))
 
@@ -196,14 +209,14 @@ cdef class {SCALAR_label}LinearOperator:
     def toLinearOperator(self):
         def matvec(x):
             if x.ndim == 1:
-                return self*x
+                return self*x.astype({SCALAR})
             elif x.ndim == 2 and x.shape[1] == 1:
                 if x.flags.c_contiguous:
                     return self*x[:, 0]
                 else:
                     y = np.zeros((x.shape[0]), dtype=x.dtype)
                     y[:] = x[:, 0]
-                    return self*y
+                    return self*y.astype({SCALAR})
             else:
                 raise NotImplementedError()
 
@@ -380,9 +393,27 @@ cdef class {SCALAR_label}TimeStepperLinearOperator({SCALAR_label}LinearOperator)
 
     def __repr__(self):
         if np.real(self.facS) >= 0:
-            return '{}*{} + {}*{}'.format(self.facM, self.M, self.facS, self.S)
+            if self.facM != 1.0:
+                if self.facS != 1.0:
+                    return '{}*{} + {}*{}'.format(self.facM, self.M, self.facS, self.S)
+                else:
+                    return '{}*{} + {}'.format(self.facM, self.M, self.S)
+            else:
+                if self.facS != 1.0:
+                    return '{} + {}*{}'.format(self.M, self.facS, self.S)
+                else:
+                    return '{} + {}'.format(self.M, self.S)
         else:
-            return '{}*{} - {}*{}'.format(self.facM, self.M, -self.facS, self.S)
+            if self.facM != 1.0:
+                if self.facS != -1.0:
+                    return '{}*{} - {}*{}'.format(self.facM, self.M, -self.facS, self.S)
+                else:
+                    return '{}*{} - {}'.format(self.facM, self.M, self.S)
+            else:
+                if self.facS != -1.0:
+                    return '{} - {}*{}'.format(self.M, -self.facS, self.S)
+                else:
+                    return '{} - {}'.format(self.M, self.S)
 
     def to_csr_linear_operator(self):
         if isinstance(self.S, {SCALAR_label}Dense_LinearOperator):

base/PyNucleus_base/plot_utils.py

@@ -217,12 +217,12 @@ def formatScientificLatex(a, useEnotation=True):
         exp = int(np.floor(np.log10(a)))
         mantissa = a/10**exp
         if useEnotation:
-            return '{:.3}\mathrm{{e}}{{{}}}'.format(mantissa, exp)
+            return '{:.3}\\mathrm{{e}}{{{}}}'.format(mantissa, exp)
         else:
             return '{:.3} \\times 10^{{{}}}'.format(mantissa, exp)
     elif abs(a) == 0:
         if useEnotation:
-            return '0.00\mathrm{{e}}{0}'
+            return '0.00\\mathrm{{e}}{0}'
         else:
             return '0.00 \\times 10^{0}'
     else:

base/PyNucleus_base/utilsFem.py

@@ -1069,6 +1069,14 @@ def process(self, override={}):
         if self.comm:
             params = self.comm.bcast(params, root=0)
         self.params = params
+
+        if params['test']:
+            import psutil
+            p = psutil.Process()
+            try:
+                p.cpu_affinity(list(range(psutil.cpu_count())))
+            except AttributeError:
+                pass
         self._timer = TimerManager(self.logger, comm=self.comm, memoryProfiling=params['showMemory'])
 
         for fun in self.processHook:

drivers/brusselator.py

@@ -27,7 +27,7 @@
 from PyNucleus_base.linear_operators import TimeStepperLinearOperator
 from PyNucleus_fem.femCy import assembleNonlinearity
 from PyNucleus_multilevelSolver import hierarchyManager
-from PyNucleus_nl import paramsForFractionalHierarchy
+from PyNucleus_nl.helpers import paramsForFractionalHierarchy
 from PyNucleus_nl.nonlocalProblems import brusselatorProblem
 from PyNucleus_base.timestepping import EulerIMEX, ARS3, koto
 import h5py

drivers/brusselatorMovie.py

@@ -8,42 +8,43 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
-from path import Path
+from pathlib import Path
 from shutil import rmtree
 import h5py
 from subprocess import Popen
 from PyNucleus_base import driver
-from PyNucleus_fem import meshNd
+from PyNucleus_fem.mesh import meshNd
 from PyNucleus_fem.DoFMaps import DoFMap
 from PyNucleus_nl.nonlocalProblems import brusselatorProblem
 
 d = driver()
-brusselatorProblem(d)
+# brusselatorProblem(d)
+d.add('inputFile', '')
 d.add('zoomIn', False)
 d.add('shading', acceptedValues=['gouraud', 'flat'])
 d.process()
 
-filename = d.identifier+'.hdf5'
+filename = d.inputFile
 resultFile = h5py.File(str(filename), 'r')
 
 mesh = meshNd.HDF5read(resultFile['mesh'])
 dm = DoFMap.HDF5read(resultFile['dm'])
 dm.mesh = mesh
 
-folder = Path('brusselatorMovie')/Path(filename).basename()
+folder = Path('brusselatorMovie')/Path(filename).name
 try:
     rmtree(str(folder))
 except:
     pass
-folder.makedirs_p()
+folder.mkdir(parents=True, exist_ok=True)
 
 if d.zoomIn:
-    folderZoom = Path('brusselatorMovie')/(Path(filename).basename()+'-zoomIn')
+    folderZoom = Path('brusselatorMovie')/(Path(filename).name+'-zoomIn')
     try:
         rmtree(str(folderZoom))
     except:
         pass
-    folderZoom.makedirs_p()
+    folderZoom.mkdir(parents=True, exist_ok=True)
 
 l = sorted([int(i) for i in resultFile['U']])
 
@@ -74,12 +75,12 @@
 resultFile.close()
 
 Popen(['mencoder', 'mf://*.png', '-mf', 'fps=10', '-o',
-       '../{}.avi'.format(Path(filename).basename().stripext()), '-ovc', 'lavc',
+       '../{}.avi'.format(Path(filename).stem), '-ovc', 'lavc',
        '-lavcopts', 'vcodec=msmpeg4v2:vbitrate=800'],
       cwd=folder).wait()
 if d.zoomIn:
     Popen(['mencoder', 'mf://*.png', '-mf', 'fps=10', '-o',
-           '../{}-zoom.avi'.format(Path(filename).basename().stripext()), '-ovc', 'lavc',
+           '../{}-zoom.avi'.format(Path(filename).stem), '-ovc', 'lavc',
            '-lavcopts', 'vcodec=msmpeg4v2:vbitrate=800'],
           cwd=folderZoom).wait()
 

drivers/runParallelGMG.py

@@ -44,11 +44,13 @@
 solver.add('smoother', 'jacobi', acceptedValues=['gauss_seidel', 'chebyshev'])
 solver.add('doPCoarsen', False)
 solver.add('maxiter', 50)
+solver.add('tolerance', 0.)
 
 d.declareFigure('residuals', default=False)
+d.declareFigure('numericalSolution')
 d.declareFigure('spSolve')
 d.declareFigure('spSolveError')
-d.declareFigure('spSolveExactSolution')
+d.declareFigure('interpolatedAnalyticSolution')
 
 params = d.process()
 
@@ -95,7 +97,7 @@
     else:
         hierarchies, connectors = paramsForMG(p.noRef,
                                               range(d.comm.size),
-                                              params, p.dim, p.element)
+                                              params, p.manifold_dim, p.element)
         connectors['input'] = {'type': inputConnector,
                                'params': {'domain': d.domain}}
 
@@ -111,10 +113,13 @@
     overlaps = hM[FINE].multilevelAlgebraicOverlapManager
     h = hM[FINE].meshLevels[-1].mesh.global_h(overlaps.comm)
     hmin = hM[FINE].meshLevels[-1].mesh.global_hmin(overlaps.comm)
-    tol = {'P1': 0.5*h**2,
-           'P2': 0.001*h**3,
-           'P3': 0.001*h**4}[d.element]
-    tol = max(tol, 2e-9)
+    if d.tolerance <= 0.:
+        tol = {'P1': 0.5*h**2,
+               'P2': 0.001*h**3,
+               'P3': 0.001*h**4}[d.element]
+        tol = max(tol, 2e-9)
+    else:
+        tol = d.tolerance
 
 # assemble rhs on finest grid
 with d.timer('Assemble rhs on finest grid'):
@@ -211,6 +216,7 @@
         residuals = solver.residuals
         A.residual_py(x, rhs, r)
         resNorm = r.norm(False)
+        numIter = max(1, numIter)
         rate.add('Rate of convergence P'+label, (resNorm/r0)**(1/numIter), tested=False if label == 'BICGSTAB' else None)
         its.add('Number of iterations P'+label, numIter, aTol=2 if label == 'BICGSTAB' else None)
         res.add('Residual norm P'+label, resNorm)
@@ -333,11 +339,16 @@
         # cg(b_global, u_global)
         # print(cg.residuals, len(cg.residuals))
 
+if p.nontrivialNullspace:
+    const = DoFMap_fine.ones()
+    x -= (x.inner(const)/const.inner(const))*const
+
 if p.L2ex:
     if p.boundaryCond:
         d.logger.warning('L2 error is wrong for inhomogeneous BCs')
     with d.timer('Mass matrix'):
         M = DoFMap_fine.assembleMass(sss_format=d.symmetric)
+
     z = DoFMap_fine.assembleRHS(p.exactSolution)
     L2err = np.sqrt(np.absolute(x.inner(M*x, True, False) -
                                 2*z.inner(x, False, True) +
@@ -379,17 +390,30 @@
      global_dm) = accumulate2global(subdomain, interfaces, DoFMap_fine, x,
                                     comm=d.comm)
     if d.isMaster:
-        from scipy.sparse.linalg import spsolve
+
+        if d.startPlot('numericalSolution'):
+            global_solution.plot()
+        if p.exactSolution and d.startPlot('interpolatedAnalyticSolution'):
+            global_dm.interpolate(p.exactSolution).plot()
+
         from numpy.linalg import norm
         A = global_dm.assembleStiffness()
+        if p.nontrivialNullspace:
+            from PyNucleus_base.linear_operators import Dense_LinearOperator
+            A = A+Dense_LinearOperator(np.ones(A.shape))
         rhs = global_dm.assembleRHS(p.rhsFun)
         if p.boundaryCond:
             global_boundaryDoFMap = global_dm.getComplementDoFMap()
             global_boundary_data = global_boundaryDoFMap.interpolate(p.boundaryCond)
             global_A_boundary = global_dm.assembleStiffness(dm2=global_boundaryDoFMap)
             rhs -= global_A_boundary*global_boundary_data
         with d.timer('SpSolver'):
-            y = spsolve(A.to_csr(), rhs)
+            directSolver = solverFactory('lu', A=A, setup=True)
+            global_solution_spsolve = global_dm.zeros()
+            directSolver(rhs, global_solution_spsolve)
+            if p.nontrivialNullspace:
+                const = global_dm.ones()
+                global_solution_spsolve -= (global_solution_spsolve.inner(const)/const.inner(const))*const
         if p.boundaryCond:
             sol_augmented, dm_augmented = global_dm.augmentWithBoundaryData(global_solution, global_boundary_data)
             global_mass = dm_augmented.assembleMass()
@@ -405,14 +429,12 @@
                                         p.L2ex))
         errsSpSolve = d.addOutputGroup('errSpSolve')
         errsSpSolve.add('L2', L2err)
-        errsSpSolve.add('2-norm', norm(global_solution-y, 2))
-        errsSpSolve.add('max-norm', np.abs(global_solution-y).max())
+        errsSpSolve.add('2-norm', norm(global_solution-global_solution_spsolve, 2))
+        errsSpSolve.add('max-norm', np.abs(global_solution-global_solution_spsolve).max())
         d.logger.info('\n'+str(errsSpSolve))
         if d.startPlot('spSolve'):
             import matplotlib.pyplot as plt
-            global_solution.plot()
-        if p.exactSolution and d.startPlot('spSolveExactSolution'):
-            global_dm.interpolate(p.exactSolution).plot()
+            global_solution_spsolve.plot()
         if d.startPlot('spSolveError'):
-            (global_solution-y).plot()
+            (global_solution-global_solution_spsolve).plot()
 d.finish()

fem/PyNucleus_fem/DoFMaps.pyx

@@ -722,7 +722,8 @@ cdef class DoFMap:
                           BOOL_t reorder=False,
                           diffusivity=None,
                           INDEX_t[::1] cellIndices=None,
-                          DoFMap dm2=None):
+                          DoFMap dm2=None,
+                          simplexQuadratureRule qr=None):
         """This function assembles the stiffness matrix for a given
         diffusivity function:
 
@@ -760,7 +761,8 @@ cdef class DoFMap:
                                  reorder,
                                  diffusivity,
                                  cellIndices,
-                                 dm2=dm2)
+                                 dm2=dm2,
+                                 qr=qr)
 
     def assembleRHS(self,
                     fun,
@@ -983,7 +985,7 @@ cdef class DoFMap:
         for i in range(numCoords):
             for k in range(dim):
                 coords[i, k] = 0.
-            for j in range(dim+1):
+            for j in range(cell.shape[0]):
                 for k in range(dim):
                     coords[i, k] += self.nodes[i, j]*cell[j, k]