Non-orthogonal case for the gradient in the advection field RTheta

simulations/geometryRTheta/diocotron/CMakeLists.txt

@@ -31,7 +31,7 @@ function (diocotron_executable PREDCORR_METHOD TIME_METHOD )
             gslx::geometry_RTheta
             gslx::interpolation_2D_rp
             gslx::io
-            gslx::advection_rp
+            gslx::advection_polar
             gslx::time_integration_rp
             gslx::quadrature
             gslx::utils

simulations/geometryRTheta/diocotron/diocotron.cpp

@@ -10,7 +10,7 @@
 #include <paraconf.h>
 #include <pdi.h>
 
-#include "bsl_advection_rp.hpp"
+#include "bsl_advection_rtheta.hpp"
 #include "bsl_predcorr.hpp"
 #include "bsl_predcorr_second_order_explicit.hpp"
 #include "bsl_predcorr_second_order_implicit.hpp"

simulations/geometryRTheta/vortex_merger/CMakeLists.txt

@@ -13,7 +13,7 @@ target_link_libraries(vortex_merger
         Eigen3::Eigen
 
 
-        gslx::advection_rp
+        gslx::advection_polar
         gslx::geometry_RTheta
         gslx::io
         gslx::interpolation_2D_rp

simulations/geometryRTheta/vortex_merger/vortex_merger.cpp

@@ -11,7 +11,7 @@
 #include <paraconf.h>
 #include <pdi.h>
 
-#include "bsl_advection_rp.hpp"
+#include "bsl_advection_rtheta.hpp"
 #include "bsl_predcorr.hpp"
 #include "bsl_predcorr_second_order_explicit.hpp"
 #include "bsl_predcorr_second_order_implicit.hpp"

src/geometryRTheta/advection/CMakeLists.txt

@@ -1,11 +1,11 @@
 # SPDX-License-Identifier: MIT
 
-add_library(advection_rp INTERFACE)
-target_include_directories(advection_rp
+add_library(advection_polar INTERFACE)
+target_include_directories(advection_polar
     INTERFACE
         "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>"
 )
-target_link_libraries(advection_rp 
+target_link_libraries(advection_polar 
     INTERFACE
         DDC::core
         Eigen3::Eigen
@@ -17,4 +17,4 @@ target_link_libraries(advection_rp
         gslx::math_tools
         gslx::utils
 )
-add_library(gslx::advection_rp ALIAS advection_rp)
+add_library(gslx::advection_polar ALIAS advection_polar)

src/geometryRTheta/advection/README.md

@@ -28,10 +28,10 @@ Y(s; s, x, y) = y.
 ```
 
 So to compute the advected function at the next time step, 
- - we compute the characteristic feet $`X(t^n; t^{n+1}, x_i, y_j)`$ and  $`Y(t^n; t^{n+1}, x_i, y_j)`$ 
+ - we compute the characteristics' feet $`X(t^n; t^{n+1}, x_i, y_j)`$ and $`Y(t^n; t^{n+1}, x_i, y_j)`$ 
  for each mesh points $`(x_i, y_j)`$ with a time integration method ; 
- - we interpolate the function $f(t = t^n)$ on the characteristic feet. 
- The proprety ensures that the interpolation gives the function at the next time step $f(t = t^{n+1})$.
+ - we interpolate the function $f(t = t^n)$ on the characteristics' feet. 
+ The property ensures that the interpolation gives the function at the next time step $f(t = t^{n+1})$.
 
 
 
@@ -47,7 +47,7 @@ There are multiple time integration methods available which are implemented in t
 
 We are listing the different schemes for this equation $`\partial_t X (t) = A_x(t, X(t),Y(t))`$. 
 
-We write $X (t) = X (t; s, x, y)$,  $X^n = X(t^n)$ and $A^n(X) = A(t^n, X)$ for a time discretisation $`\{t^n; t^n > t^{n-1},  \forall n\}_n`$. 
+We write $X (t) = X (t; s, x, y)$,  $X^n = X(t^n)$ and $A^n(X) = A(t^n, X)$ for a time discretization $`\{t^n; t^n > t^{n-1},  \forall n\}_n`$. 
 
 
 ### Explicit Euler method
@@ -112,12 +112,12 @@ compute the mesh points in the physical domain using a mapping function $\mathca
 
 This adds some steps to the advection operator, we now have to compute 
  - the mesh points in the physical domain using $\mathcal{F}$; 
- - the characteristic feet in the physical domain; 
- - the characteristic feet in the logical domain (polar grid) using $\mathcal{F}^{-1}$; 
- - then interpolate the advection function at the  characteristic feet in the logical domain. 
+ - the characteristics' feet in the physical domain; 
+ - the characteristics' feet in the logical domain (polar grid) using $\mathcal{F}^{-1}$; 
+ - then interpolate the advection function at the characteristics' feet in the logical domain. 
 
 The third step can be difficult especially if the mapping function $\mathcal{F}$ is not analytically invertible. 
-It is not impossible but the computations can be costly. 
+It is not impossible, but the computations can be costly. 
 
 
 That is why, we introduce the **pseudo-Cartesian domain**. 
@@ -130,12 +130,12 @@ We use another mapping function $\mathcal{G}$ such that:
 Then the four previous steps become
  - calculate the mesh points in the pseudo-Cartesian domain using $\mathcal{G}$; 
  - calculate the advection field $A$ in the pseudo-Cartesian domain using the Jacobian matrix of $(\mathcal{F}\circ\mathcal{G}^{-1})^{-1}$; 
- - calculate the characteristic feet in the pseudo\_Cartesian domain; 
- - calculate the characteristic feet in the logical domain (polar grid) using $\mathcal{G}^{-1}$; 
- - interpolate the advection function at the  characteristic feet in the logical domain. 
+ - calculate the characteristics' feet in the pseudo\_Cartesian domain; 
+ - calculate the characteristics' feet in the logical domain (polar grid) using $\mathcal{G}^{-1}$; 
+ - interpolate the advection function at the characteristics' feet in the logical domain. 
 
 Here, $\mathcal{G}$ is analytically invertible (we can fix  $\mathcal{G}^{-1}(x = 0, y = 0) = (r = 0, \theta = 0)$) 
-and  $`(J_{\mathcal{F}}J_{\mathcal{G}}^{-1})^{-1}`$ is well-defined. The details are given in Edoardo Zoni's article [1]. 
+and  $`(J_{\mathcal{F}}J_{\mathcal{G}}^{-1})^{-1}`$ is well-defined. The details are given in [Zoni et al. (2019)](#zoni). 
 
 
 **Remark 1:** if $\mathcal{F}$ is the circular mapping function, then the physical domain and the pseudo-Cartesian domain are the same. 
@@ -151,14 +151,14 @@ In the studied equation, the advection field is given along the physical domain
 \partial_t f + A_x \partial_x f + A_y \partial_y f = 0.
 ```
 
-The BslAdvectionRP operator can take as input the advection field along the physical domain axis or the advection field along the logical domain axis,
+The BslAdvectionRTheta operator can take as input the advection field along the physical domain axis or the advection field along the logical domain axis,
 ```math
 A = (A_x, A_y) \quad \text{or} \quad A = (A_r, A_\theta).
 ```
 
 The advection field can be computed thanks to the AdvectionFieldFinder operator. This operator returns the advection field along the physical domain axes or the advection field along the logical domain axes (see [advection\_field\_rp](./../advection_field/README.md)).
 
-* If the advection field is directly given along the physical domain axes, no treatment is needed in the BslAdvectionRP operator. 
+* If the advection field is directly given along the physical domain axes, no treatment is needed in the BslAdvectionRTheta operator. 
 
 * If the advection field is given along the logical domain axes, then we need to compute the advection field along the physical domain axes to advect in the physical domain. 
 
@@ -169,20 +169,10 @@ A = - E \wedge e_z = -\nabla \phi \wedge e_z.
 
 In [the documentation for the advection field](./../advection_field/README.md), we show that 
 ```math
-\nabla_{xy} \phi = J D_{G}^{-1}\nabla_{r\theta} \phi,
+\nabla_{xy} \phi = J \nabla_{r\theta} \phi,
 ```
 
-with *J* the Jacobian matrix and the following diagonal matrix 
-```math
-D_{G} 
-= 
-\begin{bmatrix}
-    \sqrt{g_{11}} & 0 \\
-    0 & \sqrt{g_{22}} \\
-\end{bmatrix}
-```
-
-with the metric tensor $`G = J^TJ = [g_{ij}]_{ij}`$. 
+with *J* the Jacobian matrix (and the metric tensor $`G = J^TJ = [g_{ij}]_{ij}`$). 
 
 It gives the following relation for the electric field
 ```math
@@ -191,16 +181,7 @@ It gives the following relation for the electric field
     E_y \\
 \end{bmatrix} 
 = 
-J D_{G}^{-1}
-\begin{bmatrix}
-    E_r \\
-    E_\theta \\
-\end{bmatrix}
-=
-\begin{bmatrix}
-    \frac{j_{11}}{\sqrt{g_{11}}} & \frac{j_{12}}{\sqrt{g_{22}}} \\
-    \frac{j_{21}}{\sqrt{g_{11}}} & \frac{j_{22}}{\sqrt{g_{22}}} \\
-\end{bmatrix}
+J
 \begin{bmatrix}
     E_r \\
     E_\theta \\
@@ -219,51 +200,28 @@ We deduce that
 	- Ey \\
 	E_x 
 \end{bmatrix}
-= 
-\begin{bmatrix}
-	- \frac{j_{21}}{\sqrt{g_{11}}}E_r - \frac{j_{22}}{\sqrt{g_{22}}} E_\theta\\
-	\frac{j_{11}}{\sqrt{g_{11}}}E_r + \frac{j_{12}}{\sqrt{g_{22}}} E_\theta \\
-\end{bmatrix}
-= 
-\begin{bmatrix}
-	\frac{j_{22}}{\sqrt{g_{22}}}  & - \frac{j_{21}}{\sqrt{g_{11}}} \\
-	- \frac{j_{12}}{\sqrt{g_{22}}} & \frac{j_{11}}{\sqrt{g_{11}}}
-\end{bmatrix}
-\begin{bmatrix}
-	- E_\theta \\
-	E_r
-\end{bmatrix}
-= 
-\det(JD_{G}^{-1}) ((JD_{G}^{-1})^{-1})^{T}
+= \det(J)
+J^{-T}
 \begin{bmatrix}
 	- E_\theta \\
 	E_r
 \end{bmatrix}
-= 
-\det(J) (J^{-1})^{T} det(D_{G}^{-1}) D_G 
+= \det(J)
+J^{-T}
 \begin{bmatrix}
 	A_r \\
 	A_\theta
 \end{bmatrix}.
 ```
 
-So, from the advection field along the logical domain axis, we multiply by 
-```math
-\det(J) det(D_{G}^{-1})  (J^{-1})^{T} D_G
-= 
-\begin{bmatrix}
-	\frac{j_{22}}{\sqrt{g_{22}}}  & - \frac{j_{21}}{\sqrt{g_{11}}} \\
-	- \frac{j_{12}}{\sqrt{g_{22}}} & \frac{j_{11}}{\sqrt{g_{11}}}
-\end{bmatrix},
-```
-
+So, from the advection field along the logical domain axis, we multiply by $`J^{-1}`$
 to get the advection field along the physical domain axis. 
 
 
 # Unit tests
 
-The test of the advection operator are implemented in the `tests/geometryRTheta/advection_2d_rp/` folder 
-([advection\_2d\_rp](./../../../tests/geometryRTheta/advection_2d_rp/README.md)).
+The test of the advection operator are implemented in the `tests/geometryRTheta/advection_2d_rtheta/` folder 
+([advection\_2d\_rp](./../../../tests/geometryRTheta/advection_2d_rtheta/README.md)).
 
 
 It tests: 
@@ -280,7 +238,7 @@ It tests:
  	- simulation 2: rotation of Gaussian function 
  		- $`f_0(x,y) = \exp\left( - \frac{(x- x_0)^2}{2 \sigma_x^2} - \frac{(y- y_0)^2}{2 \sigma_y^2} \right)`$, 
  		- $`A(t, x, y) = J_{\mathcal{F}_{\text{circular}}}(v_r, v_\theta)`$. 
- 	- simulation 3: decentred rotation (test given in Edoardo Zoni's article [1]). 
+ 	- simulation 3: decentred rotation (test given in [Zoni et al. (2019)](#zoni)). 
  	 	- $`f_0(x,y) = \frac{1}{2} \left( G(r_1(x,y)) + G(r_2(x,y))\right)`$,
  	 		- with 
  	 			- $`G(r) = \cos\left(\frac{\pi r}{2 a}\right)^4 * 1_{r<a}(r)`$, 
@@ -296,7 +254,7 @@ for $n = 1, 2, 4, 8,  ...$.
 
 
 # References
-[1] Edoardo Zoni, Yaman Güçlü. "Solving hyperbolic-elliptic problems on singular mapped 
+<a name="zoni"></a> [1] Edoardo Zoni, Yaman Güçlü. "Solving hyperbolic-elliptic problems on singular mapped 
 disk-like domains with the method of characteristics and spline finite elements". 
 ([https://doi.org/10.1016/j.jcp.2019.108889](https://doi.org/10.1016/j.jcp.2019.108889).)
 Journal of Computational Physics (2019).
@@ -305,9 +263,8 @@ Journal of Computational Physics (2019).
 # Contents
 
 This folder contains:
- - iadvectionrp.hpp : define the base class for advection operator (IAdvectionRP). 
- 	- bsl\_advection\_rp.hpp : define the advection operator described just before (BslAdvectionRP). 
- - spline\_foot\_finder.hpp : solve the characteristic equation thanks to a given time integration method (IFootFinder).
+ - iadvection\_rtheta.hpp : define the base class for advection operator (IAdvectionRTheta). 
+ 	- bsl\_advection\_rtheta.hpp : define the advection operator described just before (BslAdvectionRTheta). 
 
 
 

src/geometryRTheta/advection/bsl_advection_rp.hpp → src/geometryRTheta/advection/bsl_advection_rtheta.hpp

@@ -6,7 +6,7 @@
 #include "directional_tag.hpp"
 #include "geometry.hpp"
 #include "i_interpolator_2d_rp.hpp"
-#include "iadvectionrp.hpp"
+#include "iadvection_rtheta.hpp"
 #include "metric_tensor.hpp"
 #include "spline_interpolator_2d_rp.hpp"
 #include "spline_polar_foot_finder.hpp"
@@ -145,7 +145,7 @@ class BslAdvectionRTheta : public IAdvectionRTheta
      *
      * @param [in, out] allfdistribu_host
      *      A Field containing the values of the function we want to advect.
-     * @param [in] advection_field_rp
+     * @param [in] advection_field_rtheta
      *      A DConstVectorFieldRTheta containing the values of the advection field
      *      on the logical index range axis.
      * @param [in] advection_field_xy_center
@@ -158,7 +158,7 @@ class BslAdvectionRTheta : public IAdvectionRTheta
      */
     host_t<DFieldRTheta> operator()(
             host_t<DFieldRTheta> allfdistribu_host,
-            host_t<DConstVectorFieldRTheta<R, Theta>> advection_field_rp,
+            host_t<DConstVectorFieldRTheta<R, Theta>> advection_field_rtheta,
             CoordXY const& advection_field_xy_center,
             double dt) const override
     {
@@ -173,24 +173,20 @@ class BslAdvectionRTheta : public IAdvectionRTheta
         // Convert advection field on RTheta to advection field on XY
         host_t<DVectorFieldMemRTheta<X, Y>> advection_field_xy_host(grid);
 
-        MetricTensor<Mapping, CoordRTheta> metric_tensor(m_mapping);
+        InverseJacobianMatrix<Mapping, CoordRTheta> inv_jacobian_matrix(m_mapping);
 
         ddc::for_each(grid_without_Opoint, [&](IdxRTheta const irp) {
             CoordRTheta const coord_rp(ddc::coordinate(irp));
 
-            std::array<std::array<double, 2>, 2> J; // Jacobian matrix
-            m_mapping.jacobian_matrix(coord_rp, J);
-            std::array<std::array<double, 2>, 2> G; // Metric tensor
-            metric_tensor(G, coord_rp);
+            std::array<std::array<double, 2>, 2> inv_J = inv_jacobian_matrix(coord_rp);
+            double const jacobian = m_mapping.jacobian(coord_rp);
 
             ddcHelper::get<X>(advection_field_xy_host)(irp)
-                    = ddcHelper::get<R>(advection_field_rp)(irp) * J[1][1] / std::sqrt(G[1][1])
-                      + ddcHelper::get<Theta>(advection_field_rp)(irp) * -J[1][0]
-                                / std::sqrt(G[0][0]);
+                    = ddcHelper::get<R>(advection_field_rtheta)(irp) * inv_J[0][0] * jacobian
+                      + ddcHelper::get<Theta>(advection_field_rtheta)(irp) * inv_J[1][0] * jacobian;
             ddcHelper::get<Y>(advection_field_xy_host)(irp)
-                    = ddcHelper::get<R>(advection_field_rp)(irp) * -J[0][1] / std::sqrt(G[1][1])
-                      + ddcHelper::get<Theta>(advection_field_rp)(irp) * J[0][0]
-                                / std::sqrt(G[0][0]);
+                    = ddcHelper::get<R>(advection_field_rtheta)(irp) * inv_J[0][1] * jacobian
+                      + ddcHelper::get<Theta>(advection_field_rtheta)(irp) * inv_J[1][1] * jacobian;
         });
 
         ddc::for_each(Opoint_grid, [&](IdxRTheta const irp) {