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

.. only:: html

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

        :ref:`Go to the end <sphx_glr_download_api_gallery_plate_tectonics_plot_03_kelvin_cooling_earth_grid.py>`
        to download the full example code or to run this example in your browser via JupyterLite.

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

.. _sphx_glr_api_gallery_plate_tectonics_plot_03_kelvin_cooling_earth_grid.py:


Kelvin's conduction on a discrete grid
======================================

Kelvin's calculation rests on Fourier's conduction law, under which every
sinusoidal temperature mode decays at the rate ``kappa k**2``. On a periodic
grid the same mode decays at ``4 kappa sin(k dx / 2)**2 / dx**2``, which tends
to Fourier's rate as the grid is refined, and the mean temperature is
conserved because no heat leaves a periodic domain.

.. GENERATED FROM PYTHON SOURCE LINES 11-39



.. image-sg:: /api/gallery/plate_tectonics/images/sphx_glr_plot_03_kelvin_cooling_earth_grid_001.png
   :alt: Decay of one Fourier mode
   :srcset: /api/gallery/plate_tectonics/images/sphx_glr_plot_03_kelvin_cooling_earth_grid_001.png
   :class: sphx-glr-single-img





.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np

    from earthsciencekit.plate_tectonics import diffuse_temperature
    from earthsciencekit.units import convert

    n, spacing, diffusivity = 64, 1000.0, 1e-6
    x = np.arange(n) * spacing
    k = 2 * np.pi / (n * spacing)
    initial = 900 + 200 * np.cos(k * x)
    time = convert(np.linspace(0, 0.6e6, 101), "year", "s")
    result = diffuse_temperature(initial, time, diffusivity, spacing)
    np.testing.assert_allclose(result.values.mean(axis=1), 900, atol=1e-10)
    amplitude = (result.values[:, 0] - 900) / 200
    assert np.isclose(-np.log(amplitude[-1]) / time[-1], diffusivity * k**2, rtol=1e-3)

    fig, (left, right) = plt.subplots(1, 2, figsize=(9, 4))
    image = left.pcolormesh(x / 1e3, convert(time, "s", "year") / 1e3, result.values, shading="auto", cmap="inferno")
    left.set(xlabel="periodic position [km]", ylabel="elapsed time [thousand years]")
    fig.colorbar(image, ax=left, label="temperature [K]")
    right.plot(convert(time, "s", "year") / 1e3, amplitude, label="grid")
    right.plot(convert(time, "s", "year") / 1e3, np.exp(-diffusivity * k**2 * time), "--", label="Fourier rate")
    right.set(xlabel="elapsed time [thousand years]", ylabel="mode amplitude", title="Decay of one Fourier mode")
    right.legend()
    fig.tight_layout()

    plt.show()


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

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


.. _sphx_glr_download_api_gallery_plate_tectonics_plot_03_kelvin_cooling_earth_grid.py:

.. only:: html

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

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../../lite/lab/index.html?path=api/gallery/plate_tectonics/plot_03_kelvin_cooling_earth_grid.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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