.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/example_InfHorizonDirichlet.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_InfHorizonDirichlet.py>`
to download the full example code.
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_example_InfHorizonDirichlet.py:
Dirichlet condition for infinite horizon kernel
==================================================================================
.. GENERATED FROM PYTHON SOURCE LINES 15-34
The following example can be found at `examples/example_InfHorizonDirichlet.py <https://github.com/sandialabs/PyNucleus/blob/master/examples/example_InfHorizonDirichlet.py>`_ in the PyNucleus repository.
We want to solve a problem with infinite horizon fractional kernel
.. math::
\gamma(x,y) = \frac{c}{|x-y|^{d+2s}},
where :math:`c` is the usual normalization constant, and imposes an inhomogeneous Dirichlet volume condition.
We will have to make some assumption on the volume condition in order to make it computationally tractable.
Here we assume that it is zero at some distance away from the domain.
.. math::
(-\Delta)^s u &= f &&~in~ \Omega:=B_{1/2}(0),\\
u &= g &&~in~ \Omega_\mathcal{I}:=B_{1}(0)\setminus B_{1/2}(0), \\
u &= 0 &&~in~ \mathbb{R}^d \setminus B_{1}(0),
where :math:`f\equiv 1` and :math:`g` is chosen to match the known exact solution :math:`u(x)=C(1-|x|^2)_+^s` with some constant :math:`C`.
.. GENERATED FROM PYTHON SOURCE LINES 34-39
.. code-block:: Python
import numpy as np
import matplotlib.pyplot as plt
from PyNucleus import meshFactory, dofmapFactory, kernelFactory, functionFactory, solverFactory
.. GENERATED FROM PYTHON SOURCE LINES 40-41
Construct a mesh for :math:`\Omega\cup \Omega_\mathcal{I}` and set up degree of freedom maps on :math:`\Omega` and :math:`\Omega_\mathcal{I}`.
.. GENERATED FROM PYTHON SOURCE LINES 41-47
.. code-block:: Python
radius = 1.0
mesh = meshFactory('disc', n=8, radius=radius)
for _ in range(4):
mesh = mesh.refine()
.. GENERATED FROM PYTHON SOURCE LINES 48-50
Get an indicator function for :math:`\Omega`.
Subtract 1e-6 to avoid roundoff errors for nodes that are exactly on :math:`|x| = 0.5`.
.. GENERATED FROM PYTHON SOURCE LINES 50-59
.. code-block:: Python
OmegaIndicator = functionFactory('radialIndicator', 0.5*radius-1e-6)
dm = dofmapFactory('P1', mesh, OmegaIndicator)
dmBC = dm.getComplementDoFMap()
plt.figure()
dm.plot(printDoFIndices=False)
plt.figure()
dmBC.plot(printDoFIndices=False)
.. rst-class:: sphx-glr-horizontal
*
.. image-sg:: /auto_examples/images/sphx_glr_example_InfHorizonDirichlet_001.png
:alt: example InfHorizonDirichlet
:srcset: /auto_examples/images/sphx_glr_example_InfHorizonDirichlet_001.png
:class: sphx-glr-multi-img
*
.. image-sg:: /auto_examples/images/sphx_glr_example_InfHorizonDirichlet_002.png
:alt: example InfHorizonDirichlet
:srcset: /auto_examples/images/sphx_glr_example_InfHorizonDirichlet_002.png
:class: sphx-glr-multi-img
.. GENERATED FROM PYTHON SOURCE LINES 60-61
Set up the kernel, rhs and known solution.
.. GENERATED FROM PYTHON SOURCE LINES 61-67
.. code-block:: Python
s = 0.75
kernel = kernelFactory('fractional', dim=mesh.dim, s=s, horizon=np.inf)
rhs = functionFactory('constant', 1.)
uex = functionFactory('solFractional', dim=mesh.dim, s=s, radius=radius)
.. GENERATED FROM PYTHON SOURCE LINES 68-76
Assemble the linear system
.. math::
A_{\Omega,\Omega} u_{\Omega} + A_{\Omega,\Omega_\mathcal{I}} u_{\Omega_\mathcal{I}} &= f \\
u_{\Omega_\mathcal{I}} &= g
Set up a solver for :math:`A_{\Omega,\Omega}`.
.. GENERATED FROM PYTHON SOURCE LINES 76-83
.. code-block:: Python
A_OmegaOmega = dm.assembleNonlocal(kernel, matrixFormat='H2')
A_OmegaOmegaI = dm.assembleNonlocal(kernel, dm2=dmBC)
f = dm.assembleRHS(rhs)
g = g = dmBC.interpolate(uex)
solver = solverFactory('lu', A=A_OmegaOmega, setup=True)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Efficiency of assembly with 2 DoFMaps is bad.
.. GENERATED FROM PYTHON SOURCE LINES 84-85
Compute FE solution and error between FE solution and interpolation of analytic expression
.. GENERATED FROM PYTHON SOURCE LINES 85-91
.. code-block:: Python
u_Omega = dm.zeros()
solver(f-A_OmegaOmegaI*g, u_Omega)
u = u_Omega.augmentWithBoundaryData(g)
err = u.dm.interpolate(uex) - u
.. GENERATED FROM PYTHON SOURCE LINES 92-93
Plot FE solution and error
.. GENERATED FROM PYTHON SOURCE LINES 93-100
.. code-block:: Python
plt.figure()
plt.title('FE solution')
u.plot(flat=True)
plt.figure()
plt.title('error')
err.plot(flat=True)
.. rst-class:: sphx-glr-horizontal
*
.. image-sg:: /auto_examples/images/sphx_glr_example_InfHorizonDirichlet_003.png
:alt: FE solution
:srcset: /auto_examples/images/sphx_glr_example_InfHorizonDirichlet_003.png
:class: sphx-glr-multi-img
*
.. image-sg:: /auto_examples/images/sphx_glr_example_InfHorizonDirichlet_004.png
:alt: error
:srcset: /auto_examples/images/sphx_glr_example_InfHorizonDirichlet_004.png
:class: sphx-glr-multi-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
<matplotlib.collections.PolyCollection object at 0x7fbedb9fda90>
.. rst-class:: sphx-glr-timing
**Total running time of the script:** (0 minutes 10.188 seconds)
.. _sphx_glr_download_auto_examples_example_InfHorizonDirichlet.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: example_InfHorizonDirichlet.ipynb <example_InfHorizonDirichlet.ipynb>`
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: example_InfHorizonDirichlet.py <example_InfHorizonDirichlet.py>`
.. container:: sphx-glr-download sphx-glr-download-zip
:download:`Download zipped: example_InfHorizonDirichlet.zip <example_InfHorizonDirichlet.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_