# [Greenline] Warm Bubble: Least-Squares Coefficients in Tracer Advection

<!-- Add the tag for the current cycle number in the top bar -->

- Shaped by: @ChiaRuiOng , @nfarabullini
- Appetite (FTEs, weeks):
- Developers: Nikki, Rico (for later)

## Problem
<!-- The raw idea, a use case, or something we’ve seen that motivates us to work on this -->
The horizontal advection granule adopts second-order Miura advection scheme. However, the least-squares coefficients required in second-order Miura advection scheme with nonlinear local reconstruction are not ported to icon4py. A local reconstruction here means a 2D polynommial interpolation of a local patch from values at cell centers bounded by the C2E2C polygon. In mathematical notation, the interpolated value at coordinates (x,y) inside the local patch can be found by `q = Ba`, where `q` is the tracer value at neighboring cell centers and `a` is the coordinates. `A` is the least-squares matrix which is a function of distance between neighboring cell centers and the origin. We invert (pseudo-inverse if number of rows > number of unknowns) the matrix and store as pre-computed coefficients. This inverted matrix is used when computing interpolated values `f = B-1 q`,

B = (
     x1, y1, x1^2, ..., 
     x2, y2, x2^2, ..., 
     ...,
     xn, yn, xn^2, ..., 
    )
a = (a_x, a_y, a_xx, ...)
q = (q1, q2, q3, ...)

We only use the linear polynomial for the reconstruction in icon4py because only the second-order Miura scheme is currently supported.


## Appetite

<!-- Explain how much time we want to spend and how that constrains the solution -->

## Solution
<!-- The core elements we came up with, presented in a form that’s easy for people to immediately understand -->

- The least-squares coefficients are declared as `lsq_pseudoinv` in ICON.
- `dim_c`, `dim_unk`, and `wgt_exp` are namelist parameters. They are set in [mo_interpol_config.f90](https://github.com/C2SM/icon-exclaim/blob/7076be7d4d0443d9027f9e2871d530da47af31dc/src/configure_model/mo_interpol_config.f90#L351).
- The actual computation of the least-squares coefficients can be found in subroutines [lsq_compute_coeff_cell](https://github.com/C2SM/icon-exclaim/blob/b392c1a77a331015a753bc9386eea4339e994e1a/src/shr_horizontal/mo_intp_coeffs_lsq_bln.f90#L536C12-L536C34) in `mo_intp_coeffs_lsq_bln.f90`.
- The code seems to be **slightly different in Torus and icosahedral grids**. We need to port most of the computation in [lsq_compute_coeff_cell_sphere](https://github.com/C2SM/icon-exclaim/blob/b392c1a77a331015a753bc9386eea4339e994e1a/src/shr_horizontal/mo_intp_coeffs_lsq_bln.f90#L595) and [lsq_compute_coeff_cell_torus](https://github.com/C2SM/icon-exclaim/blob/b392c1a77a331015a753bc9386eea4339e994e1a/src/shr_horizontal/mo_intp_coeffs_lsq_bln.f90#L1356).
- The intermediate variables `lsq_dim_stencil`, `lsq_idx_c`, `lsq_blk_c`, `lsq_weights_c`, `lsq_qtmat_c`, `lsq_moments`, `lsq_moments_hat`, `lsq_rmat_utri_c`, `lsq_rmat_rdiag_c` do not need to be stored as interpolation fields in icon4py.
- We may also need to port [create_stencil_c3](https://github.com/C2SM/icon-exclaim/blob/b392c1a77a331015a753bc9386eea4339e994e1a/src/shr_horizontal/mo_intp_coeffs_lsq_bln.f90#L100), because it is called to compute the `C2E2C` neighbor table which is subsequently used in the construction of `lsq_pseudoinv`. **Check whether we already have the neighbor table in icon4py**.
- Namelist parameter `llsq_lin_consv` is always `False`, which is the default value, in all EXCLAIM run scripts.
- we set `llsq_svd` to `True` in our exclaim run scripts. The part of the code if `llsq_svd == False` does not need to be ported.

In summary, the steps are:
1. merge the copy-paste and do `match geometry ICO/TORUS`
2. construct neighbor table `C2E2C` (may not be necessary).
3. construct the least-squares matrix `B`.
4. invert the matrix by utilizing Python library (or simply numpy?).
5. store the inverted matrix as `lsq_pseudoinv_x` and `lsq_pseudoinv_y` in the interpolation factory.
6. construct a data test for the least-squares coefficients.

## Rabbit holes
<!-- Details about the solution worth calling out to avoid problems -->
Higher-order least squares interpolation is not required for the warm bubble experiment.

## No-gos
<!-- Anything specifically excluded from the concept: functionality or use cases we intentionally aren’t covering to fit the ## appetite or make the problem tractable -->

## Progress

- [x] merge the copy-paste and do `match geometry ICO/TORUS`
- [x] construct the least-squares matrix `B`.
- [x] invert the matrix by utilizing library
- [x] construct a data test for the least-squares coefficients.
- [ ] store the inverted matrix as `lsq_pseudoinv_x` and `lsq_pseudoinv_y` / make `lsq_pseudoinv_1` and `lsq_pseudoinv_2` factories
- [ ] merge functionality of `math/projection.py::plane_torus_closest_coordinates` with `math/helpers.py::diff_on_edges_torus`
- [ ] merge `math/projection.py::gnomonic_proj` and `gnomonic_proj_single_val`
- [ ] cleanup gridShape lat/lon + cart_x/y interface to compute_coeffs