clean ups

Dockerfile

@@ -50,3 +50,6 @@ ENV OMPI_MCA_hwloc_base_binding_policy=hwthread \
     OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1
 
 RUN python -m ipykernel install --name=PyNucleus
+
+COPY README.container.rst /README.container.rst
+RUN echo '[ ! -z "$TERM" -a -r /README.container.rst ] && cat /README.container.rst' >> /etc/bash.bashrc

PyNucleus/__init__.py

@@ -5,19 +5,51 @@
 # If you want to use this code, please refer to the README.rst and LICENSE files. #
 ###################################################################################
 
-import importlib
-import pkgutil
 import sys
 
 subpackages = {}
 __all__ = []
-for finder, name, ispkg in pkgutil.iter_modules():
-    if ispkg and name.find('PyNucleus_') == 0:
-        importName = name[len('PyNucleus_'):]
-        module = importlib.import_module(name, 'PyNucleus')
-        sys.modules['PyNucleus.'+importName] = module
-        subpackages[importName] = module
-        names = [name for name in module.__dict__ if not name.startswith('_')]
-        locals().update({name: getattr(module, name) for name in names})
-        if hasattr(module, '__all__'):
-            __all__ += module.__all__
+
+from PyNucleus_packageTools import *
+import PyNucleus_packageTools as packageTools
+
+sys.modules['PyNucleus.packageTools'] = packageTools
+subpackages['packageTools'] = packageTools
+
+from PyNucleus_base import *
+import PyNucleus_base as base
+
+sys.modules['PyNucleus.base'] = base
+subpackages['base'] = base
+__all__ += base.__all__
+
+from PyNucleus_metisCy import *
+import PyNucleus_metisCy as metisCy
+
+sys.modules['PyNucleus.metisCy'] = metisCy
+subpackages['metisCy'] = metisCy
+__all__ += metisCy.__all__
+
+from PyNucleus_fem import *
+import PyNucleus_fem as fem
+
+sys.modules['PyNucleus.fem'] = fem
+subpackages['fem'] = fem
+__all__ += fem.__all__
+
+from PyNucleus_multilevelSolver import *
+import PyNucleus_multilevelSolver as multilevelSolver
+
+sys.modules['PyNucleus.multilevelSolver'] = multilevelSolver
+subpackages['multilevelSolver'] = multilevelSolver
+__all__ += multilevelSolver.__all__
+
+try:
+    from PyNucleus_nl import *
+    import PyNucleus_nl as nl
+
+    sys.modules['PyNucleus.nl'] = nl
+    subpackages['nl'] = nl
+    __all__ += nl.__all__
+except ImportError:
+    pass

README.container.rst

@@ -0,0 +1,6 @@
+
+This is a container image for PyNucleus.
+
+The drivers and examples for PyNucleus can be found in /pynucleus/drivers and /pynucleus/examples
+
+The directory from which the container was launched on the host system is mapped to /user

base/PyNucleus_base/utilsFem.py

@@ -1822,7 +1822,9 @@ def setDriverArgs(self):
 
     @property
     def timer(self):
-        return self.driver.getTimer()
+        if self._timer is None:
+            self._timer = self.driver.getTimer()
+        return self._timer
 
     def processCmdline(self, params):
         pass

compose.yaml

@@ -3,7 +3,7 @@ version: 3
 services:
 
   # Launch with:
-  # docker compose run --profile interactive pynucleus
+  # docker compose run pynucleus
   pynucleus:
     image: ghcr.io/sandialabs/pynucleus:latest
     build: .
@@ -22,23 +22,24 @@ services:
       - /tmp/.X11-unix:/tmp/.X11-unix
       - $XAUTHORITY:/root/.Xauthority
     network_mode: host
+    hostname: pynucleus-container
 
   # Launch with:
   # docker compose up
   # Then open localhost:8889 in your browser
   pynucleus-jupyter:
     image: ghcr.io/sandialabs/pynucleus:latest
     build: .
-    command: jupyter notebook --port=8888 --no-browser --ip=0.0.0.0 --allow-root --NotebookApp.token='' --NotebookApp.password='' --notebook-dir=/pynucleus/user --KernelSpecManager.ensure_native_kernel=False --KernelSpecManager.allowed_kernelspecs=pynucleus
+    command: jupyter notebook --port=8888 --no-browser --ip=0.0.0.0 --allow-root --NotebookApp.token='' --NotebookApp.password='' --notebook-dir=/notebooks --KernelSpecManager.ensure_native_kernel=False --KernelSpecManager.allowed_kernelspecs=pynucleus
     environment:
       # expose host proxies
       - http_proxy=${http_proxy}
       - https_proxy=${https_proxy}
       - HTTP_PROXY=${HTTP_PROXY}
       - HTTPS_PROXY=${HTTPS_PROXY}
     volumes:
-      # The current directory on host gets mapped to /pynucleus/user in the container
-      - $PWD/notebooks:/pynucleus/user
+      # The notebook subdirectory on host gets mapped to /notebooks in the container
+      - $PWD/notebooks:/notebooks
     ports:
       # Expose a Jupyter notebook server from the container
       - 8889:8888

drivers/runNonlocalInterface.py

@@ -92,9 +92,9 @@
 with d.timer('assemble matrices'):
     buildDense = True
 
-    A1 = getFracLapl(domain1Mesh, dm1, nIP.kernel1, boundaryCondition=HOMOGENEOUS_NEUMANN, dense=buildDense,
+    A1 = getFracLapl(dm1, nIP.kernel1, boundaryCondition=HOMOGENEOUS_NEUMANN, dense=buildDense,
                      forceRebuild=d.forceRebuild, genKernel=d.genKernel, trySparsification=True, doSave=True)
-    A2 = getFracLapl(domain2Mesh, dm2, nIP.kernel2, boundaryCondition=HOMOGENEOUS_NEUMANN, dense=buildDense,
+    A2 = getFracLapl(dm2, nIP.kernel2, boundaryCondition=HOMOGENEOUS_NEUMANN, dense=buildDense,
                      forceRebuild=d.forceRebuild, genKernel=d.genKernel, trySparsification=True, doSave=True)
 
 # domain-domain interaction

fem/PyNucleus_fem/DoFMaps.pyx

@@ -805,11 +805,19 @@ cdef class DoFMap:
         return assembleRHSgrad(fun, self, coeff, qr)
 
     def assembleNonlocal(self, kernel, str matrixFormat='DENSE', DoFMap dm2=None, BOOL_t returnNearField=False, **kwargs):
-        """Assemble a nonlocal operator of the form
+        """Assemble a nonlocal operator.
+
+        For finite horizon kernels
+
+        .. math::
+
+           \\iint_{D \\times D} (u(x)-u(y)) (v(x) \\gamma(x, y) - v(y) \\gamma(y, x))  dy dx
+
+        and for infinite horizon kernels
 
         .. math::
 
-           \\int_D (u(x)-u(y)) \\gamma(x, y) dy
+           \\iint_{D \\times D} (u(x)-u(y)) (v(x) \\gamma(x, y) - v(y) \\gamma(y, x))  dy dx + 2\\int_{D} u(x) v(x) \\int_{D^c} \\gamma(x, y) dy dx
 
         :param kernel: The kernel function :math:`\\gamma`
 
@@ -861,11 +869,11 @@ cdef class DoFMap:
                 if isinstance(kernel, ComplexKernel):
                     from PyNucleus_nl.nonlocalAssembly import ComplexnonlocalBuilder
 
-                    builder = ComplexnonlocalBuilder(self.mesh, self, kernel, dm2=dm2, **kwargs)
+                    builder = ComplexnonlocalBuilder(self, kernel, dm2=dm2, **kwargs)
                 else:
                     from PyNucleus_nl.nonlocalAssembly import nonlocalBuilder
 
-                    builder = nonlocalBuilder(self.mesh, self, kernel, dm2=dm2, **kwargs)
+                    builder = nonlocalBuilder(self, kernel, dm2=dm2, **kwargs)
                 if matrixFormat.upper() == 'DENSE':
                     return builder.getDense()
                 elif matrixFormat.upper() == 'DIAGONAL':

fem/PyNucleus_fem/__init__.py

@@ -14,4 +14,10 @@
 from . factories import functionFactory, dofmapFactory, meshFactory
 from . pdeProblems import diffusionProblem, helmholtzProblem
 __all__ = ['functionFactory', 'dofmapFactory', 'meshFactory',
-           'diffusionProblem', 'helmholtzProblem']
+           'diffusionProblem', 'helmholtzProblem',
+           'PHYSICAL', 'INTERIOR_NONOVERLAPPING', 'INTERIOR', 'NO_BOUNDARY',
+           'DIRICHLET', 'HOMOGENEOUS_DIRICHLET',
+           'NEUMANN', 'HOMOGENEOUS_NEUMANN',
+           'NORM', 'boundaryConditions',
+           'P0_DoFMap', 'P1_DoFMap', 'P2_DoFMap', 'P3_DoFMap',
+           'str2DoFMap', 'str2DoFMapOrder', 'getAvailableDoFMaps']

metisCy/PyNucleus_metisCy/__init__.py

@@ -82,3 +82,4 @@
 from . metisCy import (OBJTYPE_CUT,
                        OBJTYPE_NODE,
                        OBJTYPE_VOL)
+__all__ = []

multilevelSolver/PyNucleus_multilevelSolver/__init__.py

@@ -17,3 +17,4 @@
 
 solverFactory.register('mg', multigrid, isMultilevelSolver=True)
 solverFactory.register('complex_mg', Complexmultigrid, isMultilevelSolver=True)
+__all__ = []

nl/PyNucleus_nl/__init__.py

@@ -31,6 +31,6 @@
                                 nonlocalMeshFactory)
 from . discretizedProblems import (discretizedNonlocalProblem,
                                    discretizedTransientProblem)
-__all__ = [twoPointFunctionFactory, fractionalOrderFactory, interactionFactory, kernelFactory, nonlocalMeshFactory,
-           fractionalLaplacianProblem, nonlocalPoissonProblem, transientFractionalProblem,
-           discretizedNonlocalProblem, discretizedTransientProblem]
+__all__ = ['twoPointFunctionFactory', 'fractionalOrderFactory', 'interactionFactory', 'kernelFactory', 'nonlocalMeshFactory',
+           'fractionalLaplacianProblem', 'nonlocalPoissonProblem', 'transientFractionalProblem',
+           'discretizedNonlocalProblem', 'discretizedTransientProblem']

nl/PyNucleus_nl/discretizedProblems.py

@@ -354,6 +354,7 @@ def __init__(self, driver, continuumProblem):
         super().__init__(driver)
         self.continuumProblem = continuumProblem
         self.addRemote(self.continuumProblem)
+        driver.addToProcessHook(self.setTimerManager)
 
     def setDriverArgs(self):
         p = self.driver.addGroup('solver')
@@ -367,6 +368,9 @@ def setDriverArgs(self):
         self.setDriverFlag('matrixFormat', acceptedValues=['H2', 'sparse', 'sparsified', 'dense'], help='matrix format', group=p)
         self.setDriverFlag('debugAssemblyTimes', False, group=p)
 
+    def setTimerManager(self, params):
+        self._timer = self.driver.getTimer().getSubManager(logging.getLogger(__name__))
+
     @generates(['meshHierarchy', 'finalMesh',
                 'dm', 'dmBC', 'dmInterior',
                 'R_interior', 'P_interior',
@@ -518,7 +522,7 @@ def buildBCoperator(self, dmInterior, dmBC,
             assemblyParams['dense'] = matrixFormat == 'dense'
             assemblyParams['matrixFormat'] = matrixFormat
             assemblyParams['tag'] = tag
-            self.A_BC = getFracLapl(dmInterior.mesh, dmInterior, dm2=dmBC, **assemblyParams)
+            self.A_BC = getFracLapl(dmInterior, dm2=dmBC, **assemblyParams)
         else:
             self.A_BC = None
 
@@ -638,7 +642,7 @@ def solve(self, b, dm, dmInterior, dmBC, P_interior, P_bc, solver, boundaryCondi
     @generates('adjointModelSolution')
     def adjointSolve(self, b, dm, dmInterior, P_interior, adjointSolver, tol, maxiter):
         uInterior = dmInterior.zeros()
-        with self.timer('solve {}'.format(self.__class__.__name__)):
+        with self.timer('solve adjoint {}'.format(self.__class__.__name__)):
             its = adjointSolver(b, uInterior)
 
         resError = (b-adjointSolver.A*uInterior).norm(False)

nl/PyNucleus_nl/helpers.py

@@ -114,8 +114,9 @@ def processBC(tag, boundaryCondition, kernel):
     return tag, zeroExterior
 
 
-def getFracLapl(mesh, DoFMap, kernel=None, rangedOpParams={}, **kwargs):
+def getFracLapl(DoFMap, kernel=None, rangedOpParams={}, **kwargs):
 
+    mesh = DoFMap.mesh
     if kernel is None and len(rangedOpParams) == 0:
         return DoFMap.assembleStiffness(dm2=kwargs.get('dm2', None))
     assert kernel is not None or 's' in rangedOpParams, (kernel, rangedOpParams)
@@ -188,7 +189,7 @@ def getFracLapl(mesh, DoFMap, kernel=None, rangedOpParams={}, **kwargs):
                 intervalOps = []
                 for s in n:
                     kernel = getFractionalKernel(mesh.dim, constFractionalOrder(s), horizon, scaling=scaling, normalized=normalized)
-                    intervalOps.append(delayedFractionalLaplacianOp(mesh, DoFMap, kernel, **kwargs))
+                    intervalOps.append(delayedFractionalLaplacianOp(DoFMap, kernel, **kwargs))
                 ops.append(intervalOps)
             A = multiIntervalInterpolationOperator(intervals, nodes, ops)
             return A
@@ -263,10 +264,10 @@ def getFracLapl(mesh, DoFMap, kernel=None, rangedOpParams={}, **kwargs):
         if dm2 is not None and matrixFormat.upper() in ('H2', 'SPARSE', 'SPARSIFIED') and DoFMap.num_boundary_dofs > 0:
             # currently not implemented
             dm, R_interior, R_bc = DoFMap.getFullDoFMap(dm2)
-            A = getFracLapl(mesh, dm, kernel, rangedOpParams={}, **kwargs)
+            A = getFracLapl(dm, kernel, rangedOpParams={}, **kwargs)
             A = R_interior*A*R_bc.transpose()
             return A
-        builder = nonlocalBuilder(mesh, DoFMap, kernel, params, zeroExterior=zeroExterior, comm=comm, logging=logging, PLogger=PLogger, dm2=dm2)
+        builder = nonlocalBuilder(DoFMap, kernel, params, zeroExterior=zeroExterior, comm=comm, logging=logging, PLogger=PLogger, dm2=dm2)
         if diagonal:
             with timer('Assemble diagonal matrix {}, zeroExterior={}'.format(kernel, zeroExterior)):
                 A = builder.getDiagonal()
@@ -343,7 +344,7 @@ def build(self, buildType):
                     self.params['assemblyComm'] = self.comm
                     self.params['assembleOnRoot'] = False
                     self.params['forceUnsymmetric'] = True
-                self.S = getFracLapl(mesh, DoFMap, **self.params)
+                self.S = getFracLapl(DoFMap, **self.params)
                 self.A = self.S
                 # if not s.symmetric:
                 #     from PyNucleus_base.linear_operators import Dense_LinearOperator
@@ -650,10 +651,9 @@ def construct(self):
 
 
 class delayedFractionalLaplacianOp(delayedConstructionOperator):
-    def __init__(self, mesh, dm, kernel, **kwargs):
+    def __init__(self, dm, kernel, **kwargs):
         super().__init__(dm.num_dofs,
                          dm.num_dofs)
-        self.mesh = mesh
         self.dm = dm
         self.kernel = kernel
         self.kwargs = kwargs
@@ -662,7 +662,7 @@ def construct(self):
         from copy import copy
         d = copy(self.kwargs)
         d.update(self.params)
-        A = getFracLapl(self.mesh, self.dm, self.kernel, **d)
+        A = getFracLapl(self.dm, self.kernel, **d)
         return A
 
 

nl/PyNucleus_nl/nonlocalAssembly.pyx

@@ -196,7 +196,7 @@ cdef class horizonCorrected(TimeStepperLinearOperator):
         if Ainf is None:
             scaling = constantTwoPoint(0.5)
             infiniteKernel = kernel.getModifiedKernel(horizon=constant(np.inf), scaling=scaling)
-            infBuilder = nonlocalBuilder(self.mesh, self.dm, infiniteKernel, zeroExterior=True, comm=self.comm, logging=self.logging)
+            infBuilder = nonlocalBuilder(self.dm, infiniteKernel, zeroExterior=True, comm=self.comm, logging=self.logging)
             self.Ainf = infBuilder.getH2()
         else:
             self.Ainf = Ainf
@@ -233,7 +233,7 @@ cdef class horizonCorrected(TimeStepperLinearOperator):
         self.kernel = kernel
 
         complementKernel = kernel.getComplementKernel()
-        builder = nonlocalBuilder(self.mesh, self.dm, complementKernel, zeroExterior=True, comm=self.comm, logging=self.logging)
+        builder = nonlocalBuilder(self.dm, complementKernel, zeroExterior=True, comm=self.comm, logging=self.logging)
         correction = builder.getH2()
 
         self.S = self.Ainf
@@ -363,7 +363,7 @@ def assembleNonlocalOperator(meshBase mesh,
                              MPI.Comm comm=None,
                              **kwargs):
     kernel = getFractionalKernel(mesh.dim, s, horizon)
-    builder = nonlocalBuilder(mesh, dm, kernel, params, zeroExterior, comm, **kwargs)
+    builder = nonlocalBuilder(dm, kernel, params, zeroExterior, comm, **kwargs)
     return builder.getDense()
 
 
@@ -487,7 +487,7 @@ def assembleNearField(list Pnear,
                       comm=None,
                       **kwargs):
     kernel = getFractionalKernel(mesh.dim, s, horizon)
-    builder = nonlocalBuilder(mesh, dm, kernel, params, zeroExterior, comm, logging=True, **kwargs)
+    builder = nonlocalBuilder(dm, kernel, params, zeroExterior, comm, logging=True, **kwargs)
     A = builder.assembleClusters(Pnear)
     return A
 

nl/PyNucleus_nl/nonlocalAssembly_{SCALAR}.pxi

@@ -858,7 +858,6 @@ cdef class {SCALAR_label}IndexManagerMixed({SCALAR_label}IndexManager):
 
 cdef class {SCALAR_label}nonlocalBuilder:
     def __init__(self,
-                 meshBase mesh,
                  DoFMap dm,
                  {SCALAR_label}Kernel kernel,
                  dict params={},
@@ -875,7 +874,6 @@ cdef class {SCALAR_label}nonlocalBuilder:
 
         self.dm = dm
         self.mesh = self.dm.mesh
-        assert self.dm.mesh == mesh
         if dm2 is not None:
             self.dm2 = dm2
             assert type(self.dm) == type(self.dm2)
@@ -889,7 +887,6 @@ cdef class {SCALAR_label}nonlocalBuilder:
         self.params = params
 
         assert isinstance(self.kernel.horizon, constant), "Need horizon to be constant."
-        assert kernel.dim == mesh.dim, "Kernel dimension must match mesh dimension"
         assert kernel.dim == dm.mesh.dim, "Kernel dimension must match dm.mesh dimension"
 
         # volume integral

tests/test_fracLapl.py

@@ -94,8 +94,8 @@ def scaling(dim, s, horizon, refinements):
     kernel2 = getFractionalKernel(mesh.dim, s, horizon, scaling=scaling)
     print(kernel1, kernel2)
     zeroExterior = not np.isfinite(horizon.value)
-    builder1 = nonlocalBuilder(mesh, dm, kernel1, zeroExterior=zeroExterior)
-    builder2 = nonlocalBuilder(mesh, dm, kernel2, zeroExterior=zeroExterior)
+    builder1 = nonlocalBuilder(dm, kernel1, zeroExterior=zeroExterior)
+    builder2 = nonlocalBuilder(dm, kernel2, zeroExterior=zeroExterior)
     A = builder1.getDense().toarray()
     B = builder2.getDense().toarray()
     assert np.allclose(A, B)
@@ -160,7 +160,7 @@ def h2(dim, s, refinements, element, errBnd, genKernel=False):
     params['eta'] = eta
     params['maxLevels'] = maxLevels
     kernel = getFractionalKernel(mesh.dim, s, constant(np.inf))
-    builder = nonlocalBuilder(mesh, DoFMap_fine, kernel, params=params, zeroExterior=True)
+    builder = nonlocalBuilder(DoFMap_fine, kernel, params=params, zeroExterior=True)
 
     A_d = np.array(builder.getDense().data)
     A_h2 = builder.getH2()