.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/example_operator_interpolation.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_example_operator_interpolation.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_example_operator_interpolation.py:


Operator interpolation
======================

.. GENERATED FROM PYTHON SOURCE LINES 15-36

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/example_operator_interpolation.py <https://github.com/sandialabs/PyNucleus/blob/master/examples/example_operator_interpolation.py>`_ in the PyNucleus repository.

.. GENERATED FROM PYTHON SOURCE LINES 36-50

.. code-block:: Python


    import logging
    from PyNucleus_base.utilsFem import TimerManager
    import numpy as np
    import matplotlib.pyplot as plt

    fmt = '{message}'
    logging.basicConfig(level=logging.INFO,
                        format=fmt,
                        style='{',
                        datefmt="%Y-%m-%d %H:%M:%S")
    logger = logging.getLogger('__main__')
    timer = TimerManager(logger)








.. GENERATED FROM PYTHON SOURCE LINES 51-53

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]`.

.. GENERATED FROM PYTHON SOURCE LINES 53-70

.. code-block:: Python


    from PyNucleus import meshFactory, dofmapFactory, kernelFactory, functionFactory, solverFactory
    from PyNucleus_nl.operatorInterpolation import admissibleSet

    # 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)








.. GENERATED FROM PYTHON SOURCE LINES 71-75

Next, we call the assembly of a nonlocal operator as before.
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.

.. GENERATED FROM PYTHON SOURCE LINES 75-78

.. code-block:: Python

    with timer('operator creation'):
        A = dm.assembleNonlocal(kernel, matrixFormat='H2')





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    operator creation in 0.126 s




.. GENERATED FROM PYTHON SOURCE LINES 79-81

Next, we choose the value of the fractional order. This needs to be within the range that we specified earlier.
We set s = 0.75.

.. GENERATED FROM PYTHON SOURCE LINES 81-83

.. code-block:: Python

    A.set(0.75)








.. GENERATED FROM PYTHON SOURCE LINES 84-87

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.

.. GENERATED FROM PYTHON SOURCE LINES 87-93

.. code-block:: Python

    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]))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    solve 1 (slow) in 0.658 s
    Solved problem for s=0.75 in 181 iterations (residual norm 9.197705043865965e-06)




.. GENERATED FROM PYTHON SOURCE LINES 94-97

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.

.. GENERATED FROM PYTHON SOURCE LINES 97-105

.. code-block:: Python


    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]))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    solve 2 (fast) in 0.14 s
    Solved problem for s=0.76 in 190 iterations (residual norm 8.897510981041587e-06)




.. GENERATED FROM PYTHON SOURCE LINES 106-108

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]`.

.. GENERATED FROM PYTHON SOURCE LINES 108-117

.. code-block:: Python


    import h5py

    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()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    save operator in 19.8 s




.. GENERATED FROM PYTHON SOURCE LINES 118-119

Next, we read the operator back in.

.. GENERATED FROM PYTHON SOURCE LINES 119-129

.. code-block:: Python


    from PyNucleus_base.linear_operators import LinearOperator
    from PyNucleus_fem.DoFMaps import DoFMap

    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()





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    load operator in 6.65 s




.. GENERATED FROM PYTHON SOURCE LINES 130-131

Finally, we set up and solve a series of linear systems with the operator we loaded.

.. GENERATED FROM PYTHON SOURCE LINES 131-144

.. code-block:: Python

    f_2 = functionFactory('constant', 2.)
    b_2 = dm_2.assembleRHS(f_2)
    for sVal in np.linspace(0.1, 0.9, 9):
        with timer('solve 3 (fast)'):
            u_2 = dm_2.zeros()
            A_2.set(sVal)
            solver = solverFactory('cg-jacobi', A=A_2, setup=True)
            solver.maxIter = 1000
            numIter = solver(b_2, u_2)
        logger.info('Solved problem for s={:.1} in {} iterations (residual norm {})'.format(A_2.get(), numIter, solver.residuals[-1]))
        u_2.plot(label='s={:.1}'.format(sVal))
    plt.legend()




.. image-sg:: /auto_examples/images/sphx_glr_example_operator_interpolation_001.png
   :alt: example operator interpolation
   :srcset: /auto_examples/images/sphx_glr_example_operator_interpolation_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    solve 3 (fast) in 0.00681 s
    Solved problem for s=0.1 in 7 iterations (residual norm 3.3399038658481327e-06)
    solve 3 (fast) in 0.0114 s
    Solved problem for s=0.2 in 13 iterations (residual norm 3.6616774310675407e-06)
    solve 3 (fast) in 0.017 s
    Solved problem for s=0.3 in 21 iterations (residual norm 7.068770471207542e-06)
    solve 3 (fast) in 0.0264 s
    Solved problem for s=0.4 in 34 iterations (residual norm 8.650357989614637e-06)
    solve 3 (fast) in 0.0422 s
    Solved problem for s=0.5 in 55 iterations (residual norm 9.299204350662863e-06)
    solve 3 (fast) in 0.0666 s
    Solved problem for s=0.6 in 90 iterations (residual norm 8.670054190291531e-06)
    solve 3 (fast) in 0.116 s
    Solved problem for s=0.7 in 146 iterations (residual norm 9.855504121491694e-06)
    solve 3 (fast) in 0.175 s
    Solved problem for s=0.8 in 234 iterations (residual norm 8.64786773747414e-06)
    solve 3 (fast) in 0.271 s
    Solved problem for s=0.9 in 361 iterations (residual norm 9.69089344852077e-06)

    <matplotlib.legend.Legend object at 0x7fbee01125d0>



.. GENERATED FROM PYTHON SOURCE LINES 145-146

.. code-block:: Python

    plt.show()








.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 28.238 seconds)


.. _sphx_glr_download_auto_examples_example_operator_interpolation.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: example_operator_interpolation.ipynb <example_operator_interpolation.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: example_operator_interpolation.py <example_operator_interpolation.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: example_operator_interpolation.zip <example_operator_interpolation.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_