From db8e7d32a2c8ce96b01eb887eeceb8b79da2fac1 Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Thu, 21 May 2026 14:15:24 +0100 Subject: [PATCH 1/2] tangent linear documentation --- .../tangent_linear_section/index.rst | 23 ++ .../tangent_linear_coding.rst | 98 ++++++++ .../tangent_linear_components.rst | 229 ++++++++++++++++++ .../tangent_linear_lin_state.rst | 16 ++ .../tangent_linear_overview.rst | 39 +++ .../tangent_linear_tests.rst | 64 +++++ .../tangent_linear_timestepping.rst | 188 ++++++++++++++ .../tl_ffsl_transport.rst | 172 +++++++++++++ 8 files changed, 829 insertions(+) create mode 100644 documentation/source/science_guide/tangent_linear_section/index.rst create mode 100644 documentation/source/science_guide/tangent_linear_section/tangent_linear_coding.rst create mode 100644 documentation/source/science_guide/tangent_linear_section/tangent_linear_components.rst create mode 100644 documentation/source/science_guide/tangent_linear_section/tangent_linear_lin_state.rst create mode 100644 documentation/source/science_guide/tangent_linear_section/tangent_linear_overview.rst create mode 100644 documentation/source/science_guide/tangent_linear_section/tangent_linear_tests.rst create mode 100644 documentation/source/science_guide/tangent_linear_section/tangent_linear_timestepping.rst create mode 100644 documentation/source/science_guide/tangent_linear_section/tl_ffsl_transport.rst diff --git a/documentation/source/science_guide/tangent_linear_section/index.rst b/documentation/source/science_guide/tangent_linear_section/index.rst new file mode 100644 index 000000000..7b63c984c --- /dev/null +++ b/documentation/source/science_guide/tangent_linear_section/index.rst @@ -0,0 +1,23 @@ +.. ----------------------------------------------------------------------------- + (c) Crown copyright Met Office. All rights reserved. + The file LICENCE, distributed with this code, contains details of the terms + under which the code may be used. + ----------------------------------------------------------------------------- + +.. _tangent_linear_section_index: + +Linear Model Section +======================= + +This describes the details of the tangent-linear model code for GungHo. + +.. toctree:: + :maxdepth: 2 + + tangent_linear_overview + tangent_linear_coding + tangent_linear_timestepping + tangent_linear_components + tl_ffsl_transport + tangent_linear_lin_state + tangent_linear_tests diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_coding.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_coding.rst new file mode 100644 index 000000000..5ef2c5016 --- /dev/null +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_coding.rst @@ -0,0 +1,98 @@ +.. _tangent_linear_coding: + +Tangent Linear Coding +===================== + +To obtain the tangent linear code, we need to differentiate the +nonlinear model code. + +First derivative +---------------- + +For the most simple example, if the nonlinear model is :math:`y = f(x)`, +then the tangent linear model is + +.. math:: \delta y = \frac{df}{dx} \delta x + +For example, if the nonlinear model code is + +:: + + y = x**2 + +then :math:`f(x)=x^2`, and the first derivative is :math:`f'(x) = 2x`. +So the linear model code is + +:: + + y = 2 * ls_x * x + +where now ``y``, and ``x`` are the perturbations, and ``ls_x`` is the +linearisation state. + +Note that we could rewrite this in other variable names as e.g. +``delta_y = 2 * x * delta_x``, but in the LFRic code we have chosen to +keep the variable names for the perturbations identical to the nonlinear +model variable names, and to introduce ``ls_`` to indicate the +linearisation state. This linearisation state is produced from a +separate nonlinear model run. + +The chain rule +-------------- + +For a more complicated expression, then the chain rule of +differentiation is required + +.. math:: + + \begin{aligned} + y & = f(g(x)) \\ + \delta y & = \frac{df}{dx} \delta x \\ + & = \frac{df}{dg} \frac{dg}{dx} \delta x + \end{aligned} + +For example, if the nonlinear model code is + +:: + + y = (sin v)**2 + +so + +.. math:: + + \begin{aligned} + f(g) & = g^2 \\ + g(x) & = sin(x)\\ + f'(g) & = 2 g = 2 sin(x) \\ + g'(x) & = cos(x) + \end{aligned} + +then the linear model code is + +:: + + y = 2 * sin(ls_v) * cos(ls_v) * v + +Multivariable calculus +---------------------- + +When dealing with more than one variable, then the rules need to be +extended to multi-variable calculus. + +If the nonlinear model is + +.. math:: y = f(x,z) + +the linear model is + +.. math:: \delta y = \frac{\partial f}{\partial x} \delta x + \frac{\partial f}{\partial z} \delta z + +and the multivariate chain rule gives + +.. math:: + + \begin{aligned} + y & = f(g_1(x_1), ..., g_n(x_n)) \\ + \delta y & = \sum_i \frac{\partial f}{\partial g_i} \frac{d g_i}{d x_i} \delta x_i + \end{aligned} diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_components.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_components.rst new file mode 100644 index 000000000..b9618fb8b --- /dev/null +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_components.rst @@ -0,0 +1,229 @@ +.. _tangent_linear_components: + +Components +===================== + +Kinetic energy gradient +----------------------- + +The nonlinear model is + +.. math:: + + \begin{aligned} + f & = \nabla (\frac{1}{2} u^2) \\ + & = g(z(u)) \\ + \text{where} \quad z(u) & = \frac{1}{2} u^2 \\ + g(z) & = \nabla(z) + \end{aligned} + +The linear model is, + +.. math:: + + \begin{aligned} + \delta f & = \frac{d g}{d z} \frac{ dz}{du} \delta u\\ + & = \nabla (u \delta u) + \end{aligned} + +Note that an alternative way to derive this is to form the difference +between the perturbed and unperturbed, and remove the small terms + +.. math:: + + \begin{aligned} + \delta f & = f(\bar{u}+ \delta u) - f(\bar{u}) \\ + & = \frac{1}{2} \nabla (\bar{u} + \delta u)^2 - \frac{1}{2} \nabla \bar{u}^2 \\ + & = \frac{1}{2} \nabla (\bar{u}^2 + 2\bar{u} \delta u + \delta u^2) - \frac{1}{2} \nabla \bar{u}^2 \\ + & \approx \nabla (\bar{u} \delta u) + \end{aligned} + +Potential temperature +--------------------- + +The nonlinear model is + +.. math:: \theta_e = \theta (m_g/m_t) + +The linear model is, + +.. math:: + + \begin{aligned} + \delta \theta_e & = \frac{\partial \theta_e}{\partial \theta} \delta \theta + + \frac{\partial \theta_e}{\partial m_g} \delta m_g + + \frac{\partial \theta_e}{\partial m_t} \delta m_t \\ + & = \frac{m_g}{m_t} \delta \theta + + \frac{\theta}{m_t} \delta m_g + - \frac{\theta m_g}{m_t^2} \delta m_t + \end{aligned} + +using + +.. math:: \frac{d (1/m_t) }{dx} = \frac{d (m_t^{-1}) }{dx} = -1 (m_t^{-2}) = - \frac{1}{m_t^2} + +This can be simplified, by substituting in the equation for the +nonlinear model :math:`\theta_e = \theta (m_g/m_t)` + +.. math:: \delta \theta_e = \theta_e \left( \frac{ \delta \theta}{\theta} + \frac{\delta m_g}{m_g} - \frac{\delta m_t}{m_t} \right) + +Logspace +-------- + +The nonlinear model is + +.. math:: + + \begin{aligned} + p & = \Pi_i |\rho_i|^{c_i} \\ + & = \Pi_i u_i(y_i(z_i(\rho_i))) \\ + & = \Pi_i u_i \\ + \text{where} \quad u_i & = y_i ^{c_i} \\ + y_i & = | \rho_i | \\ + & = (z_i) ^{1/2} \\ + z_i & = \rho_i^2 + \end{aligned} + +The derivatives are + +.. math:: + + \begin{aligned} + \frac{ \partial p}{\partial u_i} & = \Pi_{j \neq i} u_j = \Pi_{j \neq i} | \rho_j |^{c_j} \\ + \frac{ \partial u}{\partial y_i} & = c_i y_i ^{c_i -1} = c_i | \rho_i | ^{c_i -1} \\ + \frac{ \partial y}{\partial z_i} & = \frac{1}{2} z_i^{-1/2} = \frac{1}{2} \frac{1}{| \rho_i | } \\ + \frac{ \partial z}{\partial \rho_i} & = 2 \rho_i + \end{aligned} + +Therefore the linear model is, + +.. math:: + + \begin{aligned} + \delta p & = \sum_i \frac{ \partial p}{\partial u_i} \frac{ \partial u_i}{\partial y_i}\frac{ \partial y_i}{\partial z_i} \frac{ \partial z_i}{\partial \rho_i} \delta \rho_i \\ + & = \sum_i \left( \left( \Pi_{j \neq i} | \rho_j |^{c_j} \right) c_i | \rho_i | ^{c_i -1} \frac{\rho_i}{| \rho_i |} \delta \rho_i \right) + \end{aligned} + +This can be simplified, by using a rearrangement of the nonlinear model + +.. math:: + + \begin{aligned} + p & = \Pi_i |\rho_i|^{c_i} \\ + & = |\rho_i|^{c_i} \Pi_{j \neq i} | \rho_j|^{c_j} \\ + \text{so } \quad \Pi_{j \neq i} | \rho_j|^{c_j} & = \frac{p}{|\rho_i|^{c_i}} \\ + \end{aligned} + +Substituting this into the linear model, gives + +.. math:: + + \begin{aligned} + \delta p & = p \sum_i c_i \frac{\rho_i }{ | \rho_i |^2} \delta \rho_i \\ + & = p \sum_i \frac{c_i }{ \rho_i} \delta \rho_i \\ + \end{aligned} + +Newton iteration of the Equation of State +----------------------------------------- + +The equation of state is applied iteratively in the semi-implicit +timestepping. The nonlinear model rhs is + +.. math:: + + \begin{aligned} + N & = 1/E -1 + \end{aligned} + +where + +.. math:: + + \begin{aligned} + E & = \frac{p_0}{R} \frac{\Pi^\gamma}{\theta \rho} \\ + \gamma & = \frac{1-\kappa}{\kappa} + \end{aligned} + +The derivatives are: + +.. math:: + + \begin{aligned} + \frac{\partial N}{\partial E} & = - E^{-2} \\ + \frac{\partial E}{\partial \Pi} & = \frac{p_0}{R} \gamma \Pi^{\gamma-1}\theta^{-1} \rho^{-1} = E \gamma \Pi^{-1} \\ + \frac{\partial E}{\partial \rho} & = - \frac{p_0}{R} \Pi^\gamma \rho^{-2} \theta^{-1} = - E \rho^{-1} \\ + \frac{\partial E}{\partial \theta} & = - \frac{p_0}{R} \Pi^\gamma \rho^{-1} \theta^{-2} = - E \theta^{-1} + \end{aligned} + +The linear model is + +.. math:: + + \begin{aligned} + \delta N & = \frac{\partial N}{\partial E} \left( \frac{ \partial E}{\partial \Pi} \delta \Pi + \frac{ \partial E}{\partial \rho} \delta \rho + \frac{ \partial E}{\partial \theta} \delta \theta \right) \\ + & = - \frac{1}{E} \left( \frac{1-\kappa}{\kappa} \frac{\delta \Pi}{\Pi} - \frac{\delta \rho}{\rho} - \frac{\delta \theta}{\theta} \right) + \end{aligned} + +Method of Lines (MoL) and Runge-Kutta Transport +----------- + +The nonlinear model is: + +.. math:: + + \begin{aligned} + u^{(i)} & = u^n + \Delta t \sum_{k=0}^{i-1} \beta_{i,k} f(u^{(k)}) \quad i=1,\ldots,m \\ + u^{n+1} & = u^{(m)} + \end{aligned} + +In ``rk_alg_timestep``, the :math:`f(u^{(i)})` corresponds to +``rhs_prediction``, and :math:`\sum_{k=0}^{i-1} \beta_{i,k} f(u^{(k)})` +corresponds to rhs. (There is also an extra mass matrix inversion, so +they are not directly equal). + +The tangent linear model comprises two steps. First the linearisation +state is computed: + +.. math:: \bar{u}^{(i)} = \bar{u}^n + \Delta t \sum_{k=0}^{i-1} \beta_{i,k} f(\bar{u}^{(k)}) \quad i=1,\ldots,m + +then the perturbations are computed: + +.. math:: + + \begin{aligned} + u'^{(i)} & = u'^n + \Delta t \sum_{k=0}^{i-1} \beta_{i,k} f'(u'^{(k)},\bar{u}^{(k)}) \quad i=1,\ldots,m \\ + u'^{n+1} & = u'^{(m)} + \end{aligned} + +In the nonlinear model, the method of lines transport code has switches +to ensure that the transport scheme is upwinding. In the linear model, +we use the same switches as the nonlinear model. This means that the +direction is based on the linearisation state wind rather than the +perturbation wind. + +In the nonlinear model, + +:: + + direction = wind(map_w2(df) + k )*v_dot_n(df) + if ( direction > 0.0_r_def ) then + +In the linear model, + +:: + + direction = ls_wind(map_w2(df) + k )*v_dot_n(df) + if ( direction > 0.0_r_def ) then + +Similarly, where + +:: + + sign_X(unit_wind,1.0_r_def, advecting_wind) + +the corresponding tangent linear uses + +:: + + sign_X(unit_wind, 1.0_r_def, ls_advecting_wind) + diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_lin_state.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_lin_state.rst new file mode 100644 index 000000000..514478d2c --- /dev/null +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_lin_state.rst @@ -0,0 +1,16 @@ +.. _tangent_linear_lin_state: + +Linearisation state +=================== + +If the linear model is run using idealized initial conditions, then it +is possible to define an analytical linearisation state. Otherwise, the +linearisation state needs to be read and updated from file. The +linearisation state is continually updated over the linear model run +because it is not fixed in time; the linearisation state is the +nonlinear model trajectory. + +The linearisation state is read and updated from file in the same way as +the ancillaries, and the lateral-boundary-conditions. This means that +the data is read with XIOS using a time-axis and on every timestep, the +data is updated and interpolated in time. diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_overview.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_overview.rst new file mode 100644 index 000000000..67e308ae2 --- /dev/null +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_overview.rst @@ -0,0 +1,39 @@ +.. _tangent_linear_overview: + +Tangent Linear Model Overview +======================= + +In variational data assimilation, the aim is to compute the optimal +value of the initial state :math:`\mathbf{x}_0` based on the +minimization of a cost function: + +.. math:: J(\mathbf{x}_0) = (\mathbf{x}_0 - \mathbf{x}_b)^T \mathbf{B}^{-1} (\mathbf{x}_0 - \mathbf{x}_b) + \sum_t ( \mathbf{y}_t - \mathbf{H}_t N_t(\mathbf{x}_0))^T \mathbf{R}^{-1} ( \mathbf{y}_t - \mathbf{H}_t N_t(\mathbf{x}_0)) + +where :math:`t` is the time, :math:`\mathbf{x}_b` is the background +state, :math:`\mathbf{y}_t` are the observations, :math:`\mathbf{H}_t` +is the observation operator (which we can assume is linear here), and +:math:`N_t` is the nonlinear model that evolves the state +:math:`\mathbf{x}_0` from time :math:`0` to time :math:`t`. +:math:`\mathbf{B}` and :math:`\mathbf{R}` are the error covariance +matrices. + +As :math:`N` is nonlinear, there is no unique solution. Therefore, the +problem is rewritten in incremental form, by rewriting the initial state +in terms of a linearisation state plus a perturbation, +:math:`\mathbf{x}_0 = \tilde{\mathbf{x}} + \mathbf{\delta x}`, so that +on the first iteration, when :math:`\tilde{\mathbf{x}} = \mathbf{x}_b`, + +.. math:: J(\mathbf{\delta x}_0) = (\mathbf{\delta x}_0)^T B^{-1} (\mathbf{\delta x}_0) + \sum_t ( \mathbf{\delta y}_t - \mathbf{H}_t \mathbf{L}_t\mathbf{\delta x}_0)^T \mathbf{R}^{-1} ( \mathbf{\delta y}_t - \mathbf{H}_t \mathbf{L}_t \mathbf{\delta x}_0) + +where :math:`\mathbf{\delta x}_0 = \mathbf{x}_0 - \mathbf{x}_b`, and +:math:`\mathbf{\delta y}_t = \mathbf{y}_t - \mathbf{H}_t N_t (x_b)`, and +we have used the Taylor series expansion to approximate the nonlinearly +evolved perturbation, + +.. math:: N( \mathbf{x}_b + \mathbf{\delta x} ) - N(\mathbf{x}_b) \approx \mathbf{L\delta x} + +where the linear model +:math:`\mathbf{L}= \mathbf{L}(\tilde{\mathbf{x}})`, known as the tangent +linear model, is a function of the linearisation state +:math:`\tilde{\mathbf{x}}`, and is the Jacobian or first derivative of +the nonlinear model. diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_tests.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_tests.rst new file mode 100644 index 000000000..37517311e --- /dev/null +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_tests.rst @@ -0,0 +1,64 @@ +.. _tangent_linear_tests: + + +Tangent Linear Tests +==================== + +Taylor remainder convergence test +--------------------------------- + +By definition, the tangent linear model is the Jacobian of the nonlinear +model, and therefore the Taylor series expansion can be used to form a +test for the code. + +The Taylor series expansion of :math:`N(\mathbf{x})` is + +.. math:: N(\mathbf{x} + \gamma \mathbf{\delta x} ) = N(\mathbf{x}) + \gamma \mathbf{L \delta x} + \mathcal{O}(\gamma^2) + +This is used to calculate the norm of the difference between the +nonlinearly evolved perturbation +:math:`N(\mathbf{x} + \gamma \mathbf{\delta x} ) - N(\mathbf{x})` and +the linearly evolved perturbation :math:`\gamma \mathbf{L \delta x}` + +.. math:: E_n(\gamma) = || N(\mathbf{x} + \gamma \mathbf{\delta x}) - N(\mathbf{x}) - \gamma \mathbf{L \delta x} || = \mathcal{O}(\gamma^2) + +where the norm could be the L2-norm, +:math:`||\mathbf{x}|| = (\mathbf{x}^T\mathbf{x})^{1/2}`. + +The norm :math:`E_n(\gamma)` is known as the linearization error, and it +tends to zero as the size of the perturbation :math:`\gamma` tends to +zero, at order :math:`\mathcal(O)(\gamma^2)`. + +If we examine the series of :math:`E_n` for values of :math:`\gamma` +that are successively halved, we can evaluate the convergence rate. For +example, with the values of :math:`E_n(\gamma)` and :math:`E_n(2\gamma)` +we have + +.. math:: E_n(2\gamma) / E_n(\gamma) = \mathcal{O}(4 \gamma^2) / \mathcal{O}(\gamma^2) = 4 + +To create a test which gives a simple pass or fail, two values of the +sequence :math:`E_n` are chosen, and the ratio is compared to the value +:math:`4`. If the value of the ratio is :math:`4`, to within a specified +tolerance, then the code passes the test. + +Implementation +-------------- + +The Taylor remainder convergence test is implemented in the LFRic code +through the use of integration tests. In the linear app integration +tests, the convergence test is applied to small sections of code. The +process for each test is + +- Initialise the nonlinear model, linear model and linearisation state. + +- Run the nonlinear model using the linearisation state as initial + conditiosn. + +- Run the nonlinear model with a perturbed linearisation state, and run + the linear model with the perturbation, and calculate + :math:`E_n(\gamma)`. + +- Halve the size of the perturbation, and repeat. + +- Compare the ratio of the values of :math:`E_n` with :math:`4`, to a + specified tolerance. diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_timestepping.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_timestepping.rst new file mode 100644 index 000000000..98f2c5cbe --- /dev/null +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_timestepping.rst @@ -0,0 +1,188 @@ +.. _tangent_linear_timestepping: + +Iterated (Semi-) Implicit Timestepping +==================== + +Nonlinear Equations +------------------- + +The nonlinear equations are + +.. math:: + + \begin{aligned} + \frac{\partial u}{\partial t} &= - u . \nabla u - 2 \Omega \times u - \nabla \Phi - c_P \theta \nabla \Pi \\ + \frac{\partial \rho}{\partial t} & = - \nabla . (\rho u) \\ + \frac{\partial \theta}{\partial t} & = - u. \nabla \theta + \end{aligned} + +together with the nonlinear equation of state: + +.. math:: \Pi ^{\left( \frac{1-\kappa}{\kappa} \right)} = \frac{R}{p_0} \rho \theta + +Linear Equations +---------------- + +The linearized equations are: + +.. math:: + + \begin{aligned} + \frac{\partial \delta u}{\partial t} &= - \bar{u} . \nabla \delta u - \delta {u} . \nabla \bar{ u} - 2 \Omega \times \delta u - c_P \delta \theta \nabla \bar{\Pi} - c_P \bar{\theta} \nabla \delta \Pi \\ + \frac{\partial \delta \rho}{\partial t} & = - \nabla . (\bar{\rho} \delta u) - \nabla . (\delta \rho \bar{u}) \\ + \frac{\partial \delta \theta}{\partial t} & = - \bar{u}. \nabla \delta \theta - \delta u. \nabla \bar{\theta} + \end{aligned} + +together with the linear equation of state: + +.. math:: \frac{ \delta \Pi}{\bar{\Pi}} = \frac{\kappa}{1-\kappa} \left( \frac{\delta \rho}{\bar{\rho}} + \frac{\delta \theta}{\bar{\theta}} \right) + + +Iterated (Semi-) implicit Timestepping +-------------------------- + +In the iterated-implicit timestepping, the equations are solved iteratively +using a Newton method. If the full nonlinear discrete equations are +written as + +.. math:: \mathcal{R} (\mathbf{x}^{n+1}) =0 + +then this can be solved iteratively as + +.. math:: + + \begin{aligned} + \mathbf{x}^{n+1}_{(k+1)} & = \mathbf{x}^{n+1}_{(k)} + \mathbf{x}'_{(k)} \\ + \mathcal{L} (\mathbf{x}^*) \mathbf{x}'_{(k)} & = - \mathcal{R} ( \mathbf{x}^{n+1}_{(k)}) + \end{aligned} + +where :math:`\mathcal{L}` is an approximation to the Jacobian and +:math:`\mathbf{x}^*` is the reference state. + +The reference state :math:`u^*` is zero whilst the other reference state +variables are set to :math:`\mathbf{x}^n`. Relaxation parameters +:math:`\tau <1` are also introduced to give extra weight to the diagonal +terms (and hence similar to solving the resulting linear equation using +the Jacobi method or a damped Newton method). + +This gives the approximation to the Jacobian as + +.. math:: + + \mathcal{L}(\mathbf{x}^*) \mathbf{x}' = \left\{ + \begin{aligned} + u' & + \tau_u \Delta t c_P (\theta' \nabla \Pi^* + \theta^* \nabla \Pi ') \\ + \rho' & + \tau_\rho \Delta t \nabla.(\rho^* u') \\ + \theta' & + \tau_\theta \Delta t u' \nabla. \nabla \theta^* \\ + & \frac{1-\kappa}{\kappa} \frac{\Pi'}{\Pi^*} - \frac{1}{E} \frac{\rho'}{\rho^*} - \frac{1}{E} \frac{\theta'}{\theta^*} + \end{aligned} + \right. + +where + +.. math:: E = \frac{p_0}{R} \frac{\Pi^{*(1-\kappa)/\kappa}}{\theta^* \rho^*} + +We can proceed in a similar way with the linear model, so that + +.. math:: + + \begin{aligned} + \delta \mathbf{x}^{n+1}_{(k+1)} & = \delta \mathbf{x}^{n+1}_{(k)} + \delta \mathbf{x}'_{(k)} \\ + \mathcal{L} (\mathbf{x}^*) \delta \mathbf{x}'_{(k)} & = - R ( \mathbf{x}_{(k)}) \delta \mathbf{x}_{(k)} + \end{aligned} + +where :math:`R` is the linearized version of :math:`\mathcal{R}`. + +As the operator :math:`\mathcal{L}` is already linear, this does not +need to be linearized. + +Discretisation +-------------------- + +The nonlinear model evolves the linearisation state +:math:`\bar{\mathbf{x}}`. The aim is to find the state +:math:`\bar{\mathbf{x}}^n` at time :math:`t_{n+1}`. This is solved +iteratively, by letting +:math:`{\bar{\mathbf{x}}}^{n+1}_{(k+1)} = \bar{\mathbf{x}}^{n+1}_{(k)} + \bar{\mathbf{x}}'_{(k+1)}` +where :math:`k` is the iteration count. + +The increment :math:`\bar{\mathbf{x}}'_{(k+1)}` solves the linear +equation + +.. math:: \mathcal{L}( \bar{\mathbf{x}}*) \bar{\mathbf{x}}'_{(k+1)} = \mathcal{R}^n(\bar{\mathbf{x}}^n) + \mathcal{R}^{n+1}(\bar{\mathbf{x}}^{n+1}_{(k)}) + \mathcal{R}^{adv}(\bar{\mathbf{x}}^n, \bar{\mathbf{x}}^{n+1}_{(k)}) + +where :math:`\mathcal{L}( \bar{\mathbf{x}}*)` is a linear operator based +on the reference state :math:`\bar{\mathbf{x}}*`. + +The procedure is: + +- | Let :math:`\bar{\mathbf{x}}^{n+1}_{(0)} = \bar{\mathbf{x}}^n` + | The first iteration is given by the state at the previous timestep. + +- | Set + :math:`\bar{\mathbf{x}}* \equiv (\bar{\rho}*, \bar{\theta}*, \bar{\Pi}* ) = (\bar{\rho}^n, \bar{\theta}^n, \bar{\Pi}^n)` + | The reference state is given by the state at the previous timestep + (for density, potential temperature and exner pressure only). + +- | Compute :math:`\mathcal{R}^n(\bar{\mathbf{x}}^n)` + | This is the RHS based on the previous timestep. + +- Loop k=1, K + + - | Compute :math:`\mathcal{R}^{n+1}(\bar{\mathbf{x}}^{n+1}_{(k)})` + and + :math:`\mathcal{R}^{adv}(\bar{\mathbf{x}}^n, \bar{\mathbf{x}}^{n+1}_{(k)})` + | This is the RHS and advective terms based on the current timestep. + + - Solve the linear equation for :math:`\bar{\mathbf{x}}'_{(k+1)}` and + hence for :math:`\bar{\mathbf{x}}^{n+1}_{(k+1)}`. + +The linear model evolves the perturbations :math:`\mathbf{x}`. To find +the state at time :math:`t_{n+1}`, we follow the same procedure as for +the nonlinear model, but using the tangent linear versions of the right +hand side in the equation. + +Let +:math:`{\mathbf{x}}^{n+1}_{(k+1)} = \mathbf{x}^{n+1}_{(k)} + \mathbf{x}'_{(k+1)}` +where :math:`\mathbf{x}'_{(k+1)}` solves + +.. math:: \mathcal{L} ( \bar{\mathbf{x}}*) \mathbf{x}'_{(k+1)} = R^n(\mathbf{x}^n, \bar{\mathbf{x}}^n) + R^{n+1}(\mathbf{x}^{n+1}_{(k)},\bar{\mathbf{x}}^{n+1}_{(k)}) + R^{adv}(\mathbf{x}^n, \mathbf{x}^{n+1}_{(k)},\bar{\mathbf{x}}^n, \bar{\mathbf{x}}^{n+1}_{(k)}) + +where :math:`R^n` is the tangent linear of :math:`\mathcal{R}^n` etc. +The tangent linear operators are a function of both the perturbation +:math:`\mathbf{x}` and the linearisation state :math:`\bar{\mathbf{x}}`. + +The procedure is: + +- | Let :math:`\mathbf{x}^{n+1}_{(0)} = \mathbf{x}^n` + | The first iteration is given by the perturbation at the previous + timestep. + +- | Set + :math:`\bar{\mathbf{x}}* \equiv (\bar{\rho}*, \bar{\theta}*, \bar{\Pi}* ) = (\bar{\rho}^n, \bar{\theta}^n, \bar{\Pi}^n)` + | The reference state is given by the linearisation state at the + previous timestep - so that the linear operator is identical to the + that in the nonlinear model. + +- Compute :math:`R^n(\mathbf{x}^n, \bar{\mathbf{x}}^n)` + +- Loop k=1, K + + - Compute + :math:`R^{n+1}(\mathbf{x}^{n+1}_{(k)},\bar{\mathbf{x}}^{n+1}_{(k)})` + and + :math:`R^{adv}(\mathbf{x}^n, \mathbf{x}^{n+1}_{(k)},\bar{\mathbf{x}}^n, \bar{\mathbf{x}}^{n+1}_{(k)})` + + - Solve for :math:`\mathbf{x}'_{(k+1)}` and hence for + :math:`\mathbf{x}^{n+1}_{(k+1)}`. + +The linear solver :math:`\mathcal{L} ( \bar{\mathbf{x}}*)` is exactly +the same as for the nonlinear model. + +The right hand side terms :math:`R^n`, :math:`R^{n+1}` and +:math:`R^{adv}` require the linearisation state at the same iteration. +Therefore, for an exact tangent linear semi-implicit algorithm, the +nonlinear semi-implicit step needs to be run first, and the intermediate +``ls_state`` values need to be stored and used in the following tangent +linear semi-implicit step. + diff --git a/documentation/source/science_guide/tangent_linear_section/tl_ffsl_transport.rst b/documentation/source/science_guide/tangent_linear_section/tl_ffsl_transport.rst new file mode 100644 index 000000000..c3b6ec7ee --- /dev/null +++ b/documentation/source/science_guide/tangent_linear_section/tl_ffsl_transport.rst @@ -0,0 +1,172 @@ +.. _tl_ffsl_transport: + +FFSL Transport +================ + +| The non-linear flux-form semi-Lagrangian (FFSL) transport scheme solves the conservative form transport + equation in the form + + .. math:: q^{n+1} = q^{n} - \Delta{t} F(u,q^n) + + where :math:`q` is the variable we are transporting, :math:`u` is the + advecting wind, :math:`\Delta{t}` is the time step, and :math:`F` is + the divergence of the mass flux + + .. math:: F(u,q) = \frac{1}{V} \nabla \cdot \tilde{f}(u,q) + + where the flux :math:`\tilde{f}`, is computed using the chosen FFSL + scheme. + +The non-linear scheme solves the advective form transport equation in +the form + +.. math:: + + \label{eq:advective} + q^{n+1}_{adv} = \frac{q^{n} - \Delta{t} F(u,q^n)}{\sigma - \Delta{t} F(u,\sigma)} + +where :math:`\sigma=1` is the unity field, and therefore +:math:`F(u,\sigma)` is the wind divergence. The unity field plays the +role of a psuedo-density. + +Tangent Linear Transport +------------------------- + +For the linear model we split the variables into a linear state (LS) +denoted with a bar, and a perturbation denoted with a prime + +.. math:: q = \bar{q} + q' \hspace{1cm} u = \bar{u}+u' + +The conservative transport equation becomes + +.. math:: \frac{\partial \bar{q}}{\partial t} + \frac{\partial q'}{\partial t} + \nabla \cdot( (\bar{u}+u')(\bar{q}+q') ) = 0 + +The product of perturbations is assumed to be zero, and so multiplying +out the divergence term gives + +.. math:: \frac{\partial \bar{q}}{\partial t} + \frac{\partial q'}{\partial t} + \nabla \cdot( \bar{u}\bar{q} ) + \nabla \cdot( \bar{u}q' ) + \nabla \cdot( u'\bar{q} ) = 0 + +The linear state conservative transport equation is + +.. math:: \frac{\partial \bar{q}}{\partial t} + \nabla \cdot( \bar{u}\bar{q} ) = 0 + +Leaving the conservative perturbation transport equation as + +.. math:: + + \label{eq:transport_equation} + \frac{\partial q'}{\partial t} + \nabla \cdot( \bar{u}q' ) + \nabla \cdot( u'\bar{q} ) = 0 + +Flux Computation +---------------------- + +The flux computation can be split into two parts, corresponding to the +two terms in equation +`[eq:transport_equation] <#eq:transport_equation>`__. + +The first term, corresponding to :math:`\bar{u}q'`, can be computed +using standard FFSL operators, provided the scheme is linear (i.e. no +flux limiting). This is the reconstruction of the perturbation field +with the ls wind, and we denote this :math:`\tilde{f}_1`. + +The second term, corresponding to :math:`u'\bar{q}`, uses the +perturbation departure distance from the ls departure point. This is +computed as, where the subscript :math:`d` signifies the ls wind +departure cell, + +.. math:: \tilde{f}_2(u',\bar{q}) = u' \bar{q}_d + +The constant reconstruction is to ensure that this term is linear in the +perturbation wind. + +We can now use the notation + +.. math:: + + \begin{aligned} + F_1(u,q) & = \frac{1}{V} \nabla \cdot \tilde{f}_1(u,q) \\ + F_2(u,q) & = \frac{1}{V} \nabla \cdot \tilde{f}_2(u,q) + \end{aligned} + +Advective Update +-------------------- + +Following `[eq:advective] <#eq:advective>`__, if we expand the variables +into linearisation state and perturbation parts we find the flux +divergences can be written as + +.. math:: F(u,q) = F(\bar{u},\bar{q}) + F_1(\bar{u},q') + F_2(u',\bar{q}) + +and, as :math:`\sigma = \bar{\sigma}`, + +.. math:: F(u,\sigma) = F(\bar{u},\sigma) + F_2(u',\sigma) + +and the advective form becomes + +.. math:: + + \label{eq:advective_linear} + q^{n+1}_{adv} = \frac{\bar{q} + q' - \Delta{t} F(\bar{u},\bar{q}) - \Delta{t}F_1(\bar{u},q') - \Delta{t}F_2(u',\bar{q}) }{\sigma - \Delta{t} F(\bar{u},\sigma)- \Delta{t} F_2(u',\sigma)} + +Let + +.. math:: + + \begin{aligned} + \bar{N} &= \bar{q} - \Delta{t} F(\bar{u},\bar{q}) \\ + N' &= q' - \Delta{t}F_1(\bar{u},q') - \Delta{t}F_2(u',\bar{q}) \\ + \bar{D} &= \sigma - \Delta{t} F(\bar{u},\sigma) \\ + D' = - \Delta{t} F_2(u',\sigma) + \end{aligned} + +and then this can be written as + +.. math:: q^{n+1}_{adv} = \frac{\bar{N}+N'}{\bar{D}+D'} + +Noting that + +.. math:: + + \label{eq:adv_ls} + \bar{q}^{n+1}_{adv} = \frac{\bar{N}}{\bar{D}} + +then the advective perturbation field can be found from +:math:`q^{n+1}_{adv} - \bar{q}^{n+1}_{adv}`, i.e. + +.. math:: q'^{n+1}_{adv} = \frac{\bar{N}+N'}{\bar{D}+D'} - \frac{\bar{N}}{\bar{D}} + +Using the first terms of the expansion +:math:`(1+D'/\bar{D})^{-1} \approx 1 - D'/\bar{D}` we get + +.. math:: + + \begin{aligned} + q'^{n+1}_{adv}=& \frac{\bar{N}+N'}{\bar{D}}(1-D'/\bar{D}) - \frac{\bar{N}}{\bar{D}} \\ + =& \frac{N'}{\bar{D}} - \frac{\bar{N}D'}{\bar{D}\bar{D}} - \frac{N'D'}{\bar{D}\bar{D}} + \end{aligned} + +Using `[eq:adv_ls] <#eq:adv_ls>`__, and the fact products of prime +variables equals zero, this simplifies to + +.. math:: q'^{n+1}_{adv} = \frac{N' - D' \bar{q}^{n+1}_{adv}}{\bar{D}} + +Horizontal Splitting +----------------------- + +The horizontal FFSL scheme uses COSMIC/Lin-Rood splitting which has, for +the non-linear scheme, the form + +.. math:: + + \begin{aligned} + q^x &= q^{n} - \frac{\Delta}{2}F^x_{adv}(u,q^n) \\ + q^x &= q^{n} - \frac{\Delta}{2}F^y_{adv}(v,q^n) \\ + q^{n+1} &= q^{n} - \Delta{t} F^x(u,q^y)- \Delta{t} F^y(v,q^x) + \end{aligned} + +Transport Efficiency +------------------------ + +To improve efficiency, we can make the approximation that +:math:`\bar{q}^{n+1}_{adv} \approx \bar{q}^{n}`. This removes the need +to recompute :math:`\bar{q}^{n+1}_{adv}` at every step. From 8d29988826fbda40a0d74a9703f8ea460f346cbf Mon Sep 17 00:00:00 2001 From: Christine Johnson <74597224+cjohnson-pi@users.noreply.github.com> Date: Fri, 22 May 2026 12:30:57 +0100 Subject: [PATCH 2/2] Tidied --- .../tangent_linear_section/index.rst | 6 +- .../tangent_linear_coding.rst | 10 +- .../tangent_linear_components.rst | 66 +--------- ....rst => tangent_linear_ffsl_transport.rst} | 12 +- .../tangent_linear_lin_state.rst | 16 --- .../tangent_linear_mol_transport.rst | 84 ++++++++++++ .../tangent_linear_overview.rst | 4 +- .../tangent_linear_timestepping.rst | 120 +++++------------- 8 files changed, 130 insertions(+), 188 deletions(-) rename documentation/source/science_guide/tangent_linear_section/{tl_ffsl_transport.rst => tangent_linear_ffsl_transport.rst} (92%) delete mode 100644 documentation/source/science_guide/tangent_linear_section/tangent_linear_lin_state.rst create mode 100644 documentation/source/science_guide/tangent_linear_section/tangent_linear_mol_transport.rst diff --git a/documentation/source/science_guide/tangent_linear_section/index.rst b/documentation/source/science_guide/tangent_linear_section/index.rst index 7b63c984c..94bbf6bcd 100644 --- a/documentation/source/science_guide/tangent_linear_section/index.rst +++ b/documentation/source/science_guide/tangent_linear_section/index.rst @@ -9,7 +9,7 @@ Linear Model Section ======================= -This describes the details of the tangent-linear model code for GungHo. +This describes the scientific aspects of the tangent-linear model code for GungHo. .. toctree:: :maxdepth: 2 @@ -18,6 +18,6 @@ This describes the details of the tangent-linear model code for GungHo. tangent_linear_coding tangent_linear_timestepping tangent_linear_components - tl_ffsl_transport - tangent_linear_lin_state + tangent_linear_mol_transport + tangent_linear_ffsl_transport tangent_linear_tests diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_coding.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_coding.rst index 5ef2c5016..0b81cc5ed 100644 --- a/documentation/source/science_guide/tangent_linear_section/tangent_linear_coding.rst +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_coding.rst @@ -1,10 +1,14 @@ .. _tangent_linear_coding: -Tangent Linear Coding +General principles ===================== -To obtain the tangent linear code, we need to differentiate the -nonlinear model code. +To obtain the tangent linear model, we need to differentiate the +nonlinear model. There are two main ways to achieve this: differentiate the nonlinear continuous equations and then discretize, or differentiate the nonlinear discretized equations. + +Mainly the approach of differentiating the discretized equations is used here - the so-called line-by-line approach, as each line of code can be treated individually. However, some parts such as the semi-implicit solve take the approach of differentiating the nonlinear equations and then discretizing (although in fact as the solver is already linear, in practice this just means reusing the nonlinear code). + +The following examples demonstrate how the line-by-line approach is implemented, to differentiate the discretized equations. First derivative ---------------- diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_components.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_components.rst index b9618fb8b..0d2c924a7 100644 --- a/documentation/source/science_guide/tangent_linear_section/tangent_linear_components.rst +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_components.rst @@ -1,6 +1,6 @@ .. _tangent_linear_components: -Components +Right-hand-side Components ===================== Kinetic energy gradient @@ -163,67 +163,3 @@ The linear model is \delta N & = \frac{\partial N}{\partial E} \left( \frac{ \partial E}{\partial \Pi} \delta \Pi + \frac{ \partial E}{\partial \rho} \delta \rho + \frac{ \partial E}{\partial \theta} \delta \theta \right) \\ & = - \frac{1}{E} \left( \frac{1-\kappa}{\kappa} \frac{\delta \Pi}{\Pi} - \frac{\delta \rho}{\rho} - \frac{\delta \theta}{\theta} \right) \end{aligned} - -Method of Lines (MoL) and Runge-Kutta Transport ------------ - -The nonlinear model is: - -.. math:: - - \begin{aligned} - u^{(i)} & = u^n + \Delta t \sum_{k=0}^{i-1} \beta_{i,k} f(u^{(k)}) \quad i=1,\ldots,m \\ - u^{n+1} & = u^{(m)} - \end{aligned} - -In ``rk_alg_timestep``, the :math:`f(u^{(i)})` corresponds to -``rhs_prediction``, and :math:`\sum_{k=0}^{i-1} \beta_{i,k} f(u^{(k)})` -corresponds to rhs. (There is also an extra mass matrix inversion, so -they are not directly equal). - -The tangent linear model comprises two steps. First the linearisation -state is computed: - -.. math:: \bar{u}^{(i)} = \bar{u}^n + \Delta t \sum_{k=0}^{i-1} \beta_{i,k} f(\bar{u}^{(k)}) \quad i=1,\ldots,m - -then the perturbations are computed: - -.. math:: - - \begin{aligned} - u'^{(i)} & = u'^n + \Delta t \sum_{k=0}^{i-1} \beta_{i,k} f'(u'^{(k)},\bar{u}^{(k)}) \quad i=1,\ldots,m \\ - u'^{n+1} & = u'^{(m)} - \end{aligned} - -In the nonlinear model, the method of lines transport code has switches -to ensure that the transport scheme is upwinding. In the linear model, -we use the same switches as the nonlinear model. This means that the -direction is based on the linearisation state wind rather than the -perturbation wind. - -In the nonlinear model, - -:: - - direction = wind(map_w2(df) + k )*v_dot_n(df) - if ( direction > 0.0_r_def ) then - -In the linear model, - -:: - - direction = ls_wind(map_w2(df) + k )*v_dot_n(df) - if ( direction > 0.0_r_def ) then - -Similarly, where - -:: - - sign_X(unit_wind,1.0_r_def, advecting_wind) - -the corresponding tangent linear uses - -:: - - sign_X(unit_wind, 1.0_r_def, ls_advecting_wind) - diff --git a/documentation/source/science_guide/tangent_linear_section/tl_ffsl_transport.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_ffsl_transport.rst similarity index 92% rename from documentation/source/science_guide/tangent_linear_section/tl_ffsl_transport.rst rename to documentation/source/science_guide/tangent_linear_section/tangent_linear_ffsl_transport.rst index c3b6ec7ee..2589af727 100644 --- a/documentation/source/science_guide/tangent_linear_section/tl_ffsl_transport.rst +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_ffsl_transport.rst @@ -22,7 +22,6 @@ the form .. math:: - \label{eq:advective} q^{n+1}_{adv} = \frac{q^{n} - \Delta{t} F(u,q^n)}{\sigma - \Delta{t} F(u,\sigma)} where :math:`\sigma=1` is the unity field, and therefore @@ -54,15 +53,12 @@ Leaving the conservative perturbation transport equation as .. math:: - \label{eq:transport_equation} \frac{\partial q'}{\partial t} + \nabla \cdot( \bar{u}q' ) + \nabla \cdot( u'\bar{q} ) = 0 Flux Computation ---------------------- -The flux computation can be split into two parts, corresponding to the -two terms in equation -`[eq:transport_equation] <#eq:transport_equation>`__. +The flux computation can be split into two parts. The first term, corresponding to :math:`\bar{u}q'`, can be computed using standard FFSL operators, provided the scheme is linear (i.e. no @@ -91,7 +87,7 @@ We can now use the notation Advective Update -------------------- -Following `[eq:advective] <#eq:advective>`__, if we expand the variables +If we expand the variables into linearisation state and perturbation parts we find the flux divergences can be written as @@ -105,7 +101,6 @@ and the advective form becomes .. math:: - \label{eq:advective_linear} q^{n+1}_{adv} = \frac{\bar{q} + q' - \Delta{t} F(\bar{u},\bar{q}) - \Delta{t}F_1(\bar{u},q') - \Delta{t}F_2(u',\bar{q}) }{\sigma - \Delta{t} F(\bar{u},\sigma)- \Delta{t} F_2(u',\sigma)} Let @@ -127,7 +122,6 @@ Noting that .. math:: - \label{eq:adv_ls} \bar{q}^{n+1}_{adv} = \frac{\bar{N}}{\bar{D}} then the advective perturbation field can be found from @@ -145,7 +139,7 @@ Using the first terms of the expansion =& \frac{N'}{\bar{D}} - \frac{\bar{N}D'}{\bar{D}\bar{D}} - \frac{N'D'}{\bar{D}\bar{D}} \end{aligned} -Using `[eq:adv_ls] <#eq:adv_ls>`__, and the fact products of prime +Using the fact products of prime variables equals zero, this simplifies to .. math:: q'^{n+1}_{adv} = \frac{N' - D' \bar{q}^{n+1}_{adv}}{\bar{D}} diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_lin_state.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_lin_state.rst deleted file mode 100644 index 514478d2c..000000000 --- a/documentation/source/science_guide/tangent_linear_section/tangent_linear_lin_state.rst +++ /dev/null @@ -1,16 +0,0 @@ -.. _tangent_linear_lin_state: - -Linearisation state -=================== - -If the linear model is run using idealized initial conditions, then it -is possible to define an analytical linearisation state. Otherwise, the -linearisation state needs to be read and updated from file. The -linearisation state is continually updated over the linear model run -because it is not fixed in time; the linearisation state is the -nonlinear model trajectory. - -The linearisation state is read and updated from file in the same way as -the ancillaries, and the lateral-boundary-conditions. This means that -the data is read with XIOS using a time-axis and on every timestep, the -data is updated and interpolated in time. diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_mol_transport.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_mol_transport.rst new file mode 100644 index 000000000..0a24fe56a --- /dev/null +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_mol_transport.rst @@ -0,0 +1,84 @@ +.. _tangent_linear_mol_transport: + +MoL Transport +==================== + +The Method of Lines (MoL) transport scheme is a Finite Volume Scheme that uses Runge-Kutta timestepping. For the linear model, we use the notation that :math:`\bar{u}` is the linearisation state and :math:`u'` is the perturbation. + +Nonlinear MoL +---------- + +The nonlinear continuous equation is + +.. math:: \frac{\partial \theta}{\partial t} = - u. \nabla \theta + +which is discretised using an upwind scheme. On each face, + +.. math:: \begin{aligned} \frac{\partial \theta}{\partial t} + & = - u. \nabla_L \theta & \text{ if } u>0 \\ + & = - u. \nabla_R \theta & \text{ if } u<0 \\ + & = 0 & \text{ if } u=0 \\ + \end{aligned} + +where :math:`\nabla_L` uses the upwind reconstruction calculated from the left and :math:`\nabla_R` uses the upwind reconstruction calculated from the right. +As required, :math:`\frac{\partial \theta}{\partial t}` tends to 0 as u tends to 0. + +This is implemented by initialising +:math:`\frac{\partial \theta}{\partial t}` as :math:`0` and then adding +on an increment if :math:`u>0` or :math:`u<0`. + +Linearised MoL +---------- + +The linearised continuous equation is + +.. math:: \frac{\partial \theta'}{\partial t} = - \bar{u}. \nabla \theta' - u'. \nabla \bar{\theta} + +This is discretised as + +.. math:: \begin{aligned} \frac{\partial \theta'}{\partial t} + & = - \bar{u}. \nabla_L \theta' - u'. \frac{1}{2} (\nabla_L + \nabla_R) \bar{\theta} + & \text{ if } \bar{u}>0 \\ + & = - \bar{u}. \nabla_R \theta' - u'. \frac{1}{2} (\nabla_L + \nabla_R) \bar{\theta} + & \text{ if } \bar{u}<0 \\ + & = - u'. \frac{1}{2} (\nabla_L + \nabla_R) \bar{\theta} + & \text{ if } \bar{u}=0 \\ + \end{aligned} + +which, as required, tends to +:math:`- u'. \frac{1}{2} (\nabla_L + \nabla_R) \bar{\theta}` as +:math:`\bar{u}` tends to :math:`0`. + +This is implemented by initialising +:math:`\frac{\partial \theta'}{\partial t}` as +:math:`- u'. \frac{1}{2} (\nabla_L + \nabla_R) \bar{\theta}` and then +adding an increment if :math:`\bar{u}>0` or :math:`\bar{u}<0`. + +Runge-Kutta +----------- + +The nonlinear Runge Kutta scheme is: + +.. math:: + + \begin{aligned} + u^{(i)} & = u^n + \Delta t \sum_{k=0}^{i-1} \beta_{i,k} f(u^{(k)}) \quad i=1,\ldots,m \\ + u^{n+1} & = u^{(m)} + \end{aligned} + +The tangent linear Runge Kutta scheme comprises two steps. First the linearisation +state is computed: + +.. math:: \bar{u}^{(i)} = \bar{u}^n + \Delta t \sum_{k=0}^{i-1} \beta_{i,k} f(\bar{u}^{(k)}) \quad i=1,\ldots,m + +then the perturbations are computed: + +.. math:: + + \begin{aligned} + u'^{(i)} & = u'^n + \Delta t \sum_{k=0}^{i-1} \beta_{i,k} f'(u'^{(k)},\bar{u}^{(k)}) \quad i=1,\ldots,m \\ + u'^{n+1} & = u'^{(m)} + \end{aligned} + + +If ``transport_efficiency`` is set to true, then the ``ls_state`` is not updated, so that the same linearisation state is used for all iterations. diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_overview.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_overview.rst index 67e308ae2..c3f631af9 100644 --- a/documentation/source/science_guide/tangent_linear_section/tangent_linear_overview.rst +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_overview.rst @@ -1,6 +1,6 @@ .. _tangent_linear_overview: -Tangent Linear Model Overview +Overview ======================= In variational data assimilation, the aim is to compute the optimal @@ -17,7 +17,7 @@ is the observation operator (which we can assume is linear here), and :math:`\mathbf{B}` and :math:`\mathbf{R}` are the error covariance matrices. -As :math:`N` is nonlinear, there is no unique solution. Therefore, the +As :math:`N` is nonlinear, there may not be a unique solution and it is difficult to solve. Therefore, the problem is rewritten in incremental form, by rewriting the initial state in terms of a linearisation state plus a perturbation, :math:`\mathbf{x}_0 = \tilde{\mathbf{x}} + \mathbf{\delta x}`, so that diff --git a/documentation/source/science_guide/tangent_linear_section/tangent_linear_timestepping.rst b/documentation/source/science_guide/tangent_linear_section/tangent_linear_timestepping.rst index 98f2c5cbe..574002812 100644 --- a/documentation/source/science_guide/tangent_linear_section/tangent_linear_timestepping.rst +++ b/documentation/source/science_guide/tangent_linear_section/tangent_linear_timestepping.rst @@ -6,12 +6,12 @@ Iterated (Semi-) Implicit Timestepping Nonlinear Equations ------------------- -The nonlinear equations are +The continous nonlinear equations are .. math:: \begin{aligned} - \frac{\partial u}{\partial t} &= - u . \nabla u - 2 \Omega \times u - \nabla \Phi - c_P \theta \nabla \Pi \\ + \frac{\partial u}{\partial t} &= - u . \nabla u - 2 \Omega \times u - \nabla \Phi - c_p \theta \nabla \Pi \\ \frac{\partial \rho}{\partial t} & = - \nabla . (\rho u) \\ \frac{\partial \theta}{\partial t} & = - u. \nabla \theta \end{aligned} @@ -20,27 +20,10 @@ together with the nonlinear equation of state: .. math:: \Pi ^{\left( \frac{1-\kappa}{\kappa} \right)} = \frac{R}{p_0} \rho \theta -Linear Equations ----------------- - -The linearized equations are: - -.. math:: - - \begin{aligned} - \frac{\partial \delta u}{\partial t} &= - \bar{u} . \nabla \delta u - \delta {u} . \nabla \bar{ u} - 2 \Omega \times \delta u - c_P \delta \theta \nabla \bar{\Pi} - c_P \bar{\theta} \nabla \delta \Pi \\ - \frac{\partial \delta \rho}{\partial t} & = - \nabla . (\bar{\rho} \delta u) - \nabla . (\delta \rho \bar{u}) \\ - \frac{\partial \delta \theta}{\partial t} & = - \bar{u}. \nabla \delta \theta - \delta u. \nabla \bar{\theta} - \end{aligned} - -together with the linear equation of state: - -.. math:: \frac{ \delta \Pi}{\bar{\Pi}} = \frac{\kappa}{1-\kappa} \left( \frac{\delta \rho}{\bar{\rho}} + \frac{\delta \theta}{\bar{\theta}} \right) - - -Iterated (Semi-) implicit Timestepping --------------------------- +where the prognostic variables are :math:`u` the velocity vector, :math:`\rho` the density, :math:`\theta` the potential temperature, and :math:`\Pi` the Exner pressure. The constants are :math:`\Omega` the rotation vector, :math:`\Phi` the geopotential and :math:`c_p` the specific heat capacity. +Discretisation +-------------------- In the iterated-implicit timestepping, the equations are solved iteratively using a Newton method. If the full nonlinear discrete equations are written as @@ -82,107 +65,64 @@ where .. math:: E = \frac{p_0}{R} \frac{\Pi^{*(1-\kappa)/\kappa}}{\theta^* \rho^*} -We can proceed in a similar way with the linear model, so that - -.. math:: - - \begin{aligned} - \delta \mathbf{x}^{n+1}_{(k+1)} & = \delta \mathbf{x}^{n+1}_{(k)} + \delta \mathbf{x}'_{(k)} \\ - \mathcal{L} (\mathbf{x}^*) \delta \mathbf{x}'_{(k)} & = - R ( \mathbf{x}_{(k)}) \delta \mathbf{x}_{(k)} - \end{aligned} - -where :math:`R` is the linearized version of :math:`\mathcal{R}`. - -As the operator :math:`\mathcal{L}` is already linear, this does not -need to be linearized. - -Discretisation --------------------- - -The nonlinear model evolves the linearisation state -:math:`\bar{\mathbf{x}}`. The aim is to find the state -:math:`\bar{\mathbf{x}}^n` at time :math:`t_{n+1}`. This is solved -iteratively, by letting -:math:`{\bar{\mathbf{x}}}^{n+1}_{(k+1)} = \bar{\mathbf{x}}^{n+1}_{(k)} + \bar{\mathbf{x}}'_{(k+1)}` -where :math:`k` is the iteration count. - -The increment :math:`\bar{\mathbf{x}}'_{(k+1)}` solves the linear -equation - -.. math:: \mathcal{L}( \bar{\mathbf{x}}*) \bar{\mathbf{x}}'_{(k+1)} = \mathcal{R}^n(\bar{\mathbf{x}}^n) + \mathcal{R}^{n+1}(\bar{\mathbf{x}}^{n+1}_{(k)}) + \mathcal{R}^{adv}(\bar{\mathbf{x}}^n, \bar{\mathbf{x}}^{n+1}_{(k)}) - -where :math:`\mathcal{L}( \bar{\mathbf{x}}*)` is a linear operator based -on the reference state :math:`\bar{\mathbf{x}}*`. - The procedure is: -- | Let :math:`\bar{\mathbf{x}}^{n+1}_{(0)} = \bar{\mathbf{x}}^n` +- | Let :math:`\mathbf{x}^{n+1}_{(0)} = \mathbf{x}^n` | The first iteration is given by the state at the previous timestep. - | Set - :math:`\bar{\mathbf{x}}* \equiv (\bar{\rho}*, \bar{\theta}*, \bar{\Pi}* ) = (\bar{\rho}^n, \bar{\theta}^n, \bar{\Pi}^n)` + :math:`\mathbf{x}* \equiv (\rho*, \theta*, \Pi* ) = (\rho^n, \theta^n,\Pi^n)` | The reference state is given by the state at the previous timestep (for density, potential temperature and exner pressure only). -- | Compute :math:`\mathcal{R}^n(\bar{\mathbf{x}}^n)` - | This is the RHS based on the previous timestep. - - Loop k=1, K - - | Compute :math:`\mathcal{R}^{n+1}(\bar{\mathbf{x}}^{n+1}_{(k)})` - and - :math:`\mathcal{R}^{adv}(\bar{\mathbf{x}}^n, \bar{\mathbf{x}}^{n+1}_{(k)})` - | This is the RHS and advective terms based on the current timestep. + - | Compute :math:`\mathcal{R}(\mathbf{x}^n, \mathbf{x}^{n+1}_{(k)})` - - Solve the linear equation for :math:`\bar{\mathbf{x}}'_{(k+1)}` and - hence for :math:`\bar{\mathbf{x}}^{n+1}_{(k+1)}`. + - Solve the linear equation for :math:`\mathbf{x}'_{(k+1)}` and + hence for :math:`\mathbf{x}^{n+1}_{(k+1)}`. -The linear model evolves the perturbations :math:`\mathbf{x}`. To find -the state at time :math:`t_{n+1}`, we follow the same procedure as for -the nonlinear model, but using the tangent linear versions of the right -hand side in the equation. +Linear Equations +----------------------- + +We can proceed in a similar way with the linear model. Using the notation of :math:`\delta \mathbf{x}` is the perturbation and :math:`\mathbf{x}` is the linearisation state, -Let -:math:`{\mathbf{x}}^{n+1}_{(k+1)} = \mathbf{x}^{n+1}_{(k)} + \mathbf{x}'_{(k+1)}` -where :math:`\mathbf{x}'_{(k+1)}` solves +.. math:: -.. math:: \mathcal{L} ( \bar{\mathbf{x}}*) \mathbf{x}'_{(k+1)} = R^n(\mathbf{x}^n, \bar{\mathbf{x}}^n) + R^{n+1}(\mathbf{x}^{n+1}_{(k)},\bar{\mathbf{x}}^{n+1}_{(k)}) + R^{adv}(\mathbf{x}^n, \mathbf{x}^{n+1}_{(k)},\bar{\mathbf{x}}^n, \bar{\mathbf{x}}^{n+1}_{(k)}) + \begin{aligned} + \delta \mathbf{x}^{n+1}_{(k+1)} & = \delta \mathbf{x}^{n+1}_{(k)} + \delta \mathbf{x}'_{(k)} \\ + \mathcal{L} (\mathbf{x}^*) \delta \mathbf{x}'_{(k)} & = - R ( \mathbf{x}_{(k)}) \delta \mathbf{x}_{(k)} + \end{aligned} -where :math:`R^n` is the tangent linear of :math:`\mathcal{R}^n` etc. -The tangent linear operators are a function of both the perturbation -:math:`\mathbf{x}` and the linearisation state :math:`\bar{\mathbf{x}}`. +where :math:`R` is the linearized version of :math:`\mathcal{R}`. As the operator :math:`\mathcal{L}` is already linear, this does not need to be linearized. The procedure is: -- | Let :math:`\mathbf{x}^{n+1}_{(0)} = \mathbf{x}^n` +- | Let :math:`\delta \mathbf{x}^{n+1}_{(0)} = \delta \mathbf{x}^n` | The first iteration is given by the perturbation at the previous timestep. - | Set - :math:`\bar{\mathbf{x}}* \equiv (\bar{\rho}*, \bar{\theta}*, \bar{\Pi}* ) = (\bar{\rho}^n, \bar{\theta}^n, \bar{\Pi}^n)` + :math:`\mathbf{x}* \equiv (\rho*,\theta*, \Pi* ) = (\rho^n, \theta^n, \Pi^n)` | The reference state is given by the linearisation state at the previous timestep - so that the linear operator is identical to the that in the nonlinear model. -- Compute :math:`R^n(\mathbf{x}^n, \bar{\mathbf{x}}^n)` - - Loop k=1, K - Compute - :math:`R^{n+1}(\mathbf{x}^{n+1}_{(k)},\bar{\mathbf{x}}^{n+1}_{(k)})` - and - :math:`R^{adv}(\mathbf{x}^n, \mathbf{x}^{n+1}_{(k)},\bar{\mathbf{x}}^n, \bar{\mathbf{x}}^{n+1}_{(k)})` + :math:`R ( \mathbf{x}_{(k)}) \delta \mathbf{x}_{(k)}` - - Solve for :math:`\mathbf{x}'_{(k+1)}` and hence for - :math:`\mathbf{x}^{n+1}_{(k+1)}`. + - Solve for :math:`\delta \mathbf{x}'_{(k)}` using :math:`\mathcal{L} (\mathbf{x}^*) \delta \mathbf{x}'_{(k)} = R ( \mathbf{x}_{(k)}) \delta \mathbf{x}_{(k)}` and hence for + :math:`\delta \mathbf{x}^{n+1}_{(k+1)}`. -The linear solver :math:`\mathcal{L} ( \bar{\mathbf{x}}*)` is exactly -the same as for the nonlinear model. +The linear operator :math:`\mathcal{L} ( \bar{\mathbf{x}}*)` is exactly +the same as for the nonlinear model. The solution of the equation (often known as the solver) uses the same method as the nonlinear model. This could be multi-grid or Krylov iterative methods, or a combination of the both. A very exact line-by=line procedure would linearise every iteration of the solver procedure, using the nonlinear solution at every step. However, this is not necessary and the more linearise then discretize approach is taken, to form the linear equation and then solve. -The right hand side terms :math:`R^n`, :math:`R^{n+1}` and -:math:`R^{adv}` require the linearisation state at the same iteration. -Therefore, for an exact tangent linear semi-implicit algorithm, the +The right hand side term :math:`R ( \mathbf{x}_{(k)}) \delta \mathbf{x}_{(k)}` require both the perturbation :math:`\delta \mathbf{x}_{(k)}` and the linearisation state :math:`\mathbf{x}_{(k)}` at the same iteration. Therefore, for an exact tangent linear semi-implicit algorithm, the nonlinear semi-implicit step needs to be run first, and the intermediate ``ls_state`` values need to be stored and used in the following tangent linear semi-implicit step. +In practice the impact of the variation of the linearisation state through the timestep is very small. Furthermore, the evolution of the linearisation state in the linear model timestepping does not include any physics and will be at a different resolution to the outer loop nonlinear model, and hence will not be very accurate. Therefore, it is advised for both efficiency and accuracy reasons to use the ``fixed_ls`` option. This option applies the ``ls_state`` values set at the beginning of the timestep to all iterations, without recalculation. +