updates

Author: cgcgcg

README.rst

@@ -133,12 +133,8 @@ After cloning the source code, PyNucleus is installed via
 
    make
 
-The compilation of PyNucleus can be configured by modifying the file `config.yaml` in the root folder.
+The compilation of PyNucleus can be configured by modifying the file `config.yaml <https://github.com/sandialabs/PyNucleus/blob/master/config.yaml>`_ in the root folder.
 This allows for example to set paths for libraries that are installed in non-standard directories.
-The defaults are as follows:
-
-.. literalinclude:: ../config.yaml
-   :language: yaml
 
 If you want to easily modify the source code without re-installing the package every time, and editable install is available as
 

base/PyNucleus_base/linear_operators.pyx

@@ -882,6 +882,16 @@ cdef class blockOperator(LinearOperator):
                   self.blockInptrRight[j]:self.blockInptrRight[j+1]] = lo.toarray()
         return B
 
+    def isSparse(self):
+        cdef:
+            BOOL_t sparse = True
+            INDEX_t i, j
+        for i in range(self.blockShape[0]):
+            for j in range(self.blockShape[1]):
+                lo = self.subblocks[i][j]
+                sparse &= lo.isSparse()
+        return sparse
+
 
 cdef class blockDiagonalOperator(blockOperator):
     def __init__(self, list diagonalBlocks):

base/PyNucleus_base/utilsFem.py

@@ -1028,7 +1028,7 @@ def process(self, override={}):
         doTerminate = False
         if self.isMaster:
             if 'plots' in self.argGroups:
-                io = self.addGroup('input/output')
+                io = self.addGroup('plots')
                 io.add('plotFolder', '', help='folder for saving plots')
                 io.add('plotFormat', acceptedValues=['pdf', 'png', 'jpeg', 'eps', 'ps', 'svg'], help='File format for saving plots')
             try:
@@ -1450,6 +1450,12 @@ def interpret(self, s):
         m = self.regexp.match(s)
         return [p(v) for v, p in zip(m.groups(), self.params)]
 
+    def __repr__(self):
+        params = []
+        for p in self.params:
+            params.append(p.__name__)
+        return "{}({})".format(self.name, ','.join(params))
+
 
 class propertyBuilder:
     def __init__(self, baseObj, fun):
@@ -1782,7 +1788,7 @@ def interpreter(v):
             for p in parametrizedArgs:
                 if self.parametrizedArg(p).match(v):
                     return v
-            raise ArgumentTypeError()
+            raise ArgumentTypeError("\"{}\" is not in list of accepted values {} or cannot be interpreted as parametrized arg {}.".format(v, acceptedValues, [repr(self.parametrizedArg(p)) for p in parametrizedArgs]))
 
         return interpreter
 

docs/drivers.rst

@@ -0,0 +1,55 @@
+
+Drivers
+=======
+
+The code repository contains several sample problems that can be found in the `drivers <https://github.com/sandialabs/PyNucleus/blob/master/drivers>`_ subfolder.
+The drivers take several command line parameters that can be used to change the problem and the type of outputs.
+A listing of parameters can be displayed by passing the `--help` flag:
+
+.. code-block:: shell
+
+   drivers/runFractional.py --help
+
+.. program-output:: python3 ../drivers/runFractional.py --help
+
+Some of the drivers can be run in parallel using MPI, e.g.
+
+.. code-block:: shell
+
+   mpiexec -n 4 drivers/runFractional.py --domain=disc
+
+
+runFractional.py
+----------------
+
+Assembles and solves fractional Poisson problems with infinite horizon.
+
+runNonlocal.py
+----------------
+
+Assembles and solves nonlocal Poisson problems with finite horizon.
+
+runFractionalHeat.py
+--------------------
+
+Solves a fractional heat equation with infinite horizon.
+
+runNonlocalInterface.py
+-----------------------
+
+A two domain interface problem with jumps in solution and flux for finite horizon kernels.
+
+brusselator.py
+--------------
+
+Solves a fractional-order Brusselator system.
+
+runParallelGMG.py
+-----------------
+
+Assembles and solves a classical local Poisson problem using geometric multigrid.
+
+runHelmholtz.py
+-----------------
+
+Assembles and solves a classical local Helmholtz problem.

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 ``examples/example1.py``.
+The full code of this example can be found in `examples/example1.py <https://github.com/sandialabs/PyNucleus/blob/master/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 ``examples/example2.py``.
+The full code of this example can be found in `examples/example2.py <https://github.com/sandialabs/PyNucleus/blob/master/examples/example2.py>`_.
 
 PyNucleus can assemble operators of the form
 

docs/example3.py

@@ -0,0 +1,109 @@
+#!/usr/bin/env python3
+###################################################################################
+# Copyright 2021 National Technology & Engineering Solutions of Sandia,           #
+# LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the           #
+# U.S. Government retains certain rights in this software.                        #
+# If you want to use this code, please refer to the README.rst and LICENSE files. #
+###################################################################################
+
+from PyNucleus.packageTools.sphinxTools import codeRegionManager
+
+mgr = codeRegionManager()
+
+with mgr.add('preamble'):
+    ######################################################################
+    # preamble
+    import logging
+    from PyNucleus_base.utilsFem import TimerManager
+    from PyNucleus import meshFactory, dofmapFactory, kernelFactory, functionFactory, solverFactory
+    from PyNucleus_nl.operatorInterpolation import admissibleSet
+    import h5py
+    from PyNucleus_base.linear_operators import LinearOperator
+    from PyNucleus_fem.DoFMaps import DoFMap
+
+    fmt = '{message}'
+    logging.basicConfig(level=logging.INFO,
+                        format=fmt,
+                        style='{',
+                        datefmt="%Y-%m-%d %H:%M:%S")
+    logger = logging.getLogger('__main__')
+    timer = TimerManager(logger)
+
+with mgr.add('setup'):
+    # Set up a mesh and a dofmap on it.
+    mesh = meshFactory('interval', hTarget=2e-3, a=-1, b=1)
+    dm = dofmapFactory('P1', mesh)
+
+    # Construct a RHS vector and a vector for the solution.
+    f = functionFactory('constant', 1.)
+    b = dm.assembleRHS(f)
+    u = dm.zeros()
+
+    # construct a fractional kernel with fractional order in S = [0.05, 0.95]
+    s = admissibleSet([0.05, 0.95])
+    kernel = kernelFactory('fractional', s=s, dim=mesh.dim)
+
+with mgr.add('operator'):
+    ######################################################################
+    # The operator is set up to be constructed on-demand.
+    # We partition the interval S into several sub-interval and construct a Chebyshev interpolant on each sub-interval.
+    # Therefore this operation is fast.
+    with timer('operator creation'):
+        A = dm.assembleNonlocal(kernel, matrixFormat='H2')
+
+with mgr.add('firstSolve'):
+    ######################################################################
+    # Set s = 0.75.
+    A.set(0.75)
+
+    # Let's solve a system.
+    # This triggers the assembly of the operators for the matrices at the interpolation nodes of the interval that contains s.
+    # The required matrices are constructed on-demand and then stay in memory.
+    with timer('solve 1 (slow)'):
+        solver = solverFactory('cg-jacobi', A=A, setup=True)
+        solver.maxIter = 1000
+        numIter = solver(b, u)
+    logger.info('Solved problem for s={} in {} iterations (residual norm {})'.format(A.get(), numIter, solver.residuals[-1]))
+
+with mgr.add('secondSolve'):
+    ######################################################################
+    # Let's solve a second sytem for a closeby value of s.
+    # This should be faster since we no longer need to assemble any matrices.
+    with timer('solve 2 (fast)'):
+        A.set(0.76)
+        solver = solverFactory('cg-jacobi', A=A, setup=True)
+        solver.maxIter = 1000
+        numIter = solver(b, u)
+    logger.info('Solved problem for s={} in {} iterations (residual norm {})'.format(A.get(), numIter, solver.residuals[-1]))
+
+with mgr.add('saveToFile'):
+    ######################################################################
+    # We can save the operator and the DoFMap to a file.
+    # This will trigger the assembly of all matrices.
+    with timer('save operator'):
+        h5_file = h5py.File('test.hdf5', 'w')
+        A.HDF5write(h5_file.create_group('A'))
+        dm.HDF5write(h5_file.create_group('dm'))
+        h5_file.close()
+
+with mgr.add('readFromFile'):
+    ######################################################################
+    # Now we can read them back in.
+    with timer('load operator'):
+        h5_file = h5py.File('test.hdf5', 'r')
+        A_2 = LinearOperator.HDF5read(h5_file['A'])
+        dm_2 = DoFMap.HDF5read(h5_file['dm'])
+        h5_file.close()
+
+with mgr.add('thirdSolve'):
+    # Set up and solve a system with the operator we loaded.
+    f_2 = functionFactory('constant', 2.)
+    b_2 = dm_2.assembleRHS(f_2)
+    u_2 = dm_2.zeros()
+    with timer('solve 3 (fast)'):
+        A_2.set(0.8)
+        solver = solverFactory('cg-jacobi', A=A_2, setup=True)
+        solver.maxIter = 1000
+        numIter = solver(b_2, u_2)
+    logger.info('Solved problem for s={} in {} iterations (residual norm {})'.format(A_2.get(), numIter, solver.residuals[-1]))
+    ######################################################################

docs/example3.rst

@@ -0,0 +1,82 @@
+
+Example 3 - Operator interpolation
+==================================
+
+This example demostrates the construction of a family of fractional
+Laplacians parametrized by the fractional order using operator
+interpolation. This can reduce the cost compared to assembling a new
+matrix for each value.
+
+The fractional Laplacian
+
+.. math::
+
+  (-\Delta)^{s} \text{ for } s \in [s_{\min}, s_{\max}] \subset (0, 1)
+
+is approximated by
+
+.. math::
+
+   (-\Delta)^{s} \approx \sum_{m=0}^{M} \Theta_{k,m}(s) (-\Delta)^{s_{k,m}} \text{ for } s \in \mathcal{S}_{k}
+
+for a sequence of intervals :math:`\mathcal{S}_{k}` that cover :math:`[s_{\min}, s_{\max}]` and scalar coefficients :math:`\Theta_{k,m}(s)`.
+The number of intervals and interpolation nodes is picked so that the interpolation error is dominated by the discretization error.
+
+The following example can be found at `examples/example3.py <https://github.com/sandialabs/PyNucleus/blob/master/examples/example3.py>`_.
+
+We set up a mesh, a dofmap and a fractional kernel.
+Instead of specifying a single value for the fractional order, we allow a range of values :math:`[s_{\min}, s_{\max}]=[0.05, 0.95]`.
+
+.. literalinclude:: ../examples/example3.py
+   :start-after: preamble
+   :end-before: #################
+   :lineno-match:
+
+Next, we call the assembly of a nonlocal operator as before.
+Since the operator for a particular value of the fractional order will be constructed on demand, this operation is fast.
+
+.. literalinclude:: ../examples/example3.py
+   :start-after: The operator is set up to be constructed on-demand
+   :end-before: #################
+   :lineno-match:
+
+.. program-output:: python3 example3.py --finalTarget operator
+
+Next, we choose the value of the fractional order. This needs to be within the range that we specified earlier.
+We then solve a linear system involving the operator.
+
+.. literalinclude:: ../examples/example3.py
+   :start-after: Set s = 0.75
+   :end-before: #################
+   :lineno-match:
+
+.. program-output:: python3 example3.py --finalTarget firstSolve
+
+This solve is relatively slow, as it involves the assembly of the nonlocal operators that are needed for the interpolation.
+We select a different value for the fractional order that is close to the first.
+Solving a linear system with this value is faster as we have already assembled the operator needed for the interpolation.
+
+.. literalinclude:: ../examples/example3.py
+   :start-after: This should be faster
+   :end-before: #################
+   :lineno-match:
+
+.. program-output:: python3 example3.py --finalTarget secondSolve
+
+Next, we save the operator to file.
+This first triggers the assembly of all operators nescessary to represent every value in :math:`s\in[0.05,0.95]`.
+
+.. literalinclude:: ../examples/example3.py
+   :start-after: We can save the operator
+   :end-before: #################
+   :lineno-match:
+
+.. program-output:: python3 example3.py --finalTarget saveToFile
+
+Next, we read the operator back in and solve another linear system.
+
+.. literalinclude:: ../examples/example3.py
+   :start-after: Now we can read them
+   :lineno-match:
+
+.. program-output:: python3 example3.py --finalTarget thirdSolve

docs/index.rst

@@ -36,6 +36,8 @@ Features
 
 * Nonlocal assembly (1D and 2D) into dense, sparse and hierarchical matrices
 
+* Operator interpolation
+
 * Solvers/preconditioners:
 
   * LU,
@@ -70,6 +72,15 @@ Examples
 
    example1
    example2
+   example3
+
+Drivers
+-------
+
+.. toctree::
+   :maxdepth: 1
+
+   drivers
 
 
 Funding

drivers/runNonlocal.py

@@ -32,6 +32,7 @@
 ##################################################
 
 vectors = d.addOutputGroup('vectors')
+vectors.add('dm', mS.u.dm)
 vectors.add('u', mS.u)
 if mS.u_interp is not None:
     vectors.add('uEx', mS.u_interp)