diff --git a/simulations/geometryRTheta/diocotron/CMakeLists.txt b/simulations/geometryRTheta/diocotron/CMakeLists.txt index 954b22e57..f669e230f 100644 --- a/simulations/geometryRTheta/diocotron/CMakeLists.txt +++ b/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 diff --git a/simulations/geometryRTheta/diocotron/diocotron.cpp b/simulations/geometryRTheta/diocotron/diocotron.cpp index a21696cf2..fa7dfaaa0 100644 --- a/simulations/geometryRTheta/diocotron/diocotron.cpp +++ b/simulations/geometryRTheta/diocotron/diocotron.cpp @@ -10,7 +10,7 @@ #include #include -#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" diff --git a/simulations/geometryRTheta/vortex_merger/CMakeLists.txt b/simulations/geometryRTheta/vortex_merger/CMakeLists.txt index 5f0266d71..dcd1448e7 100644 --- a/simulations/geometryRTheta/vortex_merger/CMakeLists.txt +++ b/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 diff --git a/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp b/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp index c5d512731..fd2ce4de9 100644 --- a/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp +++ b/simulations/geometryRTheta/vortex_merger/vortex_merger.cpp @@ -11,7 +11,7 @@ #include #include -#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" diff --git a/src/geometryRTheta/advection/CMakeLists.txt b/src/geometryRTheta/advection/CMakeLists.txt index 53f80beab..597da08f3 100644 --- a/src/geometryRTheta/advection/CMakeLists.txt +++ b/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 "$" ) -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) diff --git a/src/geometryRTheta/advection/README.md b/src/geometryRTheta/advection/README.md index ccd85a210..95857f747 100644 --- a/src/geometryRTheta/advection/README.md +++ b/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 [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). diff --git a/src/geometryRTheta/advection/bsl_advection_rp.hpp b/src/geometryRTheta/advection/bsl_advection_rtheta.hpp similarity index 90% rename from src/geometryRTheta/advection/bsl_advection_rp.hpp rename to src/geometryRTheta/advection/bsl_advection_rtheta.hpp index 53152572a..4435c29f1 100644 --- a/src/geometryRTheta/advection/bsl_advection_rp.hpp +++ b/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 operator()( host_t allfdistribu_host, - host_t> advection_field_rp, + host_t> 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> advection_field_xy_host(grid); - MetricTensor metric_tensor(m_mapping); + InverseJacobianMatrix inv_jacobian_matrix(m_mapping); ddc::for_each(grid_without_Opoint, [&](IdxRTheta const irp) { CoordRTheta const coord_rp(ddc::coordinate(irp)); - std::array, 2> J; // Jacobian matrix - m_mapping.jacobian_matrix(coord_rp, J); - std::array, 2> G; // Metric tensor - metric_tensor(G, coord_rp); + std::array, 2> inv_J = inv_jacobian_matrix(coord_rp); + double const jacobian = m_mapping.jacobian(coord_rp); ddcHelper::get(advection_field_xy_host)(irp) - = ddcHelper::get(advection_field_rp)(irp) * J[1][1] / std::sqrt(G[1][1]) - + ddcHelper::get(advection_field_rp)(irp) * -J[1][0] - / std::sqrt(G[0][0]); + = ddcHelper::get(advection_field_rtheta)(irp) * inv_J[0][0] * jacobian + + ddcHelper::get(advection_field_rtheta)(irp) * inv_J[1][0] * jacobian; ddcHelper::get(advection_field_xy_host)(irp) - = ddcHelper::get(advection_field_rp)(irp) * -J[0][1] / std::sqrt(G[1][1]) - + ddcHelper::get(advection_field_rp)(irp) * J[0][0] - / std::sqrt(G[0][0]); + = ddcHelper::get(advection_field_rtheta)(irp) * inv_J[0][1] * jacobian + + ddcHelper::get(advection_field_rtheta)(irp) * inv_J[1][1] * jacobian; }); ddc::for_each(Opoint_grid, [&](IdxRTheta const irp) { diff --git a/src/geometryRTheta/advection/iadvectionrp.hpp b/src/geometryRTheta/advection/iadvection_rtheta.hpp similarity index 100% rename from src/geometryRTheta/advection/iadvectionrp.hpp rename to src/geometryRTheta/advection/iadvection_rtheta.hpp diff --git a/src/geometryRTheta/advection_field/README.md b/src/geometryRTheta/advection_field/README.md index acf4801ba..f0e685070 100644 --- a/src/geometryRTheta/advection_field/README.md +++ b/src/geometryRTheta/advection_field/README.md @@ -12,10 +12,14 @@ Currently, the implemented case is: The studied equation system is of the following type : ```math +\left\{ +\begin{aligned} \partial_t \rho + A\cdot\nabla \rho = 0, \\ A = E \wedge e_z, \\ E = - \nabla \phi, \\ - \nabla \cdot \nabla \phi = \rho, +\end{aligned} +\right. ``` with $`\rho`$ the density, $`\phi`$ the electrostatic potential and $`E`$ the electrical field. @@ -34,15 +38,15 @@ If a Field is given as input, it computes the spline representation (on the cros The spline representation is needed to compute the derivatives of the function $`\phi`$. If the PolarSplineMem representation is given as input, it can directly compute the derivatives of the function $`\phi`$. -Once the advection field computed, it is given as input to the BslAdvectionRP operator to advect the density $`\rho`$ function. -The BslAdvectionRP operator can handle the advection with an advection field along $`(x,y)`$ and with an advection field along $`(r,\theta)`$. -But as the BslAdvectionRP operator advects in the physical domain, it is recommend to work with the advection field along $`(x,y)`$. +Once the advection field computed, it is given as input to the BslAdvectionRTheta operator to advect the density $`\rho`$ function. +The BslAdvectionRTheta operator can handle the advection with an advection field along $`(x,y)`$ and with an advection field along $`(r,\theta)`$. +But as the BslAdvectionRTheta operator advects in the physical domain, it is recommended to work with the advection field along $`(x,y)`$. ### Advection field along the physical domain axis Thanks to the spline representation, the derivatives $`\partial_r \phi`$ and $`\partial_\theta \phi`$ are computed. -The computation of the electrical field can be ill-defined around the O-point so we treat this area separately. +The computation of the electrical field can be ill-defined around the O-point, so we treat this area separately. * If $`r > \varepsilon`$, we use ```math @@ -83,7 +87,7 @@ A = E\wedge e_z \end{bmatrix}. ``` -* If $`r \leq \varepsilon`$, we linearise. The method is detailed in Edoardo Zoni's article [1]. We use only the derivatives along $`r`$ at two linearly independent directions of $`\theta`$ : $`\theta_1`$ and $`\theta_2`$ +* If $`r \leq \varepsilon`$, we linearize. The method is detailed in [Zoni et al. (2019)](#zoni). We use only the derivatives along $`r`$ at two linearly independent directions of $`\theta`$ : $`\theta_1`$ and $`\theta_2`$ ```math \partial_r \phi (0, \theta_1) = \left[\partial_r x \partial_x \phi + \partial_r y \partial_y \phi \right] (0, \theta_1), \\ \partial_r \phi (0, \theta_2) = \left[\partial_r x \partial_x \phi + \partial_r y \partial_y \phi \right] (0, \theta_2). @@ -107,7 +111,7 @@ From these equations, we deduce the (unique) values of $`\partial_x\phi`$ and $` \end{bmatrix}. ``` -Then we compute $`E`$ at $`(x,y) = (0,0)`$ and $`(x,y) = \mathcal{F}(\varepsilon,\theta)`$ $`\forall \theta`$ (for $`\varepsilon\neq 0`$, we use the Jacobian matrix as previously) and we linearise +Then we compute $`E`$ at $`(x,y) = (0,0)`$ and $`(x,y) = \mathcal{F}(\varepsilon,\theta)`$ $`\forall \theta`$ (for $`\varepsilon\neq 0`$, we use the Jacobian matrix as previously) and we linearize ```math E_x(r, \theta) = \left( 1 - \frac{r}{\varepsilon} \right) E_x(0, \theta) + \frac{r}{\varepsilon} E_x(\varepsilon, \theta), \\ @@ -137,50 +141,36 @@ A = E\wedge e_z Firstly, the derivatives $`\partial_r \phi`$ and $`\partial_\theta \phi`$ are also computed here. #### General coordinates system -* In **general coordinates system**, the gradiant of a function is given by +* In **general coordinates system**, the gradient of a function is given by ```math -\nabla f = \sum_i \sum_j \partial_{x_i} f g^{ij} \sqrt{g_{jj}} \hat{e}_j, +\nabla f = \sum_i \sum_j \partial_{x_i} f g^{ij} e_j ``` with -* $`J`$ the Jacobian matrix associated the the mapping function of the system $`\mathcal{F}:(x_1, x_2)\mapsto(y_1,y_2)`$, +* $`J`$ the Jacobian matrix associated to the mapping function of the system $`\mathcal{F}:(x_1, x_2)\mapsto(y_1,y_2)`$, * $`G = J^T J = [g_{ij}]_{ij}`$ the metric tensor, * $`G^{-1} = [g^{ij}]_{ij}`$ the inverse metric tensor -* and $`\hat{e}_j`$ the normalized covariant vectors. +* and $`e_j`$ the unnormalized local covariant vectors. In 2D, it can be rewritten as the following matrix system ```math \nabla f = -D_{G} G^{-T} +G^{-1} \begin{bmatrix} \partial_{x_1} f \\ \partial_{x_2} f \\ \end{bmatrix} -= -\begin{bmatrix} - \sqrt{g_{11}} & 0 \\ - 0 & \sqrt{g_{22}} \\ -\end{bmatrix} -\begin{bmatrix} - g^{11} & g^{21} \\ - g^{12} & g^{22} \\ -\end{bmatrix} -\begin{bmatrix} - \partial_{x_1} f \\ - \partial_{x_2} f \\ -\end{bmatrix}. ``` -**Remark:** We can prove that $`\det(D_{G} G^{-T}) = 0`$ if $`\det(J) = 0`$ (if the coefficients of the matrix $`D_G`$ are not null). -So for an invertible matrix, we also have the relation +So, for an invertible matrix, we also have the relation ```math \begin{bmatrix} \partial_{x_1} f \\ \partial_{x_2} f \\ \end{bmatrix} = -G^{T}D_{G}^{-1} \nabla f. +G \nabla f. ``` From the relation @@ -199,21 +189,24 @@ J^{-T} we deduce the following relation for invertible case ```math -\nabla_{y_1, y_2} f -= -J D_{G}^{-1} -\nabla_{x_1, x_2} f, +\nabla_{y_1, y_2} f = J \nabla_{x_1, x_2} f, ``` -with $`\nabla_{y_1, y_2} f = [\partial_{y_1} f, \partial_{y_2} f]^T`$ and $`\nabla_{x_1, x_2} f = \sum_i \sum_j \partial_{x_i} f g^{ij} \sqrt{g_{jj}} \hat{e}_j`$. +with $`\nabla_{y_1, y_2} f = [\partial_{y_1} f, \partial_{y_2} f]^T`$ and +$`\nabla_{x_1, x_2} f = \sum_i \sum_j \partial_{x_i} f g^{ij} e_j`$. -#### Aplication to the advection field +#### Application to the advection field * In our case, we use this formula to compute the electric field along the logical axis: ```math E -= -\nabla \phi -= - D_{G} (G^{-1})^{T} += +\begin{bmatrix} + E_r \\ + E_{\theta} \\ +\end{bmatrix} += - \nabla_{r,\theta} \phi += - G^{-1} \begin{bmatrix} \partial_{r} \phi \\ \partial_{\theta} \phi \\ @@ -231,11 +224,11 @@ A \end{bmatrix}. ``` -Warning, the matrix $`(G^{-1})^{T}`$ is ill-defined for $r = 0$. +Warning, the matrix $`G^{-1}`$ is ill-defined for $r = 0$. *Example: circular mapping:* ```math -(G^{-1})^{T} +G^{-1} = \begin{bmatrix} 1 & 0 \\ @@ -243,16 +236,16 @@ Warning, the matrix $`(G^{-1})^{T}`$ is ill-defined for $r = 0$. \end{bmatrix}. ``` -In the code, the O-point is differently treated. The domain is splitted between a domain without the O-point ($`(0,\theta), \forall \theta`$) and the domain containing only the O-point. For the first domain, we compute the advection field along the logical axis as explain previously. On the second domain, we compute the unique value of the advection field along the physical axis using the linearisation done in the *Advection field along the physical domain axis* section. +In the code, the O-point is differently treated. The domain is split between a domain without the O-point ($`(0,\theta), \forall \theta`$) and the domain containing only the O-point. For the first domain, we compute the advection field along the logical axis as explain previously. On the second domain, we compute the unique value of the advection field along the physical axis using the linearization done in the [Advection field along the physical domain axis](#src_geometryRTheta_advection_field__Guiding_center_case) section. # References -[1] Edoardo Zoni, Yaman Güçlü, "Solving hyperbolic-elliptic problems on singular mapped disk-like domains with the + [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, Journal of Computational Physics, 2019. # Contents -* advection\_field\_rp.hpp : containing AdvectionFieldFinder with the advection field computation for the guiding center simulation. +* advection\_field\_rtheta.hpp : containing AdvectionFieldFinder with the advection field computation for the guiding center simulation. diff --git a/src/geometryRTheta/advection_field/advection_field_rp.hpp b/src/geometryRTheta/advection_field/advection_field_rtheta.hpp similarity index 95% rename from src/geometryRTheta/advection_field/advection_field_rp.hpp rename to src/geometryRTheta/advection_field/advection_field_rtheta.hpp index aa8c5f7c3..5b6b79600 100644 --- a/src/geometryRTheta/advection_field/advection_field_rp.hpp +++ b/src/geometryRTheta/advection_field/advection_field_rtheta.hpp @@ -68,10 +68,10 @@ * * 2- In the second case, the advection field along the logical index range axis * is computed with - * - @f$ \nabla \phi = \sum_{i,j} \partial_{x_i} f g^{ij} \sqrt{g_{jj}} \hat{e}_j@f$, + * - @f$ \nabla \phi = \sum_{i,j} \partial_{x_i} f g^{ij} e_j@f$, * - with @f$g^{ij}@f$, the coefficients of the inverse metric tensor, * - @f$g_{jj}@f$, the coefficients of the metric tensor, - * - @f$\hat{e}_j@f$, the normalized covariants vectors. + * - @f$e_j@f$, the unnormalized local covariants vectors. * * Then, we compute @f$ E = -\nabla \phi @f$ and @f$A = E \wedge e_z@f$. * @@ -347,14 +347,14 @@ class AdvectionFieldFinder * * @param[in] electrostatic_potential * The values of the solution @f$\phi@f$ of the Poisson-like equation (2). - * @param[out] advection_field_rp + * @param[out] advection_field_rtheta * The advection field on the logical axis. * @param[out] advection_field_xy_center * The advection field on the physical axis at the O-point. */ void operator()( host_t electrostatic_potential, - host_t> advection_field_rp, + host_t> advection_field_rtheta, CoordXY& advection_field_xy_center) const { IdxRangeRTheta const grid = get_idx_range(electrostatic_potential); @@ -366,7 +366,7 @@ class AdvectionFieldFinder builder(get_field(electrostatic_potential_coef), get_const_field(electrostatic_potential)); (*this)(get_field(electrostatic_potential_coef), - advection_field_rp, + advection_field_rtheta, advection_field_xy_center); } @@ -378,20 +378,20 @@ class AdvectionFieldFinder * * @param[in] electrostatic_potential_coef * The spline representation of the solution @f$\phi@f$ of the Poisson-like equation (2). - * @param[out] advection_field_rp + * @param[out] advection_field_rtheta * The advection field on the logical axis. * @param[out] advection_field_xy_center * The advection field on the physical axis at the O-point. */ void operator()( host_t electrostatic_potential_coef, - host_t> advection_field_rp, + host_t> advection_field_rtheta, CoordXY& advection_field_xy_center) const { compute_advection_field_RTheta( m_spline_evaluator, electrostatic_potential_coef, - advection_field_rp, + advection_field_rtheta, advection_field_xy_center); } @@ -402,20 +402,20 @@ class AdvectionFieldFinder * * @param[in] electrostatic_potential_coef * The polar spline representation of the solution @f$\phi@f$ of the Poisson-like equation (2). - * @param[out] advection_field_rp + * @param[out] advection_field_rtheta * The advection field on the logical axis. * @param[out] advection_field_xy_center * The advection field on the physical axis at the O-point. */ void operator()( host_t& electrostatic_potential_coef, - host_t> advection_field_rp, + host_t> advection_field_rtheta, CoordXY& advection_field_xy_center) const { compute_advection_field_RTheta( m_polar_spline_evaluator, electrostatic_potential_coef, - advection_field_rp, + advection_field_rtheta, advection_field_xy_center); } @@ -429,7 +429,7 @@ class AdvectionFieldFinder * The spline evaluator used to evaluated electrostatic_potential_coef. * @param[in] electrostatic_potential_coef * The spline representation of the solution @f$\phi@f$ of the Poisson-like equation (2). - * @param[out] advection_field_rp + * @param[out] advection_field_rtheta * The advection field on the logical axis on an index range without O-point. * @param[out] advection_field_xy_center * The advection field on the physical axis at the O-point. @@ -438,7 +438,7 @@ class AdvectionFieldFinder void compute_advection_field_RTheta( Evaluator evaluator, SplineType& electrostatic_potential_coef, - host_t> advection_field_rp, + host_t> advection_field_rtheta, CoordXY& advection_field_xy_center) const { static_assert( @@ -451,7 +451,7 @@ class AdvectionFieldFinder PolarBSplinesRTheta, ddc::NullExtrapolationRule>> && std::is_same_v>)); - IdxRangeRTheta const grid_without_Opoint = get_idx_range(advection_field_rp); + IdxRangeRTheta const grid_without_Opoint = get_idx_range(advection_field_rtheta); host_t> coords(grid_without_Opoint); ddc::for_each(grid_without_Opoint, [&](IdxRTheta const irp) { @@ -477,24 +477,18 @@ class AdvectionFieldFinder ddc::for_each(grid_without_Opoint, [&](IdxRTheta const irp) { CoordRTheta const coord_rp(ddc::coordinate(irp)); - Matrix_2x2 J; // Jacobian matrix - m_mapping.jacobian_matrix(coord_rp, J); Matrix_2x2 inv_G; // Inverse metric tensor metric_tensor.inverse(inv_G, coord_rp); - Matrix_2x2 G; // Metric tensor - metric_tensor(G, coord_rp); // E = -grad phi double const electric_field_r - = (-deriv_r_phi(irp) * inv_G[0][0] - deriv_p_phi(irp) * inv_G[1][0]) - * std::sqrt(G[0][0]); + = -deriv_r_phi(irp) * inv_G[0][0] - deriv_p_phi(irp) * inv_G[0][1]; double const electric_field_p - = (-deriv_r_phi(irp) * inv_G[0][1] - deriv_p_phi(irp) * inv_G[1][1]) - * std::sqrt(G[1][1]); + = -deriv_r_phi(irp) * inv_G[1][0] - deriv_p_phi(irp) * inv_G[1][1]; // A = E \wedge e_z - ddcHelper::get(advection_field_rp)(irp) = -electric_field_p; - ddcHelper::get(advection_field_rp)(irp) = electric_field_r; + ddcHelper::get(advection_field_rtheta)(irp) = -electric_field_p; + ddcHelper::get(advection_field_rtheta)(irp) = electric_field_r; }); // SPECIAL TREATMENT FOR THE O-POINT ===================================================== diff --git a/src/geometryRTheta/time_solver/CMakeLists.txt b/src/geometryRTheta/time_solver/CMakeLists.txt index 1948f5b44..e21d7f298 100644 --- a/src/geometryRTheta/time_solver/CMakeLists.txt +++ b/src/geometryRTheta/time_solver/CMakeLists.txt @@ -13,7 +13,7 @@ target_link_libraries(time_integration_rp gslx::interpolation_2D_rp gslx::geometry_RTheta gslx::poisson_RTheta - gslx::advection_rp + gslx::advection_polar gslx::utils gslx::timestepper gslx::advection_field_RTheta diff --git a/src/geometryRTheta/time_solver/README.md b/src/geometryRTheta/time_solver/README.md index 380a5e942..e30528d73 100644 --- a/src/geometryRTheta/time_solver/README.md +++ b/src/geometryRTheta/time_solver/README.md @@ -33,7 +33,7 @@ Advect on a half time step: 2. From $\phi^n$, we compute $E^n$ by deriving (IQNSolver); - 3. From $\rho^n \text{ and } A^n$, we compute $\rho^{n+1/2}$ by advecting (IAdvectionRP) on $\frac{dt}{2}$ with one of the selected time integration methods (ITimeStepper); + 3. From $\rho^n \text{ and } A^n$, we compute $\rho^{n+1/2}$ by advecting (IAdvectionRTheta) on $\frac{dt}{2}$ with one of the selected time integration methods (ITimeStepper); Advect on a full time step: @@ -41,7 +41,7 @@ Advect on a full time step: 5. From $\phi^{n+1/2}$, we compute $E^{n+1/2}$ by deriving (IQNSolver); - 6. From $\rho^n \text{ and } A^{n+1/2}$, we compute $\rho^{n+1}$ by advecting (IAdvectionRP) on $dt$ with one of the selected time integration methods (ITimeStepper). + 6. From $\rho^n \text{ and } A^{n+1/2}$, we compute $\rho^{n+1}$ by advecting (IAdvectionRTheta) on $dt$ with one of the selected time integration methods (ITimeStepper). ## Explicit predictor-corrector @@ -54,7 +54,7 @@ Prediction: 2. From $\phi^n$, we compute $E^n$ by deriving (IQNSolver); - 3. From $\rho^n \text{ and } A^n$, we compute $\rho^P$ by advecting (IAdvectionRP) on $dt$; + 3. From $\rho^n \text{ and } A^n$, we compute $\rho^P$ by advecting (IAdvectionRTheta) on $dt$; - We write $X^P$ the characteristic feet such that $`\partial_t X^P = A^n(X^n)`$. Correction: @@ -63,7 +63,7 @@ Correction: 5. From $\phi^{P}$, we compute $E^{P}$ by deriving (IQNSolver); - 6. From $\rho^n \text{ and } \frac{A^{P}(X^n) + A^n(X^P)}{2}$, we compute $\rho^{C} = \rho^{n+1}$ by advecting (IAdvectionRP) on $dt$. + 6. From $\rho^n \text{ and } \frac{A^{P}(X^n) + A^n(X^P)}{2}$, we compute $\rho^{C} = \rho^{n+1}$ by advecting (IAdvectionRTheta) on $dt$. @@ -78,7 +78,7 @@ Prediction: 2. From $\phi^n$, we compute $E^n$ by deriving (IQNSolver); - 3. From $\rho^n \text{ and } A^n$, we compute $\rho^P$ by advecting (IAdvectionRP) on $\frac{dt}{2}$; + 3. From $\rho^n \text{ and } A^n$, we compute $\rho^P$ by advecting (IAdvectionRTheta) on $\frac{dt}{2}$; - We write $X^P$ the characteristic feet such that it is the result of the implicit method: ```math \partial_t X^k = \frac{A^n(X^n) + A^n(X^{k-1})}{2}, \qquad X^k = X^n - \frac{dt}{2} \partial_t X^k. @@ -91,7 +91,7 @@ Correction: 5. From $\phi^{P}$, we compute $E^{P}$ by deriving (IQNSolver); - 6. From $\rho^n \text{ and } A^{P}$, we compute $\rho^{C} = \rho^{n+1}$ by advecting (IAdvectionRP) on $dt$. + 6. From $\rho^n \text{ and } A^{P}$, we compute $\rho^{C} = \rho^{n+1}$ by advecting (IAdvectionRTheta) on $dt$. - We write $X^C$ the characteristic feet such that it is the result of the implicit method: ```math \partial_t X^k = \frac{A^P(X^n) + A^P(X^{k-1})}{2}, \qquad X^k = X^n - dt \partial_t X^k. diff --git a/src/geometryRTheta/time_solver/bsl_predcorr.hpp b/src/geometryRTheta/time_solver/bsl_predcorr.hpp index 99ed46b99..a043c3697 100644 --- a/src/geometryRTheta/time_solver/bsl_predcorr.hpp +++ b/src/geometryRTheta/time_solver/bsl_predcorr.hpp @@ -9,8 +9,8 @@ #include #include -#include "advection_field_rp.hpp" -#include "bsl_advection_rp.hpp" +#include "advection_field_rtheta.hpp" +#include "bsl_advection_rtheta.hpp" #include "ddc_alias_inline_functions.hpp" #include "ddc_aliases.hpp" #include "geometry.hpp" diff --git a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp index a19031edf..fa3a59de6 100644 --- a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp +++ b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_explicit.hpp @@ -9,8 +9,8 @@ #include #include -#include "advection_field_rp.hpp" -#include "bsl_advection_rp.hpp" +#include "advection_field_rtheta.hpp" +#include "bsl_advection_rtheta.hpp" #include "ddc_alias_inline_functions.hpp" #include "ddc_aliases.hpp" #include "euler.hpp" diff --git a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp index 50f26eab7..434b67661 100644 --- a/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp +++ b/src/geometryRTheta/time_solver/bsl_predcorr_second_order_implicit.hpp @@ -9,8 +9,8 @@ #include #include -#include "advection_field_rp.hpp" -#include "bsl_advection_rp.hpp" +#include "advection_field_rtheta.hpp" +#include "bsl_advection_rtheta.hpp" #include "ddc_alias_inline_functions.hpp" #include "ddc_aliases.hpp" #include "geometry.hpp" diff --git a/tests/geometryRTheta/CMakeLists.txt b/tests/geometryRTheta/CMakeLists.txt index 69dcadcab..756a04992 100644 --- a/tests/geometryRTheta/CMakeLists.txt +++ b/tests/geometryRTheta/CMakeLists.txt @@ -8,5 +8,5 @@ if("${POISSON_2D_BUILD_TESTING}") endif() add_subdirectory(2d_spline_interpolator) add_subdirectory(quadrature) -add_subdirectory(advection_2d_rp) -add_subdirectory(advection_field_rp) +add_subdirectory(advection_2d_rtheta) +add_subdirectory(advection_field_rtheta) diff --git a/tests/geometryRTheta/README.md b/tests/geometryRTheta/README.md index 7bfee14cf..fedb4466f 100644 --- a/tests/geometryRTheta/README.md +++ b/tests/geometryRTheta/README.md @@ -2,7 +2,7 @@ The `tests/geometryRTheta` folder contains all the code describing methods which are specific to a 2D curvilinear geometry containing a singular point. It is broken up into the following sub-folders: -- [advection\_2d\_rp](./advection_2d_rp/README.md) - Tests for the advection operator and time integration methods used on 2D polar domain. +- [advection\_2d\_rtheta](./advection_2d_rtheta/README.md) - Tests for the advection operator and time integration methods used on 2D polar domain. - [2d\_spline\_interpolator](./2d_spline_interpolator/README.md) - Tests for the interpolator on 2D polar domain. diff --git a/tests/geometryRTheta/advection_2d_rp/CMakeLists.txt b/tests/geometryRTheta/advection_2d_rtheta/CMakeLists.txt similarity index 98% rename from tests/geometryRTheta/advection_2d_rp/CMakeLists.txt rename to tests/geometryRTheta/advection_2d_rtheta/CMakeLists.txt index 9b36ff3a1..791660e35 100644 --- a/tests/geometryRTheta/advection_2d_rp/CMakeLists.txt +++ b/tests/geometryRTheta/advection_2d_rtheta/CMakeLists.txt @@ -35,7 +35,7 @@ foreach(MAPPING_TYPE "CIRCULAR_MAPPING_PHYSICAL" "CZARNY_MAPPING_PHYSICAL" "CZAR Eigen3::Eigen - gslx::advection_rp + gslx::advection_polar gslx::geometry_RTheta gslx::interpolation_2D_rp gslx::io @@ -74,7 +74,7 @@ target_link_libraries(advection_ALL Eigen3::Eigen - gslx::advection_rp + gslx::advection_polar gslx::geometry_RTheta gslx::interpolation_2D_rp gslx::io diff --git a/tests/geometryRTheta/advection_2d_rp/README.md b/tests/geometryRTheta/advection_2d_rtheta/README.md similarity index 92% rename from tests/geometryRTheta/advection_2d_rp/README.md rename to tests/geometryRTheta/advection_2d_rtheta/README.md index c0a46b87e..915b2dd4f 100644 --- a/tests/geometryRTheta/advection_2d_rp/README.md +++ b/tests/geometryRTheta/advection_2d_rtheta/README.md @@ -38,17 +38,17 @@ The tests are made for different parameters which are: ## Python tests - animated\_curves.py: create `.mp4` video(s) of the advected function for the selected configuations among the 48 test ones. - - Command to launch the test in this folder: `python3 animated_curves.py ../../../build/tests/geometryRTheta/advection_2d_rp/advection_ALL` - or `python3 animated_curves.py ../../../build/tests/geometryRTheta/advection_2d_rp/` + - Command to launch the test in this folder: `python3 animated_curves.py ../../../build/tests/geometryRTheta/advection_2d_rtheta/advection_ALL` + or `python3 animated_curves.py ../../../build/tests/geometryRTheta/advection_2d_rtheta/` - display\_all\_errors\_for\_gtest.py: Google test which tests the convergence order for the 48 configurations. - - Command to launch the test in this folder: `python3 display_all_errors_for_gtest.py ../../../build/tests/geometryRTheta/advection_2d_rp/advection_ALL` + - Command to launch the test in this folder: `python3 display_all_errors_for_gtest.py ../../../build/tests/geometryRTheta/advection_2d_rtheta/advection_ALL` - display\_curves.py: display the curve of the function for 9 time steps between the initial and the final state. - - Command to launch the test in this folder: `python3 display_curves.py ../../../build/tests/geometryRTheta/advection_2d_rp/` + - Command to launch the test in this folder: `python3 display_curves.py ../../../build/tests/geometryRTheta/advection_2d_rtheta/` - display\_feet\_errors.py: compute the convergence order of the characteristic feet and display the computed and the exact feet. - - Command to launch the test in this folder: `python3 display_feet_errors.py ../../../build/tests/geometryRTheta/advection_2d_rp/` + - Command to launch the test in this folder: `python3 display_feet_errors.py ../../../build/tests/geometryRTheta/advection_2d_rtheta/` - advection\_functions.py: define all the useful functions used in the other python files. diff --git a/tests/geometryRTheta/advection_2d_rp/advection_all_tests.cpp b/tests/geometryRTheta/advection_2d_rtheta/advection_all_tests.cpp similarity index 99% rename from tests/geometryRTheta/advection_2d_rp/advection_all_tests.cpp rename to tests/geometryRTheta/advection_2d_rtheta/advection_all_tests.cpp index c0e65418b..18ee69e53 100644 --- a/tests/geometryRTheta/advection_2d_rp/advection_all_tests.cpp +++ b/tests/geometryRTheta/advection_2d_rtheta/advection_all_tests.cpp @@ -12,7 +12,7 @@ #include #include "advection_simulation_utils.hpp" -#include "bsl_advection_rp.hpp" +#include "bsl_advection_rtheta.hpp" #include "cartesian_to_circular.hpp" #include "cartesian_to_czarny.hpp" #include "circular_to_cartesian.hpp" diff --git a/tests/geometryRTheta/advection_2d_rp/advection_functions.py b/tests/geometryRTheta/advection_2d_rtheta/advection_functions.py similarity index 99% rename from tests/geometryRTheta/advection_2d_rp/advection_functions.py rename to tests/geometryRTheta/advection_2d_rtheta/advection_functions.py index 61fa0c3b6..322c43d85 100644 --- a/tests/geometryRTheta/advection_2d_rp/advection_functions.py +++ b/tests/geometryRTheta/advection_2d_rtheta/advection_functions.py @@ -1,5 +1,5 @@ """ -Define all the functions needed in the python files of the advection_2d_rp folder. +Define all the functions needed in the python files of the advection_2d_rtheta folder. """ import argparse diff --git a/tests/geometryRTheta/advection_2d_rp/advection_selected_test.cpp b/tests/geometryRTheta/advection_2d_rtheta/advection_selected_test.cpp similarity index 99% rename from tests/geometryRTheta/advection_2d_rp/advection_selected_test.cpp rename to tests/geometryRTheta/advection_2d_rtheta/advection_selected_test.cpp index 24e62131e..ebb18bacf 100644 --- a/tests/geometryRTheta/advection_2d_rp/advection_selected_test.cpp +++ b/tests/geometryRTheta/advection_2d_rtheta/advection_selected_test.cpp @@ -9,7 +9,7 @@ #include #include "advection_simulation_utils.hpp" -#include "bsl_advection_rp.hpp" +#include "bsl_advection_rtheta.hpp" #include "cartesian_to_circular.hpp" #include "cartesian_to_czarny.hpp" #include "circular_to_cartesian.hpp" diff --git a/tests/geometryRTheta/advection_2d_rp/advection_simulation_utils.hpp b/tests/geometryRTheta/advection_2d_rtheta/advection_simulation_utils.hpp similarity index 99% rename from tests/geometryRTheta/advection_2d_rp/advection_simulation_utils.hpp rename to tests/geometryRTheta/advection_2d_rtheta/advection_simulation_utils.hpp index a2860035c..f9ce53c7d 100644 --- a/tests/geometryRTheta/advection_2d_rp/advection_simulation_utils.hpp +++ b/tests/geometryRTheta/advection_2d_rtheta/advection_simulation_utils.hpp @@ -9,7 +9,7 @@ #include -#include "bsl_advection_rp.hpp" +#include "bsl_advection_rtheta.hpp" #include "directional_tag.hpp" #include "geometry.hpp" #include "l_norm_tools.hpp" diff --git a/tests/geometryRTheta/advection_2d_rp/animated_curves.py b/tests/geometryRTheta/advection_2d_rtheta/animated_curves.py similarity index 100% rename from tests/geometryRTheta/advection_2d_rp/animated_curves.py rename to tests/geometryRTheta/advection_2d_rtheta/animated_curves.py diff --git a/tests/geometryRTheta/advection_2d_rp/display_all_errors_for_gtest.py b/tests/geometryRTheta/advection_2d_rtheta/display_all_errors_for_gtest.py similarity index 99% rename from tests/geometryRTheta/advection_2d_rp/display_all_errors_for_gtest.py rename to tests/geometryRTheta/advection_2d_rtheta/display_all_errors_for_gtest.py index 2a2644808..974e1a14c 100644 --- a/tests/geometryRTheta/advection_2d_rp/display_all_errors_for_gtest.py +++ b/tests/geometryRTheta/advection_2d_rtheta/display_all_errors_for_gtest.py @@ -1,5 +1,5 @@ # Command to launch the test : -# python3 display_all_error_for_gtest.py ../../../build/tests/geometryRTheta/advection_2d_rp/advection_ALL +# python3 display_all_error_for_gtest.py ../../../build/tests/geometryRTheta/advection_2d_rtheta/advection_ALL """ Compute the space convergence order between two grid sizes diff --git a/tests/geometryRTheta/advection_2d_rp/display_curves.py b/tests/geometryRTheta/advection_2d_rtheta/display_curves.py similarity index 99% rename from tests/geometryRTheta/advection_2d_rp/display_curves.py rename to tests/geometryRTheta/advection_2d_rtheta/display_curves.py index 1ecc1c557..7837c6a56 100644 --- a/tests/geometryRTheta/advection_2d_rp/display_curves.py +++ b/tests/geometryRTheta/advection_2d_rtheta/display_curves.py @@ -1,5 +1,5 @@ # Command to launch the test : -# python3 display_curves.py ../../../build/tests/geometryRTheta/advection_2d_rp/ +# python3 display_curves.py ../../../build/tests/geometryRTheta/advection_2d_rtheta/ """ Plot the advected function at regular different time steps. diff --git a/tests/geometryRTheta/advection_2d_rp/display_feet_errors.py b/tests/geometryRTheta/advection_2d_rtheta/display_feet_errors.py similarity index 99% rename from tests/geometryRTheta/advection_2d_rp/display_feet_errors.py rename to tests/geometryRTheta/advection_2d_rtheta/display_feet_errors.py index deffe932a..93d28cce9 100644 --- a/tests/geometryRTheta/advection_2d_rp/display_feet_errors.py +++ b/tests/geometryRTheta/advection_2d_rtheta/display_feet_errors.py @@ -1,5 +1,5 @@ # Command to launch the test : -# python3 display_feet_errors.py ../../../build/tests/geometryRTheta/advection_2d_rp/advection_2d_rp_tests +# python3 display_feet_errors.py ../../../build/tests/geometryRTheta/advection_2d_rtheta/advection_2d_rtheta_tests """ Display the characteristic feet and their errors computed at the last time step of the advection diff --git a/tests/geometryRTheta/advection_2d_rp/params.yaml b/tests/geometryRTheta/advection_2d_rtheta/params.yaml similarity index 100% rename from tests/geometryRTheta/advection_2d_rp/params.yaml rename to tests/geometryRTheta/advection_2d_rtheta/params.yaml diff --git a/tests/geometryRTheta/advection_2d_rp/params.yaml.hpp b/tests/geometryRTheta/advection_2d_rtheta/params.yaml.hpp similarity index 100% rename from tests/geometryRTheta/advection_2d_rp/params.yaml.hpp rename to tests/geometryRTheta/advection_2d_rtheta/params.yaml.hpp diff --git a/tests/geometryRTheta/advection_2d_rp/test_cases.hpp b/tests/geometryRTheta/advection_2d_rtheta/test_cases.hpp similarity index 100% rename from tests/geometryRTheta/advection_2d_rp/test_cases.hpp rename to tests/geometryRTheta/advection_2d_rtheta/test_cases.hpp diff --git a/tests/geometryRTheta/advection_field_rp/CMakeLists.txt b/tests/geometryRTheta/advection_field_rtheta/CMakeLists.txt similarity index 98% rename from tests/geometryRTheta/advection_field_rp/CMakeLists.txt rename to tests/geometryRTheta/advection_field_rtheta/CMakeLists.txt index d77845be7..43c392eb5 100644 --- a/tests/geometryRTheta/advection_field_rp/CMakeLists.txt +++ b/tests/geometryRTheta/advection_field_rtheta/CMakeLists.txt @@ -32,7 +32,7 @@ function (advection_field_test_executable SIMULATION) gslx::poisson_RTheta gslx::geometry_RTheta gslx::interpolation_2D_rp - gslx::advection_rp + gslx::advection_polar gslx::time_integration_rp gslx::quadrature gslx::utils diff --git a/tests/geometryRTheta/advection_field_rp/advection_field_gtest.cpp b/tests/geometryRTheta/advection_field_rtheta/advection_field_gtest.cpp similarity index 93% rename from tests/geometryRTheta/advection_field_rp/advection_field_gtest.cpp rename to tests/geometryRTheta/advection_field_rtheta/advection_field_gtest.cpp index 29137fad6..758b2e574 100644 --- a/tests/geometryRTheta/advection_field_rp/advection_field_gtest.cpp +++ b/tests/geometryRTheta/advection_field_rtheta/advection_field_gtest.cpp @@ -15,7 +15,7 @@ #include #include -#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" @@ -181,7 +181,7 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) host_t allfdistribu_xy_alloc(grid); host_t> advection_field_exact_alloc(grid); - host_t> advection_field_rp_alloc(grid_without_Opoint); + host_t> advection_field_rtheta_alloc(grid_without_Opoint); host_t> advection_field_xy_alloc(grid); host_t> advection_field_xy_from_rp_alloc(grid); CoordXY advection_field_xy_center; @@ -193,7 +193,7 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) host_t electrostatic_potential(electrostatic_potential_alloc); host_t> advection_field_exact(advection_field_exact_alloc); - host_t> advection_field_rp(advection_field_rp_alloc); + host_t> advection_field_rtheta(advection_field_rtheta_alloc); host_t> advection_field_xy(advection_field_xy_alloc); host_t> advection_field_xy_from_rp(advection_field_xy_from_rp_alloc); @@ -220,7 +220,7 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) // Constant advection fields ************************************* advection_field_computer( electrostatic_potential, - advection_field_rp, + advection_field_rtheta, advection_field_xy_center); advection_field_computer(electrostatic_potential, advection_field_xy); @@ -245,18 +245,17 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) ddc::for_each(grid_without_Opoint, [&](IdxRTheta const irp) { CoordRTheta const coord_rp(ddc::coordinate(irp)); - std::array, 2> J; // Jacobian matrix - to_physical_mapping.jacobian_matrix(coord_rp, J); - std::array, 2> G; // Metric tensor - metric_tensor(G, coord_rp); + std::array, 2> inv_J; // inverse Jacobian matrix + to_physical_mapping.inv_jacobian_matrix(coord_rp, inv_J); + double const jacobian = to_physical_mapping.jacobian(coord_rp); // computation made in BslAdvectionRTheta operator: ddcHelper::get(advection_field_xy_from_rp)(irp) - = ddcHelper::get(advection_field_rp)(irp) * J[1][1] / std::sqrt(G[1][1]) - + ddcHelper::get(advection_field_rp)(irp) * -J[1][0] / std::sqrt(G[0][0]); + = ddcHelper::get(advection_field_rtheta)(irp) * inv_J[0][0] * jacobian + + ddcHelper::get(advection_field_rtheta)(irp) * inv_J[1][0] * jacobian; ddcHelper::get(advection_field_xy_from_rp)(irp) - = ddcHelper::get(advection_field_rp)(irp) * -J[0][1] / std::sqrt(G[1][1]) - + ddcHelper::get(advection_field_rp)(irp) * J[0][0] / std::sqrt(G[0][0]); + = ddcHelper::get(advection_field_rtheta)(irp) * inv_J[0][1] * jacobian + + ddcHelper::get(advection_field_rtheta)(irp) * inv_J[1][1] * jacobian; // compare ddcHelper::get(difference_between_fields_xy_and_rp)(irp) @@ -295,7 +294,7 @@ TEST(AdvectionFieldRThetaComputation, TestAdvectionFieldFinder) // SIMULATION | // ================================================================================================ for (int iter(0); iter < iter_nb; ++iter) { - advection_operator(allfdistribu_rp, advection_field_rp, advection_field_xy_center, dt); + advection_operator(allfdistribu_rp, advection_field_rtheta, advection_field_xy_center, dt); advection_operator(allfdistribu_xy, advection_field_xy, dt); // Check the advected functions --- diff --git a/tests/geometryRTheta/advection_field_rp/test_cases_adv_field.hpp b/tests/geometryRTheta/advection_field_rtheta/test_cases_adv_field.hpp similarity index 99% rename from tests/geometryRTheta/advection_field_rp/test_cases_adv_field.hpp rename to tests/geometryRTheta/advection_field_rtheta/test_cases_adv_field.hpp index d80c311e5..7c7271335 100644 --- a/tests/geometryRTheta/advection_field_rp/test_cases_adv_field.hpp +++ b/tests/geometryRTheta/advection_field_rtheta/test_cases_adv_field.hpp @@ -2,7 +2,7 @@ #pragma once #include -#include "../advection_2d_rp/test_cases.hpp" +#include "../advection_2d_rtheta/test_cases.hpp" #include "circular_to_cartesian.hpp" #include "ddc_aliases.hpp"