From 4a232d7e8a5d8b932618680eeae8836eae165948 Mon Sep 17 00:00:00 2001 From: Jonathan Demaeyer Date: Sun, 7 Mar 2021 12:12:40 +0100 Subject: [PATCH] First release of the user guide --- README.md | 28 +- basis/base.py | 4 + basis/fourier.py | 2 +- .../source/files/general_information.rst | 34 +- .../source/files/model/atmosphere.rst | 10 +- documentation/source/files/model/li_model.rst | 24 +- .../source/files/model/maooam_model.rst | 30 +- documentation/source/files/model/ocean.rst | 11 +- .../source/files/model/oro_model.rst | 22 +- .../source/files/model_description.rst | 5 +- .../source/files/technical/configuration.rst | 4 +- .../source/files/technical_description.rst | 4 +- documentation/source/files/user_guide.rst | 515 +++++++++++++++++- documentation/source/index.rst | 7 + functions/sparse_mul.py | 10 +- functions/tendencies.py | 2 +- inner_products/analytic.py | 12 +- inner_products/base.py | 2 +- inner_products/definition.py | 2 +- inner_products/symbolic.py | 135 ++++- params/params.py | 16 +- 21 files changed, 736 insertions(+), 143 deletions(-) diff --git a/README.md b/README.md index 8a4e3c1..d795462 100644 --- a/README.md +++ b/README.md @@ -93,7 +93,7 @@ Usage qgs can be used by editing and running the script `qgs_rp.py` and `qgs_maooam.py` found in the main folder. -For more advanced usages, please read the User Guides (TODO). +For more advanced usages, please read the [User Guides](https://qgs.readthedocs.io/en/latest/files/user_guide.html). Examples -------- @@ -130,36 +130,22 @@ Forthcoming developments * Technical mid-term developments + Dimensionally robust Parameter class operation + Windows OS support - + Symbolic inner products (using e.g. [Sympy](https://www.sympy.org/)) - - Arbitrary spatial mode basis of functions - - Automatic on-the-fly inner product calculation (numeric or analytic if possible) - - Symbolic PDE equation specification + + Numerical basis of function + Visualisation tools, e.g. based on the [movie-script](https://github.com/jodemaey/movie-script) package * Long-term development track + Active advection + True quasi-geostrophic ocean when using ocean model version - + Salinity in the ocean - + + Salinity in the ocean + + Symbolic PDE equation specification + Contributing to qgs ------------------- -If you want to contribute actively to the roadmap detailed above, please contact the authors. +If you want to contribute actively to the roadmap detailed above, please contact the main authors. In addition, if you have made changes that you think will be useful to others, please feel free to suggest these as a pull request on the [qgs Github repository](https://github.com/Climdyn/qgs). -A review of your pull request will follow with possibly suggestions of changes before merging it in the master branch. -Please consider the following guidelines before submitting: -* Before submitting a pull request, double check that the branch to be merged contains only changes you wish to add to the master branch. This will save time in reviewing the code. -* For any changes to the core model files, please run the tests found in the folder [model_test](./model_test) to ensure that the model tensors are still valid. -* For substantial additions of code, including a test case in the folder [model_test](./model_test) is recommended. -* Please do not make changes to existing test cases. -* Please document the new functionalities in the documentation. Code addition without documentation addition will not be accepted. -The documentation is done with [sphinx](https://www.sphinx-doc.org/en/master/) and follows the Numpy conventions. Please take a look to the actual code to get an idea about how to document the code. -* If your addition can be considered as a tool not directly related to the core of the model, please develop it in the toolbox folder. -* The team presently maintaining qgs is not working full-time on it, so please be patient as the review of the submission may take some time. - -For more information about git, Github and the pull request framework, a good source of information is the [contributing guide](https://mitgcm.readthedocs.io/en/latest/contributing/contributing.html) of the [MITgcm](https://github.com/MITgcm/MITgcm). - +More information and guidance about how to do a pull request for qgs can be found in the documentation [here](https://qgs.readthedocs.io/en/latest/files/general_information.html#contributing-to-qgs). Other atmospheric models in Python ---------------------------------- diff --git a/basis/base.py b/basis/base.py index 1dc9fdf..26b41be 100644 --- a/basis/base.py +++ b/basis/base.py @@ -19,6 +19,10 @@ .. _abstract base class: https://docs.python.org/3/glossary.html#term-abstract-base-class """ + +# TODO: define setters and init arguments +# TODO: define NumericBasis class + import sys from abc import ABC diff --git a/basis/fourier.py b/basis/fourier.py index c1cd711..902ada2 100644 --- a/basis/fourier.py +++ b/basis/fourier.py @@ -21,7 +21,7 @@ from sympy import symbols, sin, cos, sqrt _x, _y = symbols('x y') -_n = symbols('n') +_n = symbols('n', real=True, nonnegative=True) class ChannelFourierBasis(SymbolicBasis): diff --git a/documentation/source/files/general_information.rst b/documentation/source/files/general_information.rst index 41f3115..735e2ae 100644 --- a/documentation/source/files/general_information.rst +++ b/documentation/source/files/general_information.rst @@ -87,9 +87,7 @@ or :: python qgs_maooam.py -For more advanced usages, please read the User Guide. - -.. TODO: add more details and a link to it when the User Guide is ready. +For more advanced usages, please read the :ref:`files/user_guide:User guide`. Examples -------- @@ -129,12 +127,7 @@ Forthcoming developments + Dimensionally robust Parameter class operation + Windows OS support - + Symbolic inner products (using e.g. `Sympy`_) - - - Arbitrary spatial mode basis of functions - - Automatic on-the-fly inner product calculation (numeric or analytic if possible) - - Symbolic PDE equation specification - + + Numerical basis of function + Visualisation tools, e.g. based on the `movie-script`_ package * Long-term development track @@ -142,11 +135,12 @@ Forthcoming developments + Active advection + True quasi-geostrophic ocean when using ocean model version + Salinity in the ocean + + Symbolic PDE equation specification Contributing to qgs ------------------- -If you want to contribute actively to the roadmap detailed above, please contact the authors. +If you want to contribute actively to the roadmap detailed above, please contact the main authors. In addition, if you have made changes that you think will be useful to others, please feel free to suggest these as a pull request on the `qgs Github repository `_. @@ -154,7 +148,7 @@ A review of your pull request will follow with possibly suggestions of changes b Please consider the following guidelines before submitting: * Before submitting a pull request, double check that the branch to be merged contains only changes you wish to add to the master branch. This will save time in reviewing the code. -* For any changes to the core model files, please check your submission by :ref:`files/general_information:Running the tests` found in the folder `model_test <../../../../model_test>`_ to ensure that the model tensors are still valid. Please do not make changes to existing test cases. +* For any changes to the core model files, please check your submission by running the tests found in the folder `model_test <../../../../model_test>`_ to ensure that the model tensors are still valid (see the section :ref:`files/user_guide:5. Developers information` of the :ref:`files/user_guide:User guide`). Please do not make changes to existing test cases. * For substantial additions of code, including a test case (using `unittest`_) in the folder `model_test <../../../../model_test>`_ is recommended. * Please document the new functionalities in the documentation. Code addition without documentation addition will not be accepted. The documentation is done with `sphinx`_ and follows the Numpy conventions. Please take a look to the actual code to get an idea about how to document the code. * If your addition can be considered as a tool not directly related to the core of the model, please develop it in the toolbox folder. @@ -162,23 +156,6 @@ Please consider the following guidelines before submitting: For more information about git, Github and the pull request framework, a good source of information is the `contributing guide `_ of the `MITgcm `_. -Running the tests ------------------ - -.. TODO: move this to the user guide later. - -The model core tensors can be tested by running `pytest`_: :: - - pytest - -This will run all the tests and return a report. The test cases are written using `unittest`_. Additionally, test cases can be executed separately by running: :: - - python -m unittest model_test/test_name.py - -E.g., testing the MAOOAM inner products can be done by running: :: - - python -m unittest model_test/test_inner_products.py - Reporting issues with the software and getting support ------------------------------------------------------ @@ -217,5 +194,4 @@ References .. _beta-plane: https://en.wikipedia.org/wiki/Beta_plane .. _sparse: https://sparse.pydata.org/ .. _sphinx: https://www.sphinx-doc.org/en/master/ -.. _pytest: https://docs.pytest.org/en/stable/ .. _unittest: https://docs.python.org/3/library/unittest.html diff --git a/documentation/source/files/model/atmosphere.rst b/documentation/source/files/model/atmosphere.rst index 7d2ae6f..4604002 100644 --- a/documentation/source/files/model/atmosphere.rst +++ b/documentation/source/files/model/atmosphere.rst @@ -27,11 +27,11 @@ the vertical velocity :math:`\omega = \text{d}p/\text{d}t`, read & = +k'_d \nabla^2 (\psi^1_{\rm a}-\psi^3_{\rm a}) - \quad \ \frac{f_0}{\Delta p} \omega \nonumber \\ where :math:`\nabla^2` is the horizontal Laplacian. -The Coriolis parameter :math:`f` is linearized around a value :math:`f_0` (:attr:`~params.params.ScaleParams.f0`) evaluated at -latitude :math:`\phi_0` (:attr:`~params.params.ScaleParams.phi0_npi`), :math:`f = f_0 + \beta y`, with -:math:`\beta=\text{d}f/\text{d}y` (:attr:`~params.params.ScaleParams.beta`). The parameter :math:`k'_d` -(:attr:`~params.params.AtmosphericParams.kdp`) quantify the friction between the two atmospheric layers, -and :math:`\Delta p = 500` hPa (:attr:`~params.params.ScaleParams.deltap`) is the pressure difference between the atmospheric layers. +The Coriolis parameter :math:`f` is linearized around a value :math:`f_0` (:attr:`~.params.ScaleParams.f0`) evaluated at +latitude :math:`\phi_0` (:attr:`~.params.ScaleParams.phi0_npi`), :math:`f = f_0 + \beta y`, with +:math:`\beta=\text{d}f/\text{d}y` (:attr:`~.params.ScaleParams.beta`). The parameter :math:`k'_d` +(:attr:`~.params.AtmosphericParams.kdp`) quantify the friction between the two atmospheric layers, +and :math:`\Delta p = 500` hPa (:attr:`~.params.ScaleParams.deltap`) is the pressure difference between the atmospheric layers. Finally, :math:`J` is the Jacobian :math:`J(S, G) = \partial_x S\, \partial_y G - \partial_y S\, \partial_x G`. diff --git a/documentation/source/files/model/li_model.rst b/documentation/source/files/model/li_model.rst index bddaeb4..a6545df 100644 --- a/documentation/source/files/model/li_model.rst +++ b/documentation/source/files/model/li_model.rst @@ -3,7 +3,7 @@ Model with an orography and heat exchanges ========================================== In the :ref:`files/model/oro_model:Model with an orography and a temperature profile`, the radiative equilibrium temperature field at the middle of the atmosphere -is specified by a given profile :math:`\theta^\star` (:attr:`~params.params.AtmosphericTemperatureParams.thetas`) and the system relaxes to +is specified by a given profile :math:`\theta^\star` (:attr:`~.params.AtmosphericTemperatureParams.thetas`) and the system relaxes to this profile due to the `Newtonian cooling`_. In Li et al. :cite:`li-LHHBD2018`, another scheme for the temperature is proposed, based on the @@ -51,8 +51,8 @@ All the modes of this model version are expanded on the set of Fourier modes :ma and as in MAOOAM, the fields, parameters and variables are non-dimensionalized -by dividing time by :math:`f_0^{-1}` (:attr:`~params.params.ScaleParams.f0`), distance by -the characteristic length scale :math:`L` (:attr:`~params.params.ScaleParams.L`), pressure by the difference :math:`\Delta p` (:attr:`~params.params.ScaleParams.deltap`), +by dividing time by :math:`f_0^{-1}` (:attr:`~.params.ScaleParams.f0`), distance by +the characteristic length scale :math:`L` (:attr:`~.params.ScaleParams.L`), pressure by the difference :math:`\Delta p` (:attr:`~.params.ScaleParams.deltap`), temperature by :math:`f_0^2 L^2/R`, and streamfunction by :math:`L^2 f_0`. As a result of this non-dimensionalization, the fields :math:`\theta_{\rm a}` and :math:`\delta T_{\rm a}` can be identified: :math:`2 \theta_{\rm a} \equiv \delta T_{\rm a}`. @@ -68,15 +68,15 @@ The equations of the system of ordinary differential equations for this model th \dot\delta T_{{\rm g},i} & = & - \left(\lambda'_{\rm g}+ s_{B,{\rm g}}\right) \, \delta T_{{\rm g},i} + \left(2 \,\lambda'_{\rm g} + s_{B,{\rm a}}\right) \, \theta_{{\rm a},i} + C'_{{\rm g},i} where the parameters values have been replaced by their non-dimensional ones and we have also defined -:math:`G = - L^2/L_R^2` (:attr:`~params.params.QgParams.G`), -:math:`\lambda'_{{\rm a}} = \lambda/(\gamma_{\rm a} f_0)` (:attr:`~params.params.QgParams.Lpa`), -:math:`\lambda'_{{\rm g}} = \lambda/(\gamma_{\rm g} f_0)` (:attr:`~params.params.QgParams.Lpgo`), -:math:`S_{B,{\rm a}} = 8\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm a} f_0)` (:attr:`~params.params.QgParams.LSBpa`), -:math:`S_{B,{\rm g}} = 2\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm a} f_0)` (:attr:`~params.params.QgParams.LSBpgo`), -:math:`s_{B,{\rm a}} = 8\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm g} f_0)` (:attr:`~params.params.QgParams.sbpa`), -:math:`s_{B,{\rm g}} = 4\,\sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm g} f_0)` (:attr:`~params.params.QgParams.sbpgo`), -:math:`C'_{{\rm a},i} = R C_{{\rm a},i} / (2 \gamma_{\rm a} L^2 f_0^3)` (:attr:`~params.params.QgParams.Cpa`), -:math:`C'_{{\rm g},i} = R C_{{\rm g},i} / (\gamma_{\rm g} L^2 f_0^3)` (:attr:`~params.params.QgParams.Cpgo`). +:math:`G = - L^2/L_R^2` (:attr:`~.params.QgParams.G`), +:math:`\lambda'_{{\rm a}} = \lambda/(\gamma_{\rm a} f_0)` (:attr:`~.params.QgParams.Lpa`), +:math:`\lambda'_{{\rm g}} = \lambda/(\gamma_{\rm g} f_0)` (:attr:`~.params.QgParams.Lpgo`), +:math:`S_{B,{\rm a}} = 8\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm a} f_0)` (:attr:`~.params.QgParams.LSBpa`), +:math:`S_{B,{\rm g}} = 2\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm a} f_0)` (:attr:`~.params.QgParams.LSBpgo`), +:math:`s_{B,{\rm a}} = 8\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm g} f_0)` (:attr:`~.params.QgParams.sbpa`), +:math:`s_{B,{\rm g}} = 4\,\sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm g} f_0)` (:attr:`~.params.QgParams.sbpgo`), +:math:`C'_{{\rm a},i} = R C_{{\rm a},i} / (2 \gamma_{\rm a} L^2 f_0^3)` (:attr:`~.params.QgParams.Cpa`), +:math:`C'_{{\rm g},i} = R C_{{\rm g},i} / (\gamma_{\rm g} L^2 f_0^3)` (:attr:`~.params.QgParams.Cpgo`). The coefficients :math:`a_{i,j}`, :math:`g_{i, j, m}`, :math:`b_{i, j, m}` and :math:`c_{i, j}` are the inner products of the Fourier modes :math:`F_i`: diff --git a/documentation/source/files/model/maooam_model.rst b/documentation/source/files/model/maooam_model.rst index 81a9269..b17f408 100644 --- a/documentation/source/files/model/maooam_model.rst +++ b/documentation/source/files/model/maooam_model.rst @@ -25,9 +25,9 @@ The evolution equations for the atmospheric and oceanic streamfunctions defined \frac{\partial}{\partial t} \left( \nabla^2 \psi_\text{o} - \frac{\psi_\text{o}}{L_\text{R}^2} \right) + J(\psi_\text{o}, \nabla^2 \psi_\text{o}) + \beta \frac{\partial \psi_\text{o}}{\partial x} & = -r \nabla^2 \psi_\text{o} +\frac{C}{\rho_{\rm o} h} \nabla^2 (\psi^3_\text{a}-\psi_\text{o}).\nonumber -where :math:`\rho_{\rm o}` is the density of the ocean's water and :math:`h` is the depth of its layer (:attr:`~params.params.OceanicParams.h`). +where :math:`\rho_{\rm o}` is the density of the ocean's water and :math:`h` is the depth of its layer (:attr:`~.params.OceanicParams.h`). The rightmost term of the last equation represents the impact of the wind stress on the ocean, and is modulated -by the coefficient of the mechanical ocean-atmosphere coupling, :math:`d = C/(\rho_{\rm o} h)` (:attr:`~params.params.OceanicParams.d`). +by the coefficient of the mechanical ocean-atmosphere coupling, :math:`d = C/(\rho_{\rm o} h)` (:attr:`~.params.OceanicParams.d`). As we did for the :ref:`files/model/oro_model:Model with an orography and a temperature profile`, we rewrite these equations in terms of the `barotropic`_ streamfunction :math:`\psi_{\rm a}` and `baroclinic`_ streamfunction :math:`\theta_{\rm a}`: @@ -57,7 +57,7 @@ where :math:`\gamma_\text{a}` (:attr:`.AtmosphericTemperatureParams.gamma`) and (:attr:`.OceanicTemperatureParams.gamma`) are the heat capacities of the atmosphere and the active ocean layer. :math:`\lambda` (:attr:`~.AtmosphericTemperatureParams.hlambda`) is the heat transfer coefficient at the ocean-atmosphere interface. -:math:`\sigma` (:attr:`~params.params.AtmosphericParams.sigma`) is the static stability of the atmosphere, taken to be constant. +:math:`\sigma` (:attr:`~.params.AtmosphericParams.sigma`) is the static stability of the atmosphere, taken to be constant. The quartic terms represent the long-wave radiation fluxes between the ocean, the atmosphere, and outer space, with :math:`\epsilon_\text{a}` (:attr:`~.AtmosphericTemperatureParams.eps`) the emissivity of the grey-body atmosphere and @@ -120,7 +120,7 @@ Again, :math:`x` and :math:`y` are the horizontal adimensionalized coordinates d To easily manipulate these functions and the coefficients of the fields expansion, we number the basis functions along increasing values of :math:`H_{\rm o}` and then :math:`P_{\rm o}`. It allows to write the set as :math:`\left\{ \phi_i(x,y); 1 \leq i \leq n_\text{o}\right\}` where :math:`n_{\mathrm{o}}` -(:attr:`~params.params.QgParams.nmod` [1]) is the number of modes of the spectral expansion in the ocean. +(:attr:`~.params.QgParams.nmod` [1]) is the number of modes of the spectral expansion in the ocean. For example, the model derived in :cite:`mao-VDDG2015` can be specified by setting :math:`H_{\rm o} \in \{1,4\}`; :math:`P_{\rm o} \in \{1,2\}` and the set of basis functions is @@ -205,8 +205,8 @@ Ordinary differential equations ------------------------------- The fields, parameters and variables are non-dimensionalized -by dividing time by :math:`f_0^{-1}` (:attr:`~params.params.ScaleParams.f0`), distance by -the characteristic length scale :math:`L` (:attr:`~params.params.ScaleParams.L`), pressure by the difference :math:`\Delta p` (:attr:`~params.params.ScaleParams.deltap`), +by dividing time by :math:`f_0^{-1}` (:attr:`~.params.ScaleParams.f0`), distance by +the characteristic length scale :math:`L` (:attr:`~.params.ScaleParams.L`), pressure by the difference :math:`\Delta p` (:attr:`~.params.ScaleParams.deltap`), temperature by :math:`f_0^2 L^2/R`, and streamfunction by :math:`L^2 f_0`. As a result of this non-dimensionalization, the fields :math:`\theta_{\rm a}` and :math:`\delta T_{\rm a}` can be identified: :math:`2 \theta_{\rm a} \equiv \delta T_{\rm a}`. @@ -225,15 +225,15 @@ The ordinary differential equations of the truncated model are: \dot\delta T_{{\rm o},i} & = & - \sum_{j,m = 1}^{n_{\mathrm{o}}} \, O_{i,j,k} \, \psi_{{\rm o},j} \, \delta T_{{\rm o},k} - \left(\lambda'_{\rm o}+ s_{B,{\rm o}}\right) \, \delta T_{{\rm o},i} + \left(2 \,\lambda'_{\rm o} + s_{B,{\rm a}}\right) \, \sum_{j=1}^{n_{\mathrm{a}}} \, W_{i,j} \, \theta_{{\rm a},j} + \sum_{j=1}^{n_{\mathrm{a}}} \, W_{i,j} \, C'_{{\rm o},j} where the parameters values have been replaced by their non-dimensional ones and we have also defined -:math:`G = - L^2/L_R^2` (:attr:`~params.params.QgParams.G`), -:math:`\lambda'_{{\rm a}} = \lambda/(\gamma_{\rm a} f_0)` (:attr:`~params.params.QgParams.Lpa`), -:math:`\lambda'_{{\rm o}} = \lambda/(\gamma_{\rm o} f_0)` (:attr:`~params.params.QgParams.Lpgo`), -:math:`S_{B,{\rm a}} = 8\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm a} f_0)` (:attr:`~params.params.QgParams.LSBpa`), -:math:`S_{B,{\rm o}} = 2\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm a} f_0)` (:attr:`~params.params.QgParams.LSBpgo`), -:math:`s_{B,{\rm a}} = 8\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm o} f_0)` (:attr:`~params.params.QgParams.sbpa`), -:math:`s_{B,{\rm o}} = 4\,\sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm o} f_0)` (:attr:`~params.params.QgParams.sbpgo`), -:math:`C'_{{\rm a},i} = R C_{{\rm a},i} / (2 \gamma_{\rm a} L^2 f_0^3)` (:attr:`~params.params.QgParams.Cpa`), -:math:`C'_{{\rm o},i} = R C_{{\rm o},i} / (\gamma_{\rm o} L^2 f_0^3)` (:attr:`~params.params.QgParams.Cpgo`). +:math:`G = - L^2/L_R^2` (:attr:`~.params.QgParams.G`), +:math:`\lambda'_{{\rm a}} = \lambda/(\gamma_{\rm a} f_0)` (:attr:`~.params.QgParams.Lpa`), +:math:`\lambda'_{{\rm o}} = \lambda/(\gamma_{\rm o} f_0)` (:attr:`~.params.QgParams.Lpgo`), +:math:`S_{B,{\rm a}} = 8\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm a} f_0)` (:attr:`~.params.QgParams.LSBpa`), +:math:`S_{B,{\rm o}} = 2\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm a} f_0)` (:attr:`~.params.QgParams.LSBpgo`), +:math:`s_{B,{\rm a}} = 8\,\epsilon_{\rm a}\, \sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm o} f_0)` (:attr:`~.params.QgParams.sbpa`), +:math:`s_{B,{\rm o}} = 4\,\sigma_B \, T_{{\rm a},0}^3 / (\gamma_{\rm o} f_0)` (:attr:`~.params.QgParams.sbpgo`), +:math:`C'_{{\rm a},i} = R C_{{\rm a},i} / (2 \gamma_{\rm a} L^2 f_0^3)` (:attr:`~.params.QgParams.Cpa`), +:math:`C'_{{\rm o},i} = R C_{{\rm o},i} / (\gamma_{\rm o} L^2 f_0^3)` (:attr:`~.params.QgParams.Cpgo`). The coefficients :math:`a_{i,j}`, :math:`g_{i, j, m}`, :math:`b_{i, j, m}` and :math:`c_{i, j}` are the inner products of the Fourier modes :math:`F_i`: diff --git a/documentation/source/files/model/ocean.rst b/documentation/source/files/model/ocean.rst index e404ff0..b05db14 100644 --- a/documentation/source/files/model/ocean.rst +++ b/documentation/source/files/model/ocean.rst @@ -2,7 +2,7 @@ Oceanic component ================= -The oceanic component is a shallow-water active oceanic layer superimposed on a deep ocean layer at rest. +The oceanic component is a `shallow-water`_ active oceanic layer superimposed on a deep ocean layer at rest. The dynamics is given by the reduced-gravity quasi-geostrophic vorticity equation. Therefore, the equation of motion for the streamfunction :math:`\psi_\text{o}` of the ocean @@ -14,10 +14,10 @@ layer reads :cite:`oc-P2011` :cite:`oc-DDV2016` \frac{\psi_\text{o}}{L_\text{R}^2} \right) + J(\psi_\text{o}, \nabla^2 \psi_\text{o}) + \beta \frac{\partial \psi_\text{o}}{\partial x} = -r \nabla^2 \psi_\text{o}. -:math:`L_\text{R} = \sqrt{g' \, h }/ f_0` (:attr:`~params.params.QgParams.LR`) is the `reduced Rossby deformation radius`_ -where :math:`g'` (:attr:`~params.params.OceanParams.gp`) is the reduced gravity, :math:`h` is the depth of the layer (:attr:`~params.params.OceanicParams.h`), -and :math:`f_0` is the Coriolis parameter (:attr:`~params.params.ScaleParams.f0`). -:math:`r` (:attr:`~params.params.OceanicParams.r`) is the friction at the bottom of the active ocean layer. +:math:`L_\text{R} = \sqrt{g' \, h }/ f_0` (:attr:`~.params.QgParams.LR`) is the `reduced Rossby deformation radius`_ +where :math:`g'` (:attr:`~.params.OceanParams.gp`) is the reduced gravity, :math:`h` is the depth of the layer (:attr:`~.params.OceanicParams.h`), +and :math:`f_0` is the Coriolis parameter (:attr:`~.params.ScaleParams.f0`). +:math:`r` (:attr:`~.params.OceanicParams.r`) is the friction at the bottom of the active ocean layer. References ---------- @@ -28,3 +28,4 @@ References .. _MAOOAM: https://github.com/Climdyn/MAOOAM .. _reduced Rossby deformation radius: https://en.wikipedia.org/wiki/Rossby_radius_of_deformation +.. _shallow-water: https://en.wikipedia.org/wiki/Shallow_water_equations diff --git a/documentation/source/files/model/oro_model.rst b/documentation/source/files/model/oro_model.rst index c4a416d..46e391a 100644 --- a/documentation/source/files/model/oro_model.rst +++ b/documentation/source/files/model/oro_model.rst @@ -27,7 +27,7 @@ It transforms the equations in :ref:`files/model/atmosphere:Atmospheric componen \frac{\partial}{\partial t} \left( \nabla^2 \psi^3_{\rm a} \right) + J(\psi^3_{\rm a}, \nabla^2 \psi^3_{\rm a}) + J(\psi^3_{\rm a}, f_0 \, h/H_{\rm a}) + \beta \frac{\partial \psi^3_{\rm a}}{\partial x} & = +k'_d \nabla^2 (\psi^1_{\rm a}-\psi^3_{\rm a}) - k_d \nabla^2 \psi^3_{\rm a} - \frac{f_0}{\Delta p} \omega \nonumber \\ -where :math:`H_{\rm a}` is the characteristic depth of each fluid layer and the parameter :math:`k_d` (:attr:`~params.params.AtmosphericParams.kd`) quantifies the friction +where :math:`H_{\rm a}` is the characteristic depth of each fluid layer and the parameter :math:`k_d` (:attr:`~.params.AtmosphericParams.kd`) quantifies the friction between the atmosphere and the surface. :math:`h` is the topographic height field. Mid-layer equations and the thermal wind relation @@ -48,8 +48,8 @@ On the other hand, the thermodynamic equation governing the mean `potential temp 2 \, \frac{\partial}{\partial t} T_{\rm a} + J(\psi_{\rm a}, 2 T_{\rm a}) = - \frac{\sigma}{H_{\rm a}} \,\omega + 2 h_d \, (T^\ast - T_{\rm a}) -where :math:`h_d` (:attr:`~params.params.AtmosphericTemperatureParams.hd`) is the `Newtonian cooling`_ coefficient indicating the tendency to return to the equilibrium temperature profile :math:`T^\ast`. -:math:`\sigma` (:attr:`~params.params.AtmosphericParams.sigma`) is the static stability of the atmosphere, taken to be constant. +where :math:`h_d` (:attr:`~.params.AtmosphericTemperatureParams.hd`) is the `Newtonian cooling`_ coefficient indicating the tendency to return to the equilibrium temperature profile :math:`T^\ast`. +:math:`\sigma` (:attr:`~.params.AtmosphericParams.sigma`) is the static stability of the atmosphere, taken to be constant. The thermal wind relation .. math:: @@ -82,14 +82,14 @@ The fields are projected on Fourier modes respecting these boundary conditions: with integer values of :math:`M`, :math:`H`, :math:`P`. :math:`x` and :math:`y` are the horizontal adimensionalized coordinates, rescaled -by dividing the dimensional coordinates by the characteristic length :math:`L` (:attr:`~params.params.ScaleParams.L`). -The model's domain is then defined by :math:`(0 \leq x \leq \frac{2\pi}{n}, 0 \leq y \leq \pi)`, with :math:`n` (:attr:`~params.params.ScaleParams.n`) the aspect ratio -between its meridional and zonal extents :math:`L_y` (:attr:`~params.params.ScaleParams.L_y`) and :math:`L_x` (:attr:`~params.params.ScaleParams.L_x`). +by dividing the dimensional coordinates by the characteristic length :math:`L` (:attr:`~.params.ScaleParams.L`). +The model's domain is then defined by :math:`(0 \leq x \leq \frac{2\pi}{n}, 0 \leq y \leq \pi)`, with :math:`n` (:attr:`~.params.ScaleParams.n`) the aspect ratio +between its meridional and zonal extents :math:`L_y` (:attr:`~.params.ScaleParams.L_y`) and :math:`L_x` (:attr:`~.params.ScaleParams.L_x`). To easily manipulate these functions and the coefficients of the fields expansion, we number the basis functions along increasing values of :math:`M= H` and then :math:`P`. It allows to write the set as :math:`\left\{ F_i(x,y); 1 \leq i \leq n_\text{a}\right\}` where :math:`n_{\mathrm{a}}` -(:attr:`~params.params.QgParams.nmod` [0]) is the number of modes of the spectral expansion. +(:attr:`~.params.QgParams.nmod` [0]) is the number of modes of the spectral expansion. For example, with :math:`M=H=1` and :math:`P \in \{1,2\}`, one obtains the spectral truncation used by :cite:`om-CS1980`. The model derived in :cite:`om-RP1982` extended this set by two blocks of two functions each, and the @@ -142,16 +142,16 @@ the vertical velocity :math:`\omega(x,y)` also have to be decomposed into the ei h(x,y) & = & \sum_{i=1}^{n_{\mathrm{a}}} \, h_i \, F_i(x,y) \\ \omega(x,y) & = & \sum_{i=1}^{n_{\mathrm{a}}} \, \omega_i \, F_i(x,y) . -These fields can be specified in the model by setting the (non-dimensional) vectors :attr:`~params.params.GroundParams.hk` -and :attr:`~params.params.AtmosphericTemperatureParams.thetas`. Note that :math:`h` is scaled by the characteristic height :math:`H_{\rm a}` of each layer, +These fields can be specified in the model by setting the (non-dimensional) vectors :attr:`~.params.GroundParams.hk` +and :attr:`~.params.AtmosphericTemperatureParams.thetas`. Note that :math:`h` is scaled by the characteristic height :math:`H_{\rm a}` of each layer, and :math:`\theta^\star` is scaled by :math:`A f_0^2 L^2` (see section below). Ordinary differential equations ------------------------------- The fields, parameters and variables are non-dimensionalized -by dividing time by :math:`f_0^{-1}` (:attr:`~params.params.ScaleParams.f0`), distance by -the characteristic length scale :math:`L` (:attr:`~params.params.ScaleParams.L`), pressure by the difference :math:`\Delta p` (:attr:`~params.params.ScaleParams.deltap`), +by dividing time by :math:`f_0^{-1}` (:attr:`~.params.ScaleParams.f0`), distance by +the characteristic length scale :math:`L` (:attr:`~.params.ScaleParams.L`), pressure by the difference :math:`\Delta p` (:attr:`~.params.ScaleParams.deltap`), temperature by :math:`A f_0^2 L^2`, and streamfunction by :math:`L^2 f_0`. As stated above, a result of this non-dimensionalization is that the field :math:`\theta_{\rm a}` is identified with :math:`T_{\rm a}`: :math:`\theta_{\rm a} \equiv T_{\rm a}`. diff --git a/documentation/source/files/model_description.rst b/documentation/source/files/model_description.rst index eeff166..b18a592 100644 --- a/documentation/source/files/model_description.rst +++ b/documentation/source/files/model_description.rst @@ -7,7 +7,7 @@ models the dynamics of a 2-layer `quasi-geostrophic`_ (QG) channel atmosphere on a `beta-plane`_, coupled to a simple land or `shallow-water`_ ocean component. -It currently has three modes of operation: +It currently has three main modes of operation: * A QG atmosphere coupled to a land surface with topography through friction and with a simple thermal relaxation toward a climatological temperature. See :ref:`files/model/oro_model:Model with an orography and a temperature profile` for a description of this model version. @@ -42,6 +42,9 @@ It currently has three modes of operation: `doi:10.1007/s13351-018-8012-y `_ The shallow-water ocean can in principle be used as a stand-alone model, but this is not implemented for the moment. + +The modes over which the equations are projected can be modified, leading to other model version. +See :ref:`files/examples/VSPD:Recovering the result of Vannitsem, Solé-Pomies and De Cruz (2019)` for an example. More developments are yet to come, see the :ref:`files/general_information:Forthcoming developments` section. Components description diff --git a/documentation/source/files/technical/configuration.rst b/documentation/source/files/technical/configuration.rst index 1990a0a..d7022e1 100644 --- a/documentation/source/files/technical/configuration.rst +++ b/documentation/source/files/technical/configuration.rst @@ -2,9 +2,7 @@ Model configuration =================== -We describe here the modules used to configure the model. Please read the User Guide to learn how to use them. - -.. TODO: add more details and a link to it when the User Guide is ready. +We describe here the modules used to configure the model. Please read the :ref:`files/user_guide:User guide` to learn how to use them. .. automodule:: params.params :show-inheritance: diff --git a/documentation/source/files/technical_description.rst b/documentation/source/files/technical_description.rst index 50c879c..55fceef 100644 --- a/documentation/source/files/technical_description.rst +++ b/documentation/source/files/technical_description.rst @@ -3,7 +3,7 @@ Code Description ================ Presently, the `ordinary differential equations`_ (ODEs) of the qgs model are at most bilinear -in their variables :math:`\eta_i` (:math:`1\leq i\leq` :attr:`~params.params.QgParams.ndim`). +in their variables :math:`\eta_i` (:math:`1\leq i\leq` :attr:`~.params.QgParams.ndim`). This system of ODEs can therefore be expressed as the sum of a constant, a matrix multiplication, and a tensor contraction: @@ -58,7 +58,7 @@ Computational flow The computational flow is as follows: -1. The parameters are specified by instantiating a :class:`~params.params.QgParams` . +1. The parameters are specified by instantiating a :class:`~.params.QgParams` . 2. The inner products are computed and stored in :class:`~inner_products.analytic.AtmosphericInnerProducts` and :class:`~inner_products.analytic.OceanicInnerProducts` objects. 3. The tensor of the tendencies terms are computed in a :class:`~tensors.qgtensor.QgsTensor` object. 4. The functions :obj:`~functions.tendencies.create_tendencies` create Numba optimized functions that return the tendencies and the Jacobian matrix. diff --git a/documentation/source/files/user_guide.rst b/documentation/source/files/user_guide.rst index faa109a..70b410d 100644 --- a/documentation/source/files/user_guide.rst +++ b/documentation/source/files/user_guide.rst @@ -2,4 +2,517 @@ User guide ========== -Todo \ No newline at end of file +This guide explains how the qgs framework can be used and configured to obtain the various model version. + +1. Rationale behind qgs +------------------------ + +The purpose of qgs is to provide a Python callable representing the model's tendencies to the user, e.g. to provide +a function :math:`\boldsymbol{f}` (and also its Jacobian matrix) that can be integrated over time: + +.. math:: \dot{\boldsymbol{x}} = \boldsymbol{f}(t, \boldsymbol{x}) + +to obtain the model's trajectory. This callable can of course also be used to perform any other kind of model analysis. + +This guide provides: + +1. The different ways to configure qgs in order to obtain the function :math:`\boldsymbol{f}` are explained in the section :ref:`files/user_guide:2. Configuration of qgs`. +2. Some specific ways to configure qgs :ref:`files/user_guide:3. Using user-defined symbolic basis`. +3. Examples of usages of the model's tendencies function :math:`\boldsymbol{f}` are given in the section :ref:`files/user_guide:4. Using qgs (once configured)`. + +2. Configuration of qgs +----------------------- + +The computational flow to compute the function :math:`\boldsymbol{f}` of the model's tendencies +:math:`\dot{\boldsymbol{x}} = \boldsymbol{f}(t, \boldsymbol{x})` is explained in +:ref:`files/technical_description:Computational flow` and sketched below: + + +.. figure:: figures/compuflow.png + :scale: 70% + :align: center + + Sketch of the computational flow. + +This flow can be done step by step by the user, or can be automatically performed by the functions :meth:`~functions.tendencies.create_tendencies`. +This function takes a :class:`~.params.QgParams` parameters object and return the function :math:`\boldsymbol{f}` and its Jacobian matrix. +Optionally, it can also return the byproducts of the tendencies generation process: objects containing the inner products +between the model's spatial modes and the tensors of the model's tendencies terms. + +This section of the user guide explains how to configure the :class:`~.params.QgParams` object to obtain the desired model version and +model's parameters. + +2.1 Two different computational methods and three different geophysical components +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Two different methods are available to configure the models, related to two different methods of computing the inner products between +the functions of the spatial basis of qgs: + +1. **Symbolic:** It is the main way to specify the basis. The user can specify an arbitrary base of functions on which the model's partial differential equations will be projected. Depending on the dimension of the model version and on the computing power available, it might take a while to initialize. +2. **Analytic:** A legacy way of specifying the model version, in which the basis of functions are composed only of Fourier modes described in :ref:`files/model/maooam_model:Coupled ocean-atmosphere model (MAOOAM)` and is thus limited in term of modelling capacity. It is nevertheless the fastest way since it relies on analytic formulas derived in :cite:`user-DDV2016` to compute the basis functions. + +There are also three different geophysical components presently available in qgs: + +1. **Atmosphere:** This component is mandatory and provides a two-layer `quasi-geostrophic`_ atmosphere. See :ref:`files/model/atmosphere:Atmospheric component` for more details. +2. **Ocean:** This component provides a `shallow-water`_ active oceanic layer. See :ref:`files/model/ocean:Oceanic component` for more details. +3. **Ground:** This component provides a simple model for the ground (orography + heat exchange). It is activated by default if only the atmospheric component is defined, but then with only the orography (no heat exchange). + +The components needed by the user and their parameters have to be defined by instantiating a :class:`~.params.QgParams` object. +Creating this object, initializing these components and setting the parameters of the model is the subject of the next sections. + +2.2 Initializing qgs +^^^^^^^^^^^^^^^^^^^^ + +The model initialization first step requires the creation of a :class:`~.params.QgParams` object: + +.. code:: ipython3 + + from params.params import QgParams + model_parameters = QgParams() + +This object contains basically all the information needed by qgs to construct the inner products and the tendencies tensor of the model, +needed to finally produces the model's function :math:`\boldsymbol{f}`. + +The different components required by the user need then to be specified, by providing information about the basis of functions used to project +the partial differential equations of the model. As said before, two methods are available: + +2.2.1 The symbolic method +"""""""""""""""""""""""""" + +With this method, the user has to provide directly the basis of functions of each component with symbolic function expressions, using `Sympy`_. +This has to be done using a :class:`~basis.base.SymbolicBasis` object, which is basically a list of Sympy functions. + +The user can construct his own basis (see below) or use the various built-in Fourier basis provided with qgs: :class:`~basis.fourier.ChannelFourierBasis` or :class:`~basis.fourier.BasinFourierBasis`. +In the latter case, convenient constructor functions have been defined to help the user get the Fourier basis: :meth:`~basis.fourier.contiguous_basin_basis` and :meth:`~basis.fourier.contiguous_channel_basis`. +These functions create `contiguous` Fourier basis for two different kind of boundary conditions (a channel or a closed basin) shown on the first figure in :ref:`files/model/maooam_model:Coupled ocean-atmosphere model (MAOOAM)`. + +.. note:: + + A `contiguous` Fourier basis means here that the Fourier modes are all present in the model up to a given maximum wavenumber in each direction (`zonal and meridional`_). + Hence one has only to specify the maximum wavenumbers (and the model's domain aspect ratio) to these constructor functions. One can also create non-`contiguous` Fourier basis by specifying wavenumbers explicitely at + the :class:`~basis.fourier.ChannelFourierBasis` or :class:`~basis.fourier.BasinFourierBasis` instantiation. + +Once constructed, the basis has to be provided to the :class:`~.params.QgParams` object by using dedicated methods: :meth:`~.params.QgParams.set_atmospheric_modes`, :meth:`~.params.QgParams.set_oceanic_modes` and :meth:`~.params.QgParams.set_ground_modes`. +With the constructor functions, one can activate the mandatory atmospheric layer by typing + +.. code:: ipython3 + + from basis.fourier import contiguous_channel_basis + basis = contiguous_channel_basis(2, 2, 1.5) + model_parameters.set_atmospheric_modes(basis) + +where we have defined a channel Fourier basis up to wavenumber 2 in both directions and an aspect ratio of :math:`1.5`. + +.. note:: + + Please note that the aspect ratio of the basis object provided to qgs is not very important, because it is superseded by the aspect ratio sets in the :class:`~.params.QgParams` object. + +To activate the ocean or the ground components, the user has simply to use the method :meth:`~.params.QgParams.set_oceanic_modes` and :meth:`~.params.QgParams.set_ground_modes`. +Note that providing a oceanic basis of functions automatically deactivate the ground component, and vice-versa. + +Finally, since the `MAOOAM`_ Fourier basis are used frequently in qgs, convenient methods of the :class:`~.params.QgParams` object allow one to create easily these basis inside this object +(without the need to create them externally and then pass them to the qgs parameters object). These are the methods :meth:`~.params.QgParams.set_atmospheric_channel_fourier_modes`, :meth:`~.params.QgParams.set_oceanic_basin_fourier_modes` and :meth:`~.params.QgParams.set_ground_channel_fourier_modes`. +For instance, the effect obtained with the 3 previous lines of code (activating the atmosphere) can also be obtained by typing: + +.. code:: ipython3 + + model_parameters.set_atmospheric_channel_fourier_modes(2, 2, mode='symbolic') + +These convenient methods can also initialize qgs with another method (called `analytic`) and which is described in the next section. + +.. warning:: + + If you initialize one component with the symbolic method, all the other component that you define must be initialized with the same method. + +2.2.2 The analytic method +"""""""""""""""""""""""""" + +Computing the inner products of the symbolic functions defined with `Sympy`_ **can be very resources consuming**, therefore if the basis +of functions that you intend to use are the ones described in :ref:`files/model/maooam_model:Coupled ocean-atmosphere model (MAOOAM)`, you might be interested to use +the analytic method, which uses the analytic formula for the inner products given in :cite:`user-DDV2016`. This initializing mode can simply be used by using the +convenient methods of the :class:`~.params.QgParams` object: :meth:`~.params.QgParams.set_atmospheric_channel_fourier_modes`, :meth:`~.params.QgParams.set_oceanic_basin_fourier_modes` and :meth:`~.params.QgParams.set_ground_channel_fourier_modes`. + +For instance, to initialize a channel atmosphere with up to wavenumber 2 in both directions, one can simply write: + +.. code:: ipython3 + + model_parameters.set_atmospheric_channel_fourier_modes(2, 2, mode='analytic') + +Note that it is the default mode, so removing the `mode` argument will result in the same behavior. + +.. warning:: + + If you initialize one component with the analytic method, all the other component that you define must be initialized with the same method. + +2.3 Changing the default parameters of qgs +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +This section deals with how to change the parameters of qgs. As stated in the :ref:`files/technical/configuration:The model's parameters module` section of the :ref:`files/references:References`, +there are seven types of parameters arranged in classes: + +* :class:`~.params.ScaleParams` contains the model scale parameters. +* :class:`~.params.AtmosphericParams` contains the atmospheric dynamical parameters. +* :class:`~.params.AtmosphericTemperatureParams` containing the atmosphere's temperature and heat-exchange parameters. +* :class:`~.params.OceanicParams` contains the oceanic dynamical parameters. +* :class:`~.params.OceanicTemperatureParams` contains the ocean's temperature and heat-exchange parameters. +* :class:`~.params.GroundParams` contains the ground dynamical parameters (e.g. orography). +* :class:`~.params.GroundTemperatureParams` contains the ground's temperature and heat-exchange parameters. + +These parameters classes are regrouped into the global structure :class:`~.params.QgParams` and are accessible through the attributes: + +* :attr:`~.params.QgParams.scale_params` for :class:`~.params.ScaleParams`. +* :attr:`~.params.QgParams.atmospheric_params` for :class:`~.params.AtmosphericParams`. +* :attr:`~.params.QgParams.atemperature_params` for :class:`~.params.AtmosphericTemperatureParams`. +* :attr:`~.params.QgParams.oceanic_params` for :class:`~.params.OceanicParams`. +* :attr:`~.params.QgParams.otemperature_params` for :class:`~.params.OceanicTemperatureParams`. +* :attr:`~.params.QgParams.ground_params` for :class:`~.params.GroundParams`. +* :attr:`~.params.QgParams.otemperature_params` for :class:`~.params.GroundTemperatureParams`. + +The parameters inside these structures can be changed by passing a dictionnary of the new values to the :meth:`~.params.QgParams.set_params` method. For example, if one want to change the +Coriolis parameter :math:`f_0` and the static stability of the atmosphere :math:`\sigma`, one has to write: + +.. code:: ipython3 + + model_parameters.set_params({'f0': 1.195e-4, 'sigma':0.14916}) + +where :obj:`model_parameters` is an instance of the :class:`~.params.QgParams` class. This method will find where the parameters are stored and will peform the +subsitution. However, some parameters may not have an unique name, for instance there is a parameter :attr:`T0` for the stationary solution :math:`T_0` of the 0-th order temperature for both the +atmosphere and the ocean. In this case, one need to find to which parameter class the parameter belong, and then call the :meth:`set_params` of the corresponding object. +For example, changing the parameter :attr:`~.params.AtmosphericTemperatureParams.T0` in the atmosphere can be done with: + +.. code:: ipython3 + + model_parameters.atemperature_params.set_params({'T0': 280.}) + + +Finally, some specific methods allow to setup expansion [#expansion]_ coefficients. +Presently these are the :attr:`.AtmosphericTemperatureParams.set_thetas`, :attr:`.AtmosphericTemperatureParams.set_insolation`, +:attr:`.OceanicTemperatureParams.set_insolation`, :attr:`.GroundTemperatureParams.set_insolation` +and :attr:`.GroundParams.set_orography` methods. For example, to activate the Newtonian cooling, one has to write: + +.. code:: ipython3 + + model_parameters.atemperature_params.set_thetas(0.1, 0) + +which indicates that the first component [#component]_ of the radiative equilibrium mean temperature should be equal to :math:`0.1`. + +.. note:: + + Using both the atmospheric Newtonian cooling coefficients with :attr:`.AtmosphericTemperatureParams.set_thetas` and the heat exchange scheme :attr:`.AtmosphericTemperatureParams.set_insolation` + together doesn't make so much sense. Using the Newtonian cooling scheme is useful when one wants to use the atmospheric model alone, while using the heat exchange scheme is useful when the atmosphere is + connected to another component lying beneath it (ocean or ground). + +Similarly, one activate the orography by typing: + +.. code:: ipython3 + + model_parameters.ground_params.set_orography(0.2, 1) + +We refer the reader to the description of these methods for more details (just click on the link above to get there). + +Once your model is configured, you can review the list of parameters by calling the method :meth:`.QgParams.print_params`: + +.. code:: ipython3 + + model_parameters.print_params() + +2.4 Creating the tendencies function +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Once you have configured your :class:`.QgParams`, it is very simple to obtain the model's tendencies :math:`\boldsymbol{f}` and the its +Jacobian matrix :math:`\boldsymbol{\mathrm{Df}}` by using the function :func:`.create_tendencies` and running: + +.. code:: ipython3 + + from functions.tendencies import create_tendencies + + f, Df = create_tendencies(model_parameters) + +The function :meth:`f` hence produced can be used to generate the model's trajectories. +See the section :ref:`files/user_guide:4. Using qgs (once configured)` for the possible usage. + +2.5 Saving your model +^^^^^^^^^^^^^^^^^^^^^^^ + +The simplest way to save your model is to `pickle`_ the functions generating the model's tendencies and the Jacobian matrix. +Hence, using the same name as in the previous section, one can type: + +.. code:: ipython3 + + import pickle + + # saving the model + model={'f': f, 'Df': Df, 'parameters': model_parameters} + + with open('model.pickle', "wb") as ff: + pickle.dump(model, ff, pickle.HIGHEST_PROTOCOL) + +and it can be loaded again by typing + +.. code:: ipython3 + + from params.params import QgParams + + # loading the model + with open('model.pickle', "rb") as ff: + model = pickle.load(ff) + + f = model['f'] + model_parameters = model['parameters'] + +.. warning:: + + Due to several different possible reasons, loading models saved previously on another machine may not work. + The only thing to do is then to recompute the model tendencies with the loaded model parameters (using the function + :func:`.create_tendencies`. In this case, it is better to save only the model parameters: + + .. code:: ipython3 + + import pickle + + # saving the model + + with open('model_parameters.pickle', "wb") as ff: + pickle.dump(model_parameters, ff, pickle.HIGHEST_PROTOCOL) + +It is also possible to save the inner products and/or the tensor storing the terms of the model's tendencies. For instance, the function +:func:`.create_tendencies` allows to obtain these information: + +.. code:: ipython3 + + f, Df, inner_products, tensor = create_tendencies(model_parameters, return_inner_products=True, return_qgtensor=True) + +The objects :class:`.QgParams`, the inner products, and the object :class:`.QgsTensor` hence obtained can be saved using `pickle`_ or the built-in +:meth:`save_to_file` methods (respectively :meth:`.QgParams.save_to_file`, :meth:`.AtmosphericInnerProducts.save_to_file` and :meth:`.QgsTensor.save_to_file`). + +Using these object, it is possible to reconstruct by hand the model's tendencies (see the section :ref:`files/user_guide:3.2 A more involved example: Manually setting the basis and the inner products definition` for an example). + +3. Using user-defined symbolic basis +-------------------------------------- + +3.1 A simple example +^^^^^^^^^^^^^^^^^^^^^ + +In this simple example, we are going to create an atmosphere-ocean coupled model like in :cite:`user-VDDG2015`, but with some atmospheric modes missing. + +First we create the parameters object: + +.. code:: ipython3 + + from params.params import QgParams + model_parameters = QgParams({'n': 1.5}) + +and we create a :class:`.ChannelFourierBasis` with all the modes up to wavenumber 2 in both except the ones with wavenumbers 1 and 2 in respectively the :math:`x` and :math:`y` direction: + +.. code:: ipython3 + + from basis.fourier import ChannelFourierBasis + atm_basis = ChannelFourierBasis(np.array([[1,1], + [2,1], + [2,2]]),1.5) + model_parameters.set_atmospheric_modes(atm_basis) + +Finally we add the same ocean as in :cite:`user-VDDG2015`: + +.. code:: ipython3 + + model_parameters.set_oceanic_basin_fourier_modes(2,4,mode='symbolic') + +The last step is to set the parameters according to your needs (as seen in the section :ref:`files/user_guide:2.3 Changing the default parameters of qgs`). + +The model hence configured can be passed to the function creating the model's tendencies :math:`\boldsymbol{f}`, as detailed in the section :ref:`files/user_guide:2.4 Creating the tendencies function`. + +3.2 A more involved example: Manually setting the basis and the inner products definition +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. warning:: + + This initialization method is not yet well-defined in the model. It builds the model block by block to construct an ad-hoc model version. + +In this section, we describe how to setup an user-defined basis for one of the model's component. We will do it for the ocean, but +the approach is similar for the other components. We will project the ocean equations on four modes proposed in :cite:`user-P2011`: + +.. math:: + + \tilde\phi_1(x,y) & = & 2\, e^{-\alpha x} \, \sin(\frac{n}{2} x)\, \sin(y), \nonumber \\ + \tilde\phi_2(x,y) & = & 2\, e^{-\alpha x} \, \sin(n x)\, \sin(y), \nonumber \\ + \tilde\phi_3(x,y) & = & 2\, e^{-\alpha x} \, \sin(\frac{n}{2} x)\, \sin(2 y), \nonumber \\ + \tilde\phi_4(x,y) & = & 2\, e^{-\alpha x} \, \sin(n x)\, \sin(2 y), \nonumber \\ + +and connected it to the channel atmosphere defined in the sections above, using a symbolic base of functions. +First we create the parameters object and the atmosphere: + +.. code:: ipython3 + + from params.params import QgParams + model_parameters = QgParams({'n': 1.5}) + model_parameters.set_atmospheric_channel_fourier_modes(2, 2, mode="symbolic") + +Then we create a :class:`.SymbolicBasis` object: + +.. code:: ipython3 + + from basis.base import SymbolicBasis + ocean_basis = SymbolicBasis() + +We must then specify the function of the basis using `Sympy`_: + +.. code:: ipython3 + + from sympy import symbols, sin, cos, sqrt, exp + x, y = symbols('x y') # x and y coordinates on the model's spatial domain + n, al = symbols('n al') # aspect ratio and alpha coefficients + for i in range(1, 3): + for j in range(1, 3): + ocean_basis.functions.append(2 * exp(- al * x) * sin(j * n * x / 2) * sin(i * y)) + +We then set the value of the parameter :math:`\alpha` to a certain value (here :math:`\alpha=1`). Please note that the +:math:`\alpha` is then an extrinsic parameter of the model that you have to specify through a substitution: + +.. code:: ipython3 + + ocean_basis.substitutions.append(('al', 1.)) + +The base of function hence defined needs to be passed to the model's parameter object: + +.. code:: ipython3 + + model_parameters.set_oceanic_modes(ocean_basis) + +and the user can set the parameters according to it needs (as seen in the previous section). + +Additionally, for these particular basis function, special inner products needs to be defined instead of the standard one proposed in :ref:`files/model/maooam_model:Coupled ocean-atmosphere model (MAOOAM)`. +We consider thus as in :cite:`user-P2011` and :cite:`user-VD2014` the following inner product: + +.. math:: + + (S, G) = \frac{n}{2\pi^2}\int_0^\pi\int_0^{2\pi/n} e^{2 \alpha x} \, S(x,y)\, G(x,y)\, \mathrm{d} x \, \mathrm{d} y + +This can be specified in the model by creating a user-defined subclass of the :class:`.SymbolicInnerProductDefinition` class defining +the expression of the inner product: + +.. code:: ipython3 + + from sympy import pi + class UserInnerProductDefinition(StandardSymbolicInnerProductDefinition): + + def symbolic_inner_product(self, S, G): + """Function defining the inner product to be computed symbolically: + :math:`(S, G) = \\frac{n}{2\\pi^2}\\int_0^\\pi\\int_0^{2\\pi/n} e^{2 \\alpha x} \\, S(x,y)\\, G(x,y)\\, \\mathrm{d} x \\, \\mathrm{d} y`. + + Parameters + ---------- + S: Sympy expression + Left-hand side function of the product. + G: Sympy expression + Right-hand side function of the product. + + Returns + ------- + Sympy expression + The result of the symbolic integration + """ + expr = (n / (2 * pi ** 2)) * exp(2 * al * x) * S * G + return self.integrate_over_domain(self.optimizer(expr)) + + +and passing it to an :class:`.OceanicSymbolicInnerProducts` object: + +.. code:: ipython3 + + ip = UserInnerProductDefinition() + + from inner_products.symbolic import AtmosphericSymbolicInnerProducts, OceanicSymbolicInnerProducts + + aip = AtmosphericSymbolicInnerProducts(model_parameters) + oip = OceanicSymbolicInnerProducts(model_parameters, inner_product_definition=ip) + + +This will compute the inner products and may take a certain time (depending on your number of cpus available). +Once computed, the corresponding tenedencies must then be created manually, first by creating the :class:`.QgsTensor` object: + +.. code:: ipython3 + + from tensors.qgtensor import QgsTensor + + aotensor = QgsTensor(model_parameters, aip, oip) + +and then finally creating the Python-`Numba`_ callable for the model's tendencies :math:`\boldsymbol{f}`: + +.. code:: ipython3 + + from numba import njit + coo = aotensor.tensor.coords.T + val = aotensor.tensor.data + + @njit + def f(t, x): + xx = np.concatenate((np.full((1,), 1.), x)) + xr = sparse_mul3(coo, val, xx, xx) + + return xr[1:] + +This conclude this section, the function :meth:`f` hence produced can be used to generate the model's trajectories. +See the following section for the possible usage. + +4. Using qgs (once configured) +--------------------------------- + +Once the function :math:`\boldsymbol{f}` giving the model's tendencies has been obtained, it is possible to use it with +the qgs built-in integrator to obtain model's trajectories: + +.. code:: python3 + + from integrators.integrator import RungeKuttaIntegrator + import numpy as np + + integrator = RungeKuttaIntegrator() + integrator.set_func(f) + + ic = np.random.rand(model_parameters.ndim)*0.1 # create a vector of initial conditions with the same dimension as the model + integrator.integrate(0., 3000000., 0.1, ic=ic, write_steps=1) + time, trajectory = integrator.get_trajectories() + +Note that it is also possible to use other integrator available on the market, see for instance the :ref:`files/examples/diffeq_example:Example of DiffEqPy usage`. + +Other examples of usage will be added to the present section soon. + +5. Developers information +------------------------- + +5.1 Running the test +^^^^^^^^^^^^^^^^^^^^^ + +The model core tensors can be tested by running `pytest`_ in the main folder: :: + + pytest + +This will run all the tests and return a report. The test cases are written using `unittest`_. Additionally, test cases can be executed separately by running: :: + + python -m unittest model_test/test_name.py + +E.g., testing the MAOOAM inner products can be done by running: :: + + python -m unittest model_test/test_inner_products.py + +References +----------- + +.. bibliography:: model/ref.bib + :labelprefix: USER- + :keyprefix: user- + +.. rubric:: Footnotes + +.. [#expansion] Generally in term of the specified basis functions. +.. [#component] The component corresponding to the first basis function of the atmopshere. + +.. _quasi-geostrophic: https://en.wikipedia.org/wiki/Quasi-geostrophic_equations +.. _shallow-water: https://en.wikipedia.org/wiki/Shallow_water_equations +.. _Sympy: https://www.sympy.org/ +.. _zonal and meridional: https://en.wikipedia.org/wiki/Zonal_and_meridional_flow +.. _MAOOAM: https://github.com/Climdyn/MAOOAM +.. _pytest: https://docs.pytest.org/en/stable/ +.. _unittest: https://docs.python.org/3/library/unittest.html +.. _Numba: https://numba.pydata.org/ +.. _pickle: https://docs.python.org/3.8/library/pickle.html diff --git a/documentation/source/index.rst b/documentation/source/index.rst index 54adfa1..42a0914 100644 --- a/documentation/source/index.rst +++ b/documentation/source/index.rst @@ -35,6 +35,13 @@ About Part of the code comes from the Python `MAOOAM`_ implementation by Maxime Tondeur and Jonathan Demaeyer. qgs is licensed under the `MIT`_ license: +**Please cite the code description article if you use (a part of) this software for a publication:** + +* Demaeyer J., De Cruz, L. and Vannitsem, S. , (2020). qgs: A flexible Python framework of reduced-order multiscale climate models. + *Journal of Open Source Software*, **5**\(56), 2597, `https://doi.org/10.21105/joss.02597 `_. + +Please consult the qgs `code repository `_ for updates. + .. code-block:: none The MIT License (MIT) diff --git a/functions/sparse_mul.py b/functions/sparse_mul.py index c1fa5b0..d939bb8 100644 --- a/functions/sparse_mul.py +++ b/functions/sparse_mul.py @@ -29,14 +29,14 @@ def sparse_mul3(coo, value, vec_a, vec_b): value: ~numpy.ndarray(float) A 1D array of shape (n_elems,), a list of value in the tensor vec_a: ~numpy.ndarray(float) - The vector :math:`a_j` to contract the tensor with. Must be of shape (:attr:`~params.params.QgParams.ndim` + 1,). + The vector :math:`a_j` to contract the tensor with. Must be of shape (:attr:`~.params.QgParams.ndim` + 1,). vec_b: ~numpy.ndarray(float) - The vector :math:`b_k` to contract the tensor with. Must be of shape (:attr:`~params.params.QgParams.ndim` + 1,). + The vector :math:`b_k` to contract the tensor with. Must be of shape (:attr:`~.params.QgParams.ndim` + 1,). Returns ------- ~numpy.ndarray(float) - The vector :math:`v_i`, of shape (:attr:`~params.params.QgParams.ndim` + 1,). + The vector :math:`v_i`, of shape (:attr:`~.params.QgParams.ndim` + 1,). """ res = np.zeros_like(vec_a) n_elems = coo.shape[0] @@ -65,12 +65,12 @@ def sparse_mul2(coo, value, vec): value: ~numpy.ndarray(float) A 1D array of shape (n_elems,), a list of value in the tensor vec: ~numpy.ndarray(float) - The vector :math:`a_k` to contract the tensor with. Must be of shape (:attr:`~params.params.QgParams.ndim` + 1,). + The vector :math:`a_k` to contract the tensor with. Must be of shape (:attr:`~.params.QgParams.ndim` + 1,). Returns ------- ~numpy.ndarray(float) - The matrix :math:`A_{i,j}`, of shape (:attr:`~params.params.QgParams.ndim` + 1, :attr:`~params.params.QgParams.ndim` + 1). + The matrix :math:`A_{i,j}`, of shape (:attr:`~.params.QgParams.ndim` + 1, :attr:`~.params.QgParams.ndim` + 1). """ res = np.zeros((len(vec), len(vec))) diff --git a/functions/tendencies.py b/functions/tendencies.py index 9a081b3..e7941cd 100644 --- a/functions/tendencies.py +++ b/functions/tendencies.py @@ -33,7 +33,7 @@ def create_tendencies(params, return_inner_products=False, return_qgtensor=False Parameters ---------- - params: ~params.params.QgParams + params: ~.params.QgParams The parameters fully specifying the model configuration. return_inner_products: bool If True, return the inner products of the model. Default to False. diff --git a/inner_products/analytic.py b/inner_products/analytic.py index 044f993..dd706be 100644 --- a/inner_products/analytic.py +++ b/inner_products/analytic.py @@ -50,7 +50,7 @@ class AtmosphericAnalyticInnerProducts(AtmosphericInnerProducts): Parameters ---------- - params: None or ~params.params.QgParams or list, optional + params: None or ~.params.QgParams or list, optional An instance of model's parameters object or a list in the form [aspect_ratio, ablocks, natm]. If a list is provided, `aspect_ratio` is the aspect ratio of the domain, `ablocks` is a spectral blocks detailing the model's atmospheric modes :math:`x`-and :math:`y`-wavenumber as an array of shape (natm, 2), and `natm` is @@ -74,7 +74,7 @@ class AtmosphericAnalyticInnerProducts(AtmosphericInnerProducts): stored: bool Indicate if the inner products are stored in the object. atmospheric_wavenumbers: ~numpy.ndarray(WaveNumber) - An array of shape (:attr:`~params.params.QgParams.nmod` [0], ) of the wavenumber object of each mode. + An array of shape (:attr:`~.params.QgParams.nmod` [0], ) of the wavenumber object of each mode. """ def __init__(self, params=None, stored=True): @@ -435,7 +435,7 @@ class OceanicAnalyticInnerProducts(OceanicInnerProducts): Parameters ---------- - params: None or ~params.params.QgParams or list, optional + params: None or ~.params.QgParams or list, optional An instance of model's parameters object or a list in the form [aspect_ratio, oblocks, noc]. If a list is provided, `aspect_ratio` is the aspect ratio of the domain, `ablocks` is a spectral blocks detailing the model's oceanic modes :math:`x`-and :math:`y`-wavenumber as an array of shape (noc, 2), and `noc` is @@ -455,7 +455,7 @@ class OceanicAnalyticInnerProducts(OceanicInnerProducts): stored: bool Indicate if the inner products are stored in the object. oceanic_wavenumbers: ~numpy.ndarray(WaveNumber) - An array of shape (:attr:`~params.params.QgParams.nmod` [1], ) of the wavenumber object of each mode. + An array of shape (:attr:`~.params.QgParams.nmod` [1], ) of the wavenumber object of each mode. """ @@ -687,7 +687,7 @@ class GroundAnalyticInnerProducts(GroundInnerProducts): Parameters ---------- - params: None or ~params.params.QgParams or list, optional + params: None or ~.params.QgParams or list, optional An instance of model's parameters object or a list in the form [aspect_ratio, gblocks, ngr]. If a list is provided, `aspect_ratio` is the aspect ratio of the domain, `gblocks` is a spectral blocks detailing the model's oceanic modes :math:`x`-and :math:`y`-wavenumber as an array of shape (ngr, 2), and `ngr` is @@ -707,7 +707,7 @@ class GroundAnalyticInnerProducts(GroundInnerProducts): stored: bool Indicate if the inner products are stored in the object. ground_wavenumbers: ~numpy.ndarray(WaveNumber) - An array of shape (:attr:`~params.params.QgParams.nmod` [1], ) of the wavenumber object of each mode. + An array of shape (:attr:`~.params.QgParams.nmod` [1], ) of the wavenumber object of each mode. """ diff --git a/inner_products/base.py b/inner_products/base.py index a18b43a..9bda36d 100644 --- a/inner_products/base.py +++ b/inner_products/base.py @@ -8,7 +8,7 @@ Description of the classes -------------------------- - The two classes computing and holding the inner products of the basis functions are: + The three classes computing and holding the inner products of the basis functions are: * :class:`AtmosphericInnerProducts` * :class:`OceanicInnerProducts` diff --git a/inner_products/definition.py b/inner_products/definition.py index 7aa8ed6..6946247 100644 --- a/inner_products/definition.py +++ b/inner_products/definition.py @@ -15,7 +15,7 @@ from sympy.simplify.fu import TR8, TR10 from sympy import diff, integrate, symbols, pi -_n = symbols('n') +_n = symbols('n', real=True, nonnegative=True) _x, _y = symbols('x y') diff --git a/inner_products/symbolic.py b/inner_products/symbolic.py index 1100e3c..d9149fe 100644 --- a/inner_products/symbolic.py +++ b/inner_products/symbolic.py @@ -30,10 +30,11 @@ from inner_products.definition import StandardSymbolicInnerProductDefinition from sympy import symbols -_n = symbols('n') +_n = symbols('n', real=True, nonnegative=True) _x, _y = symbols('x y') -# TODO: Add warnings if trying to connect analytic and symbolic inner products together +# TODO: - Add warnings if trying to connect analytic and symbolic inner products together +# - Switch to numerical integration of the inner products if the symbolic one is to long (to be done when NumericBasis is ready) class AtmosphericSymbolicInnerProducts(AtmosphericInnerProducts): @@ -42,15 +43,24 @@ class AtmosphericSymbolicInnerProducts(AtmosphericInnerProducts): Parameters ---------- - params: ~params.params.QgParams or list + params: ~.params.QgParams or list An instance of model's parameters object or a list in the form [aspect_ratio, atmospheric_basis, basis, oog, oro_basis]. If a list is provided, `aspect_ratio` is the aspect ratio of the domain, `atmospheric_basis` is a SymbolicBasis with the modes of the atmosphere, and `ocean_basis` is either `None` or a SymbolicBasis object with the modes of the ocean or the ground. Finally `oog` indicates if it is an ocean or a ground component that is connected, by setting it to `ocean` or to 'ground', and in this latter case, `oro_basis` indicates on which basis the orography is developed. + stored: bool, optional + Indicate if the inner product must be stored or computed on the fly. Default to `True` inner_product_definition: None or InnerProductDefinition, optional The definition of the inner product being used. If `None`, use the canonical StandardInnerProductDefinition object. Default to `None`. + interaction_inner_product_definition: None or InnerProductDefinition, optional + The definition of the inner product being used for the interaction with the other components, i.e. to compute the inner products with the other component base of funcitons. + If `None`, use the `inner_product_definition` provided. + Default to `None`. + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. Attributes ---------- @@ -62,14 +72,18 @@ class AtmosphericSymbolicInnerProducts(AtmosphericInnerProducts): Object holding the symbolic modes of the ocean (or `None` if there is no ocean). connected_to_ocean: bool Indicate if the atmosphere is connected to an ocean. + stored: bool + Indicate if the inner product must be stored or computed on the fly. ip: InnerProductDefinition Object defining the inner product. + iip: InnerProductDefinition + Object defining the interaction inner product. subs: list(tuple) List of 2-tuples containing the substitutions to be made with the functions after the inner products symbolic computation. """ - def __init__(self, params=None, stored=True, inner_product_definition=None, num_threads=None): + def __init__(self, params=None, stored=True, inner_product_definition=None, interaction_inner_product_definition=None, num_threads=None): AtmosphericInnerProducts.__init__(self) @@ -117,6 +131,11 @@ def __init__(self, params=None, stored=True, inner_product_definition=None, num_ else: self.ip = inner_product_definition + if interaction_inner_product_definition is None: + self.iip = self.ip + else: + self.iip = interaction_inner_product_definition + self.stored = stored if stored: self.compute_inner_products(num_threads) @@ -142,6 +161,9 @@ def connect_to_ocean(self, ocean_basis, num_threads=None): ---------- ocean_basis: SymbolicBasis or OceanicSymbolicInnerProducts Basis of function of the ocean or a symbolic oceanic inner products object containing the basis. + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. """ if isinstance(ocean_basis, OceanicSymbolicInnerProducts): ocean_basis = ocean_basis.oceanic_basis @@ -163,7 +185,7 @@ def connect_to_ocean(self, ocean_basis, num_threads=None): self._s = sp.zeros((self.natm, noc), dtype=float, format='dok') # d inner products - args_list = [[(i, j), self.ip.ip_lap, (self._F(i), self._phi(j))] for i in range(self.natm) + args_list = [[(i, j), self.iip.ip_lap, (self._F(i), self._phi(j))] for i in range(self.natm) for j in range(noc)] result = pool.map(apply, args_list) @@ -173,7 +195,7 @@ def connect_to_ocean(self, ocean_basis, num_threads=None): .subs(self.oceanic_basis.substitutions)) # s inner products - args_list = [[(i, j), self.ip.symbolic_inner_product, (self._F(i), self._phi(j))] for i in range(self.natm) + args_list = [[(i, j), self.iip.symbolic_inner_product, (self._F(i), self._phi(j))] for i in range(self.natm) for j in range(noc)] result = pool.map(apply, args_list) @@ -197,6 +219,9 @@ def connect_to_ground(self, ground_basis, orographic_basis, num_threads=None): orographic_basis: str String to select which component basis modes to use to develop the orography in series. Can be either 'atmospheric' or 'ground'. + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. """ if isinstance(ground_basis, GroundSymbolicInnerProducts): @@ -223,7 +248,7 @@ def connect_to_ground(self, ground_basis, orographic_basis, num_threads=None): self._s = sp.zeros((self.natm, ngr), dtype=float, format='dok') # s inner products - args_list = [[(i, j), self.ip.symbolic_inner_product, (self._F(i), self._phi(j))] for i in range(self.natm) + args_list = [[(i, j), self.iip.symbolic_inner_product, (self._F(i), self._phi(j))] for i in range(self.natm) for j in range(ngr)] result = pool.map(apply, args_list) @@ -233,7 +258,7 @@ def connect_to_ground(self, ground_basis, orographic_basis, num_threads=None): .subs(self.ground_basis.substitutions)) # gh inner products - args_list = [[(i, j, k), self.ip.ip_jac, (self._F(i), self._F(j), self._phi(k))] for i in range(self.natm) + args_list = [[(i, j, k), self.iip.ip_jac, (self._F(i), self._F(j), self._phi(k))] for i in range(self.natm) for j in range(self.natm) for k in range(ngr)] result = pool.map(apply, args_list) @@ -249,6 +274,14 @@ def connect_to_ground(self, ground_basis, orographic_basis, num_threads=None): pool.terminate() def compute_inner_products(self, num_threads=None): + """Function computing and storing all the inner products at once. + + Parameters + ---------- + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. + """ self._a = sp.zeros((self.natm, self.natm), dtype=float, format='dok') self._u = sp.zeros((self.natm, self.natm), dtype=float, format='dok') @@ -373,7 +406,7 @@ def gh(self, i, j, k): if self.stored and self._gh is not None: return self._gh[i, j, k] else: - res = self.ip.ip_jac(self._F(i), self._F(j), self._phi(k)) + res = self.iip.ip_jac(self._F(i), self._F(j), self._phi(k)) return float(res.subs(self.subs).subs(self.atmospheric_basis.substitutions).subs(extra_subs)) else: @@ -392,7 +425,7 @@ def s(self, i, j): if self.stored and self._s is not None: return self._s[i, j] else: - res = self.ip.symbolic_inner_product(self._F(i), self._phi(j)) + res = self.iip.symbolic_inner_product(self._F(i), self._phi(j)) return float(res.subs(self.subs).subs(self.atmospheric_basis.substitutions).subs(extra_subs)) else: return 0 @@ -410,7 +443,7 @@ def d(self, i, j): if self.stored and self._d is not None: return self._d[i, j] else: - res = self.ip.ip_lap(self._F(i), self._phi(j)) + res = self.iip.ip_lap(self._F(i), self._phi(j)) return float(res.subs(self.subs).subs(self.atmospheric_basis.substitutions).subs(extra_subs)) else: return 0 @@ -422,14 +455,23 @@ class OceanicSymbolicInnerProducts(OceanicInnerProducts): Parameters ---------- - params: ~params.params.QgParams or list + params: ~.params.QgParams or list An instance of model's parameters object or a list in the form [aspect_ratio, ocean_basis, atmospheric_basis]. If a list is provided, `aspect_ratio` is the aspect ratio of the domain, `ocean_basis` is a SymbolicBasis object with the modes of the ocean, and `atmospheric_basis` is either a SymbolicBasis with the modes of the atmosphere or `None` if there is no atmosphere. + stored: bool, optional + Indicate if the inner product must be stored or computed on the fly. Default to `True` inner_product_definition: None or InnerProductDefinition, optional The definition of the inner product being used. If `None`, use the canonical StandardInnerProductDefinition object. Default to `None`. + interaction_inner_product_definition: None or InnerProductDefinition, optional + The definition of the inner product being used for the interaction with the other components, i.e. to compute the inner products with the other component base of funcitons. + If `None`, use the `inner_product_definition` provided. + Default to `None`. + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. Attributes ---------- @@ -441,13 +483,17 @@ class OceanicSymbolicInnerProducts(OceanicInnerProducts): Object holding the symbolic modes of the atmosphere (or `None` if there is no atmosphere). connected_to_atmosphere: bool Indicate if the ocean is connected to an atmosphere. + stored: bool + Indicate if the inner product must be stored or computed on the fly. ip: InnerProductDefinition Object defining the inner product. + iip: InnerProductDefinition + Object defining the interaction inner product. subs: list(tuple) List of 2-tuples containing the substitutions to be made with the functions after the inner products symbolic computation. """ - def __init__(self, params=None, inner_product_definition=None, stored=True, num_threads=None): + def __init__(self, params=None, stored=True, inner_product_definition=None, interaction_inner_product_definition=None, num_threads=None): OceanicInnerProducts.__init__(self) @@ -476,6 +522,11 @@ def __init__(self, params=None, inner_product_definition=None, stored=True, num_ else: self.ip = inner_product_definition + if interaction_inner_product_definition is None: + self.iip = self.ip + else: + self.iip = interaction_inner_product_definition + self.stored = stored if stored: self.compute_inner_products(num_threads) @@ -498,6 +549,9 @@ def connect_to_atmosphere(self, atmosphere_basis, num_threads=None): ---------- atmosphere_basis: SymbolicBasis or AtmosphericSymbolicInnerProducts Basis of function of the atmosphere or a symbolic atmospheric inner products object containing the basis. + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. """ if isinstance(atmosphere_basis, AtmosphericSymbolicInnerProducts): @@ -515,7 +569,7 @@ def connect_to_atmosphere(self, atmosphere_basis, num_threads=None): self._W = sp.zeros((self.noc, natm), dtype=float, format='dok') # K inner products - l = [[(i, j), self.ip.ip_lap, (self._phi(i), self._F(j))] for i in range(self.noc) + l = [[(i, j), self.iip.ip_lap, (self._phi(i), self._F(j))] for i in range(self.noc) for j in range(natm)] result = pool.map(apply, l) @@ -525,7 +579,7 @@ def connect_to_atmosphere(self, atmosphere_basis, num_threads=None): .subs(self.atmospheric_basis.substitutions)) # W inner products - l = [[(i, j), self.ip.symbolic_inner_product, (self._phi(i), self._F(j))] for i in range(self.noc) + l = [[(i, j), self.iip.symbolic_inner_product, (self._phi(i), self._F(j))] for i in range(self.noc) for j in range(natm)] result = pool.map(apply, l) @@ -540,6 +594,14 @@ def connect_to_atmosphere(self, atmosphere_basis, num_threads=None): pool.terminate() def compute_inner_products(self, num_threads=None): + """Function computing and storing all the inner products at once. + + Parameters + ---------- + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. + """ self._M = sp.zeros((self.noc, self.noc), dtype=float, format='dok') self._U = sp.zeros((self.noc, self.noc), dtype=float, format='dok') @@ -652,7 +714,7 @@ def K(self, i, j): """Function to commpute the forcing of the ocean by the atmosphere: :math:`K_{i,j} = (\phi_i, \\nabla^2 F_j)`.""" if self.connected_to_atmosphere: if not self.stored: - res = self.ip.ip_lap(self._phi(i), self._F(j)) + res = self.iip.ip_lap(self._phi(i), self._F(j)) return float(res.subs(self.subs).subs(self.oceanic_basis.substitutions).subs(self.atmospheric_basis.substitutions)) else: return self._K[i, j] @@ -663,7 +725,7 @@ def W(self, i, j): """Function to compute the short-wave radiative forcing of the ocean: :math:`W_{i,j} = (\phi_i, F_j)`.""" if self.connected_to_atmosphere: if not self.stored: - res = self.ip.symbolic_inner_product(self._phi(i), self._F(j)) + res = self.iip.symbolic_inner_product(self._phi(i), self._F(j)) return float(res.subs(self.subs).subs(self.atmospheric_basis.substitutions).subs(self.oceanic_basis.substitutions)) else: return self._W[i, j] @@ -677,14 +739,23 @@ class GroundSymbolicInnerProducts(GroundInnerProducts): Parameters ---------- - params: ~params.params.QgParams or list + params: ~.params.QgParams or list An instance of model's parameters object or a list in the form [aspect_ratio, ground_basis, atmospheric_basis]. If a list is provided, `aspect_ratio` is the aspect ratio of the domain, `ground_basis` is a SymbolicBasis object with the modes of the ground, and `atmospheric_basis` is either a SymbolicBasis with the modes of the atmosphere or `None` if there is no atmosphere. + stored: bool, optional + Indicate if the inner product must be stored or computed on the fly. Default to `True` inner_product_definition: None or InnerProductDefinition, optional The definition of the inner product being used. If `None`, use the canonical StandardInnerProductDefinition object. Default to `None`. + interaction_inner_product_definition: None or InnerProductDefinition, optional + The definition of the inner product being used for the interaction with the other components, i.e. to compute the inner products with the other component base of funcitons. + If `None`, use the `inner_product_definition` provided. + Default to `None`. + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. Attributes ---------- @@ -696,13 +767,17 @@ class GroundSymbolicInnerProducts(GroundInnerProducts): Object holding the symbolic modes of the atmosphere (or `None` if there is no atmosphere). connected_to_atmosphere: bool Indicate if the ground is connected to an atmosphere. + stored: bool + Indicate if the inner product must be stored or computed on the fly. ip: InnerProductDefinition Object defining the inner product. + iip: InnerProductDefinition + Object defining the interaction inner product. subs: list(tuple) List of 2-tuples containing the substitutions to be made with the functions after the inner products symbolic computation. """ - def __init__(self, params=None, inner_product_definition=None, stored=True, num_threads=None): + def __init__(self, params=None, stored=True, inner_product_definition=None, interaction_inner_product_definition=None, num_threads=None): GroundInnerProducts.__init__(self) @@ -731,7 +806,12 @@ def __init__(self, params=None, inner_product_definition=None, stored=True, num_ else: self.ip = inner_product_definition + if interaction_inner_product_definition is None: + self.iip = self.ip + else: + self.iip = interaction_inner_product_definition self.stored = stored + if stored: self.compute_inner_products(num_threads) @@ -753,6 +833,9 @@ def connect_to_atmosphere(self, atmosphere_basis, num_threads=None): ---------- atmosphere_basis: SymbolicBasis or AtmosphericSymbolicInnerProducts Basis of function of the atmosphere or a symbolic atmospheric inner products object containing the basis. + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. """ if isinstance(atmosphere_basis, AtmosphericSymbolicInnerProducts): atmosphere_basis = atmosphere_basis.atmospheric_basis @@ -769,7 +852,7 @@ def connect_to_atmosphere(self, atmosphere_basis, num_threads=None): self._W = sp.zeros((self.ngr, natm), dtype=float, format='dok') # W inner products - l = [[(i, j), self.ip.symbolic_inner_product, (self._phi(i), self._F(j))] for i in range(self.ngr) + l = [[(i, j), self.iip.symbolic_inner_product, (self._phi(i), self._F(j))] for i in range(self.ngr) for j in range(natm)] result = pool.map(apply, l) @@ -783,6 +866,14 @@ def connect_to_atmosphere(self, atmosphere_basis, num_threads=None): pool.terminate() def compute_inner_products(self, num_threads=None): + """Function computing and storing all the inner products at once. + + Parameters + ---------- + num_threads: int or None, optional + Number of threads to use to compute the symbolic inner products. If `None` use all the cpus available. + Default to `None`. + """ self._U = sp.zeros((self.ngr, self.ngr), dtype=float, format='dok') @@ -802,7 +893,7 @@ def compute_inner_products(self, num_threads=None): @property def ngr(self): - """Number of oceanic modes.""" + """Number of ground modes.""" return len(self.ground_basis) # !-----------------------------------------------------! @@ -861,7 +952,7 @@ def W(self, i, j): """Function to compute the short-wave radiative forcing of the ocean: :math:`W_{i,j} = (\phi_i, F_j)`.""" if self.connected_to_atmosphere: if not self.stored: - res = self.ip.symbolic_inner_product(self._phi(i), self._F(j)) + res = self.iip.symbolic_inner_product(self._phi(i), self._F(j)) return float(res.subs(self.subs).subs(self.atmospheric_basis.substitutions).subs(self.ground_basis.substitutions)) else: return self._W[i, j] diff --git a/params/params.py b/params/params.py index 957df47..a719022 100644 --- a/params/params.py +++ b/params/params.py @@ -408,6 +408,8 @@ def set_insolation(self, value, pos=None): Indicate in which component to set the `value`. """ + # TODO: - check for the dimensionality of the arguments + if isinstance(value, (float, int)) and pos is not None and self.C is not None: self.C[pos] = Parameter(value, units='[W][m^-2]', scale_object=self._scale_params, description="spectral component "+str(pos+1)+" of the short-wave radiation of the atmosphere", @@ -445,6 +447,8 @@ def set_thetas(self, value, pos=None): Indicate in which component to set the `value`. """ + # TODO: - check for the dimensionality of the arguments + if isinstance(value, (float, int)) and pos is not None and self.thetas is not None: self.thetas[pos] = Parameter(value, scale_object=self._scale_params, description="spectral components "+str(pos+1)+" of the temperature profile", @@ -622,7 +626,7 @@ def __init__(self, scale_params, dic=None): self.set_params(dic) - def set_orography(self, value, pos=None): + def set_orography(self, value, pos=None, basis="atmospheric"): """Function to define the spectral decomposition of the orography profile :math:`h_k` (:attr:`~.GroundParams.hk`). @@ -633,8 +637,16 @@ def set_orography(self, value, pos=None): If an iterable is provided, create a vector of spectral decomposition parameters corresponding to it. pos: int, optional Indicate in which component to set the `value`. + basis: str, optional + Indicate which basis should be used to decompose the orography. Can be either `atmospheric`, `oceanic` or `ground`. + Default to `atmospheric`. """ + # TODO: - check for the dimensionality of the arguments + # - check that inner products are symbolic if basis is not 'atmospheric' + + self.orographic_basis = basis + if isinstance(value, (float, int)) and pos is not None and self.hk is not None: self.hk[pos] = Parameter(value, scale_object=self._scale_params, description="spectral components "+str(pos+1)+" of the orography", @@ -713,6 +725,8 @@ def set_insolation(self, value, pos=None): Indicate in which component to set the `value`. """ + # TODO: - check for the dimensionality of the arguments + if isinstance(value, (float, int)) and pos is not None and self.C is not None: self.C[pos] = Parameter(value, units='[W][m^-2]', scale_object=self._scale_params, description="spectral component "+str(pos+1)+" of the short-wave radiation of the ground",