Tricks section of the documentation (#46)

Author: jodemaey

README.md

@@ -4,6 +4,8 @@
 >
 > -- *Michael Gambon - Layer Cake (2004)*
 
+[![Documentation Status](https://img.shields.io/badge/docs-passing-green.svg)](https://climdyn.github.io/LayerCake/)
+
 ## General Information
 
 LayerCake is a framework to design models based on systems of partial differential equations (PDEs), 
@@ -69,5 +71,5 @@ developed.
 
 ## Contributing
 
-LayerCake is in betaa development phase, bug reports and tests of the features are welcome.
+LayerCake is in beta development phase, bug reports and tests of the features are welcome.
 Please simply raise an issue on [Github](https://github.com/Climdyn/LayerCake/issues).

documentation/source/files/general_information.rst

@@ -68,8 +68,8 @@ Examples
 A few examples are available in the `examples <../../../../examples>`_ folder. More examples will be provided as the code is
 developed.
 
-Forthcoming development
------------------------
+Contributing
+------------
 
-LayerCake is in alpha development phase, with many new functionalities planned over the next few months.
-Stay tuned !
+LayerCake is in beta development phase, bug reports and tests of the features are welcome.
+Please simply raise an issue on `Github <https://github.com/Climdyn/LayerCake/issues>`_.

documentation/source/files/technical/utils.rst

@@ -19,3 +19,6 @@ Module with various submodules with useful functions used throughout the code.
 
 .. automodule:: layercake.utils.integration
     :members:
+
+.. automodule:: layercake.utils.matrix
+    :members:

documentation/source/files/tricks.rst

@@ -0,0 +1,192 @@
+
+Various tricks
+==============
+
+On this page, we list several tricks on how to do this or that, related to particular aspects
+of the modelling based on systems of partial differential equations.
+
+Spatially varying parameters
+----------------------------
+
+It is possible in LayerCake to add parameters varying as a function of the coordinates.
+To illustrate this, we are going to modify a spatially uniform parameter in one of the
+`examples <../../../../examples/>`_, namely the
+`Reinhold & Pierrehumber model <https://github.com/Climdyn/LayerCake/blob/main/examples/atmospheric/baroclinic_one_layer.py>`_ one.
+
+In this model, two terms are representing the friction of the bottom layer with the ground:
+
+* :math:`-\frac{k_d}{2} \nabla^2 (\psi - \theta)` for the equation involving the barotropic streamfunction :math:`\psi`
+* :math:`+\frac{k_d}{2} \nabla^2 (\psi - \theta)` for the equation involving the baroclinic streamfunction :math:`\theta`
+
+where the parameter controlling the friction is :math:`k_d`. The full equations of the model can be found `here <https://qgs.readthedocs.io/en/latest/files/model/oro_model.html#mid-layer-equations-and-the-thermal-wind-relation>`_.
+
+Running the `LayerCake code <../../../../examples/atmospheric/baroclinic_one_layer.py>`_  produces the following figure:
+
+.. figure:: tricks/RP1982.png
+    :scale: 100%
+    :align: center
+
+which is a section of the model's attractor in the :math:`\psi_2, \psi_3` plane.
+
+Now we can modify the original model's equations by making the parameter :math:`k_d` varies as
+
+.. math::
+
+    k_{d,0} + 2 \, k_{d,2} \cos(n x) \sin(y)
+
+This can be done by using the :class:`~layercake.arithmetic.terms.operations.ProductOfTerms` objects to do the product of
+a :class:`~layercake.arithmetic.terms.linear.LinearTerm` term involving a :class:`~layercake.variables.field.ParameterField` field representing
+the spatial variation of the parameter, and the Laplacian terms representing :math:`\nabla^2 (\psi - \theta)`.
+It can be implemented with a few lines added to the `Reinhold & Pierrehumber model code <https://github.com/Climdyn/LayerCake/blob/main/examples/atmospheric/baroclinic_one_layer.py>`_:
+
+.. code:: ipython3
+
+    dfriction = ProductOfTerms(LinearTerm(Dk), OperatorTerm(psi, Laplacian, atmospheric_basis.coordinate_system, sign=-1))
+    dofriction = ProductOfTerms(LinearTerm(Dk), OperatorTerm(theta, Laplacian, atmospheric_basis.coordinate_system))
+    barotropic_equation.add_rhs_terms([dfriction, dofriction])
+
+in the barotropic equation and
+
+.. code:: ipython3
+
+    dfriction = ProductOfTerms(LinearTerm(Dk), OperatorTerm(psi, Laplacian, atmospheric_basis.coordinate_system, sign=1))
+    dofriction = ProductOfTerms(LinearTerm(Dk), OperatorTerm(theta, Laplacian, atmospheric_basis.coordinate_system, sign=-1))
+    baroclinic_equation.add_rhs_terms([dfriction, dofriction])
+
+in the baroclinic equation, where :code:`Dk` is the :class:`~layercake.variables.field.ParameterField` object defining the spatial
+variation of :math:`k_d`
+
+.. code:: ipython3
+
+    # Variable bottom friction
+    dk = np.zeros(len(atmospheric_basis))
+    dk[1] = 0.1 * kdp_deriv
+    Dk = ParameterField('D_k', u'D_k', dk, atmospheric_basis, inner_products_definition)
+
+with :math:`k_{d,2} = D_{k,1} = 0.1 \, k_d` (note the difference of index because of the Python-specific indexing starting from zero).
+The equations LaTeX representation now clearly shows the new terms:
+
+.. figure:: tricks/mod_eqs.png
+    :align: center
+
+with the appearance of the :math:`D_k \nabla^2` terms.
+
+The impact of this spatial variation of the bottom friction coefficient on the model's dynamics
+is clearly visible on the 2-dimensional section of the attractor:
+
+.. figure:: tricks/RP1982mod.png
+    :scale: 100%
+    :align: center
+
+(Compare with the first figure above.)
+
+Usage of mathematical expressions in the PDEs
+---------------------------------------------
+
+Mathematical function can appear in some models PDEs. In general, these are functions of the model's domain coordinates.
+Modelling terms involving mathematical functions (of the coordinates) can be implemented in LayerCake using the
+:class:`~layercake.arithmetic.symbolic.expressions.Expression` class. In this section, we recycle the example of the
+previous one by modifying again the friction with the bottom layer in
+the `Reinhold & Pierrehumber model <https://github.com/Climdyn/LayerCake/blob/main/examples/atmospheric/baroclinic_one_layer.py>`_:
+
+.. math::
+
+    k_{d} \to k_{d} + F(y) = k_d  + k_{d,y} \, \sin(y)
+
+which means that the friction is higher in the higher latitude.
+
+To implement this, we simply need to create an :class:`~layercake.arithmetic.symbolic.expressions.Expression` object to
+represent :math:`F(y)` as follow:
+
+.. code:: ipython3
+
+    # Variable bottom friction
+    from sympy import sin
+    kd_deriv_y = Parameter(0.1 * kd_deriv, symbol=Symbol('k_dy'), latex='k_{d,y}')
+    F = Expression(kd_deriv_y.symbol * sin(y), expression_parameters=[kd_deriv_y,], latex=r'k_{d,y} \sin(y)')
+
+where we set :math:`k_{d,y} = 0.1 \, k_d`.
+Then the variable friction with the bottom can be added like in the previous section, involving the expression
+:code:`F` as prefactor:
+
+.. code:: ipython3
+
+    dfriction = OperatorTerm(psi, Laplacian, atmospheric_basis.coordinate_system, prefactor=F, sign=-1)
+    dofriction = OperatorTerm(theta, Laplacian, atmospheric_basis.coordinate_system, prefactor=F)
+    barotropic_equation.add_rhs_terms([dfriction, dofriction])
+
+    dfriction = OperatorTerm(psi, Laplacian, atmospheric_basis.coordinate_system, prefactor=F, sign=1)
+    dofriction = OperatorTerm(theta, Laplacian, atmospheric_basis.coordinate_system, prefactor=F, sign=-1)
+    baroclinic_equation.add_rhs_terms([dfriction, dofriction])
+
+and this result in the following equations in LayeCake:
+
+.. figure:: tricks/mod_eqs_siny.png
+    :align: center
+
+where the terms involving :math:`\sin(y)` are appearing.
+
+Constructing this model and then integrating it then yield:
+
+.. figure:: tricks/RP1982mod_siny.png
+    :scale: 100%
+    :align: center
+
+(Again, you can compare with the first figure above.)
+
+Free-threading
+--------------
+
+Since Python 3.14, a `free-threaded version of Python is available <https://docs.python.org/3/howto/free-threading-python.html>`_.
+LayerCake has been tested for it and seems so far to run smoothly.
+
+This feature is particularly useful if you try to derive complicate, big models with a
+large number of modes. LayerCake will then use multiple threads to perform |Sympy| symbolic evaluations,
+while the integration of the inner products will still be done using multiple processes.
+When being used in a free-threading environment, this is the default behavior, but this can be controlled by environment
+variables:
+
+* The :code:`LAYERCAKE_PARALLEL_METHOD` environment variable defines how the Sympy symbolic evaluations are done. It can take two different values:
+    + :code:`threads`: the evaluations will be done using threads
+    + :code:`processes`: the evaluations will be done using processes. In some complicate cases, it might lead to wrong answers or even crash. This mode is thus not recommended unless you know what you are doing.
+
+  If this environment variable is not defined, then LayerCake default behavior is to use threads.
+* The :code:`LAYERCAKE_PARALLEL_INTEGRATION` environment variable controls the Sympy symbolic integration parallelization. If set to :code:`none`, the parallelization will be deactivated.
+  Otherwise, it will parallelized using processes.
+  If this environment variable is not defined, then LayerCake default behavior is to parallelize using processes.
+
+For example,
+
+.. code:: bash
+
+    LAYERCAKE_PARALLEL_METHOD=processes python examples/atmospheric/barotropic_one_layer.py
+
+will launch the `barotropic one layer <https://github.com/Climdyn/LayerCake/blob/main/examples/atmospheric/barotropic_one_layer.py>`_
+script with Sympy symbolic evaluation being done using processes.
+
+.. warning::
+
+    Launching the previous line with symbolic evaluation being done in processes with a free-threaded Python will result in
+    a crash. Use the free-threaded Python only if :code:`LAYERCAKE_PARALLEL_METHOD` is set to :code:`threads`.
+
+Installing the free-threaded version of Python can be done using Anaconda, by typing:
+
+.. code:: bash
+
+    conda env create -f environment_freethreading.yml
+
+which will create a conda environment :code:`layercake_ft`.
+Upon activation:
+
+.. code:: bash
+
+    conda activate layercake_ft
+
+any Python code will be run in `free-threading` mode.
+
+.. warning::
+
+    Python free-threading mode is still somewhat experimental, and
+    the obtained models and results must be
+    scrutinized with care and double-checked.
+

documentation/source/files/tricks/RP1982.png

@@ -1,13 +1,13 @@
 User guide
 ==========
 
-This guide explains how the LayerCake framework can be used to transform a set of two-dimensional partial differential
+This guide explains how the LayerCake framework can be used to transform a set of two-dimensional partial differential equations (PDEs)
 
 .. math::
 
-    \partial_t \mathcal{F}^{\mathrm LHS}_i (\psi_1, \ldots, \psi_N) = \mathcal{F}^{\mathrm{RHS}}_i (\psi_1, \ldots, \psi_N) \qquad , \quad i = 1,\ldots,N
+    \partial_t \mathcal{F}^{\mathrm LHS}_i \left[\psi_1, \ldots, \psi_N\right] = \mathcal{F}^{\mathrm{RHS}}_i \left[\psi_1, \ldots, \psi_N\right] \qquad , \quad i = 1,\ldots,N
 
-equations (PDEs) defined on a particular domain into a system of ordinary differential equations (ODEs)
+defined on a particular domain into a system of ordinary differential equations (ODEs)
 with an automated `Galerkin method`_. This method projects all the fields :math:`\psi_j` on given function basis :math:`\phi_{j,k}`:
 
 .. math::
@@ -197,7 +197,7 @@ If one type :code:`atmospheric_basis` in a terminal after defining it way above,
     [sqrt(2)*cos(y), 2*sin(y)*cos(n*x), 2*sin(y)*sin(n*x), sqrt(2)*cos(2*y), 2*sin(2*y)*cos(n*x), 2*sin(2*y)*sin(n*x), 2*sin(y)*cos(2*n*x), 2*sin(y)*sin(2*n*x), 2*sin(2*y)*cos(2*n*x), 2*sin(2*y)*sin(2*n*x)]
 
 We also need to create an inner product definition so that LayerCake knows how you want your PDEs to be projected.
-In general, using the :class:`~layercake.inner_products.definition.StandardSymbolicInnerProductDefinition` definition is sufficient for most model:
+In general, using the :class:`~layercake.inner_products.definition.StandardSymbolicInnerProductDefinition` definition is sufficient for most models:
 
 .. code:: ipython3
 
@@ -809,6 +809,11 @@ This functionality is thus particularly useful to get the tendencies and Jacobia
 involving the model's parameters. For example, this enables the study of the models with bifurcation analysis
 tools such as AUTO-07p :cite:`user-doedel2007` and auto-AUTO :cite:`user-DH2025`.
 
+.. note::
+
+    As a final note, you can also take a look at the other examples in the `examples <../../../examples/>`_ folder.
+    You can also find specific LayerCake tricks on the :ref:`files/tricks:Various tricks` page.
+
 
 References
 ----------