.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/example_Neumann.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_Neumann.py>`
to download the full example code.
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_example_Neumann.py:
Neumann condition for finite horizon kernel
==============================================================================
.. GENERATED FROM PYTHON SOURCE LINES 15-35
The following example can be found at `examples/example_Neumann.py <https://github.com/sandialabs/PyNucleus/blob/master/examples/example_Neumann.py>`_ in the PyNucleus repository.
.. math::
\gamma(x,y) = c(\delta) \chi_{B_\delta(x)}(y),
where :math:`c(\delta)` is the usual normalization constant.
.. math::
\int_{\Omega\cup\Omega_\mathcal{I}} (u(x)-u(y))\gamma(x,y) &= f(x) &&~in~ \Omega:=B_{1}(0),\\
\int_{\Omega\cup\Omega_\mathcal{I}} (u(x)-u(y))\gamma(x,y) &= g(x) &&~in~ \Omega_\mathcal{I}:=B_{1+\delta}(0)\setminus B_{1}(0),
where :math:`f\equiv 2` and
.. math::
g(x)=c(\delta) \left[ |x| \left((1+\delta-|x|)^2-\delta^2\right) + \frac{1}{3} \left((1+\delta-|x|)^3+\delta^3\right)\right].
The exact solution is :math:`u(x)=C-x^2` where :math:`C` is an arbitrary constant.
.. GENERATED FROM PYTHON SOURCE LINES 35-42
.. code-block:: Python
import numpy as np
import matplotlib.pyplot as plt
from PyNucleus import (nonlocalMeshFactory, dofmapFactory, kernelFactory,
functionFactory, solverFactory, NEUMANN, NO_BOUNDARY)
from PyNucleus_base.linear_operators import Dense_LinearOperator
.. GENERATED FROM PYTHON SOURCE LINES 43-44
Set up kernel, load $f$, the analytic solution and the flux data :math:`g`.
.. GENERATED FROM PYTHON SOURCE LINES 44-57
.. code-block:: Python
kernel = kernelFactory('constant', dim=1, horizon=0.4)
load = functionFactory('constant', 2.)
analyticSolution = functionFactory('Lambda', lambda x: -x[0]**2)
def fluxFun(x):
horizon = kernel.horizonValue
dist = 1+horizon-abs(x[0])
assert dist >= 0
return 2*kernel.scalingValue * (abs(x[0]) * (dist**2-horizon**2) + 1./3. * (dist**3+horizon**3))
flux = functionFactory('Lambda', fluxFun)
.. GENERATED FROM PYTHON SOURCE LINES 58-59
Construct a degree of freedom map for the entire mesh
.. GENERATED FROM PYTHON SOURCE LINES 59-67
.. code-block:: Python
mesh, nI = nonlocalMeshFactory('interval', kernel=kernel, boundaryCondition=NEUMANN)
for _ in range(3):
mesh = mesh.refine()
dm = dofmapFactory('P1', mesh, NO_BOUNDARY)
dm
.. rst-class:: sphx-glr-script-out
.. code-block:: none
P1 DoFMap with 57 DoFs and 0 boundary DoFs.
.. GENERATED FROM PYTHON SOURCE LINES 68-69
The second return value of the nonlocal mesh factory contains indicator functions:
.. GENERATED FROM PYTHON SOURCE LINES 69-75
.. code-block:: Python
dm.interpolate(nI['domain']).plot(label='domain')
dm.interpolate(nI['boundary']).plot(label='local boundary')
dm.interpolate(nI['interaction']).plot(label='interaction')
plt.legend()
.. image-sg:: /auto_examples/images/sphx_glr_example_Neumann_001.png
:alt: example Neumann
:srcset: /auto_examples/images/sphx_glr_example_Neumann_001.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
<matplotlib.legend.Legend object at 0x7fbee00cf380>
.. GENERATED FROM PYTHON SOURCE LINES 76-77
Assemble the RHS vector by using the load on the domain :math:`\Omega` and the flux function on the interaction domain :math:`\Omega_\mathcal{I}`
.. GENERATED FROM PYTHON SOURCE LINES 77-81
.. code-block:: Python
A = dm.assembleNonlocal(kernel)
b = dm.assembleRHS(load*nI['domain'] + flux*(nI['interaction']+nI['boundary']))
.. GENERATED FROM PYTHON SOURCE LINES 82-83
Solve the linear system. Since it is singular (shifts by a constant form the nullspace) we augment the system and then project out the zero mode.
.. GENERATED FROM PYTHON SOURCE LINES 83-86
.. code-block:: Python
u = dm.zeros()
.. GENERATED FROM PYTHON SOURCE LINES 87-88
Augment the system
.. GENERATED FROM PYTHON SOURCE LINES 88-92
.. code-block:: Python
correction = Dense_LinearOperator(np.ones(A.shape))
solver = solverFactory('lu', A=A+correction, setup=True)
solver(b, u)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
1
.. GENERATED FROM PYTHON SOURCE LINES 93-94
project out the nullspace component
.. GENERATED FROM PYTHON SOURCE LINES 94-97
.. code-block:: Python
const = dm.ones()
u = u - (u.inner(const)/const.inner(const))*const
.. GENERATED FROM PYTHON SOURCE LINES 98-99
Interpolate the exact solution and project out zero mode as well.
.. GENERATED FROM PYTHON SOURCE LINES 99-103
.. code-block:: Python
uex = dm.interpolate(analyticSolution)
uex = uex - (uex.inner(const)/const.inner(const))*const
.. GENERATED FROM PYTHON SOURCE LINES 104-107
.. code-block:: Python
u.plot(label='numerical', marker='x')
uex.plot(label='analytic')
plt.legend()
.. image-sg:: /auto_examples/images/sphx_glr_example_Neumann_002.png
:alt: example Neumann
:srcset: /auto_examples/images/sphx_glr_example_Neumann_002.png
:class: sphx-glr-single-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
<matplotlib.legend.Legend object at 0x7fbedfcdc910>
.. rst-class:: sphx-glr-timing
**Total running time of the script:** (0 minutes 0.132 seconds)
.. _sphx_glr_download_auto_examples_example_Neumann.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: example_Neumann.ipynb <example_Neumann.ipynb>`
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: example_Neumann.py <example_Neumann.py>`
.. container:: sphx-glr-download sphx-glr-download-zip
:download:`Download zipped: example_Neumann.zip <example_Neumann.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_