diff --git a/src/geometryRTheta/README.md b/src/geometryRTheta/README.md index d30a588b3..e746e1409 100644 --- a/src/geometryRTheta/README.md +++ b/src/geometryRTheta/README.md @@ -8,6 +8,6 @@ The `geometryRTheta` folder contains all the code describing methods which are s - [interpolation](./interpolation/README.md) - Code describing interpolation methods on 2D polar domain. - - Code describing the polar Poisson solver. +- [poisson](./poisson/README.md) - Code describing the polar Poisson solver. diff --git a/src/geometryRTheta/geometry/geometry.hpp b/src/geometryRTheta/geometry/geometry.hpp index d01439347..ac013e3d3 100644 --- a/src/geometryRTheta/geometry/geometry.hpp +++ b/src/geometryRTheta/geometry/geometry.hpp @@ -171,6 +171,23 @@ using Spline2D = ddc::Chunk; using Spline2DSpan = ddc::ChunkSpan; using Spline2DView = ddc::ChunkSpan; +/** + * @brief Tag the polar B-splines decomposition of a function. + * + * Store the polar B-splines coefficients of the function. + */ +using SplinePolar = PolarSpline; + +/** + * @brief Tag type of an element of a couple of a B-splines in the first + * dimension and in the second dimension. + */ +using IDimBSpline2D = ddc::DiscreteElement; +/** + * @brief Tag type of an element of polar B-splines. + */ +using IDimPolarBspl = ddc::DiscreteElement; + template diff --git a/src/geometryRTheta/poisson/CMakeLists.txt b/src/geometryRTheta/poisson/CMakeLists.txt index 4ed9abd06..540e174bd 100644 --- a/src/geometryRTheta/poisson/CMakeLists.txt +++ b/src/geometryRTheta/poisson/CMakeLists.txt @@ -1,24 +1,18 @@ # SPDX-License-Identifier: MIT -add_library("poisson_RTheta" STATIC - nullpoissonsolver.cpp -) - -target_compile_features("poisson_RTheta" - PUBLIC - cxx_std_17 -) -target_include_directories("poisson_RTheta" - PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}" +add_library(poisson_RTheta INTERFACE) +target_include_directories(poisson_RTheta + INTERFACE + "$" ) - -target_link_libraries("poisson_RTheta" - PUBLIC +target_link_libraries(poisson_RTheta INTERFACE DDC::DDC sll::splines gslx::geometry_RTheta gslx::speciesinfo + gslx::utils + gslx::interpolation_2D_rp ) add_library("gslx::poisson_RTheta" ALIAS "poisson_RTheta") diff --git a/src/geometryRTheta/poisson/README.md b/src/geometryRTheta/poisson/README.md new file mode 100644 index 000000000..1a2e7aa65 --- /dev/null +++ b/src/geometryRTheta/poisson/README.md @@ -0,0 +1,137 @@ +# Polar Poisson solver + + +## The Poisson equation + +(For more details, see Emily Bourne's thesis "Non-Uniform Numerical Schemes for the Modelling of Turbulence +in the 5D GYSELA Code". December 2022.) + +The equation we are solving here is + +$$L\phi = - \nabla \cdot (\alpha \nabla \phi) + \beta \phi = \rho$$ + +with the boundary condition $\phi = 0$, on $\partial \Omega$. + + +To solve this equation, the PolarSplineFEMPoissonSolver uses a finite element method on the B-splines. + +### B-splines FEM + +#### The B-splines +We introduce a basis of B-splines $\{B_l\}_l$, cross-product of two 1D B-splines basis: + +$$\{B_l\}_l = \{b_{i,r} b_{j,\theta}\}_{i n_{n,\theta} +j}.$$ + + +#### Treatment of the O-point + +To conserve $\mathcal{C}^k$ property at the O-point, the B-splines basis is modified. + +The $\frac{(k+1)(k+2)}{2}$ first elements in the 1D B-splines basis in the $r$ dimension are removed, and +replaced by 3 functions, linear combinations of these removed B-splines. +(See (2.34) in Emily Bourne's thesis.) + +The linear composition uses as coefficients, the Berstein polynomials defined from barycentric coordinates. + +The B-splines with the treatment of the O-point, are called PolarBSplines and are written $\{\hat{B}_l\}_l$. + + +#### Weak formulation +The Poisson equation is solved by solving its weak form: + +$$\int_{\Omega} \lbrack \beta(r) \phi(r,\theta) \hat{B}_l(r,\theta) + \alpha(r) \nabla \phi(r,\theta) \cdot \nabla \hat{B}_l(r,\theta) \rbrack |det(J_{\mathcal{F}}(r,\theta))| dr d\theta = \int_{\Omega} \rho(r,\theta) \hat{B}_l(r,\theta) |det(J_{\mathcal{F}}(r,\theta))| dr d\theta$$ + +with $\{\hat{B}_l\}_l$ the B-splines basis. + + +Written as matrix form, it gives + +$$(M + S) \hat{\phi} = M \hat{\rho},$$ + +with + * the mass matrix, $M_{i,j} = \int_{\Omega} \beta(r,\theta) \hat{B}_i(r,\theta)\hat{B}_j(r,\theta) |det(J_{\mathcal{F}}(r,\theta))| dr d\theta$, + * the stiffness matrix, $S_{i,j} = \int_{\Omega} \alpha(r) \left[\sum_{\xi_1\in[r,\theta]} \sum_{\xi_2\in[r,\theta]} \partial_{\xi_1}\hat{B}_i(r,\theta) \partial_{\xi_2}\hat{B}_j(r,\theta) g^{\xi_1, \xi_2}\right] |det(J_{\mathcal{F}}(r,\theta))| dr d\theta$ + * with $g^{\xi_1, \xi_2}$, the scalar product between the unit vector in $\xi_1$ and the unit vector in $\xi_2$. + * the solution vector, $\hat{\phi}_i = \sum_{l = 0}^{n_{b,\theta}(n_{b,r} -2) +3} \phi_l \hat{B}_l(r_i,\theta_i)$, + * and the rhs vector, $\hat{\rho}_j = \sum_{l = 0}^{n_{b,\theta}(n_{b,r} -2) +3} \rho_l \hat{B}_l(r_j,\theta_j)$. + +So we compute the solution B-splines coefficients $\{\phi_l\}_l$ by solving this matrix equation. + + + + + + +## Evaluation of electric field + +(See for more details, Edoardo Zoni's article, https://doi.org/10.1016/j.jcp.2019.108889 .) + +We write in this section $\nabla_{r,\theta}$ the gradient in the logical coordinates $(r, \theta)$ +and $\nabla_{x,y}$ the gradient in the physical coordinates $(x, y)$. + +Coupled with the Vlasov equation +$$\partial_t \rho - E_y \partial_x \rho + E_x \partial_y\rho = 0,$$ + +the VlasovPoissonSolver also computes the electric field $E$ in the physical domain from the solution of PolarSplineFEMPoissonSolver. + +The electric field is given by $E = -\nabla \phi$. The solution of the Poisson solver is defined +on the logical domain, so we can easily compute $\nabla_{r,\theta} \phi$. From that, we use the Jacobian +matrix of the mapping $\mathcal{F}$: + +$$ E_r e_r + E_\theta e_\theta = -\nabla_{r,\theta} \phi, $$ + +$$ E_x e_x + E_y e_y = J_{\mathcal{F}}(r,\theta)) (E_r e_r + E_\theta e_\theta), $$ + +with $e_i = \partial_{x_i}x$ the unnormalized local contravariant base. + +However the inverse Jacobian matrix $J_{\mathcal{F}}$ can be ill-defined at the O-point. In +Edoardo Zoni's article, they suggest to linearize around the O-point: + +for $r < \varepsilon$, +$$E(r, \theta) = \left( 1 - \frac{r}{\varepsilon} \right) E(0, \theta) + \frac{r}{\varepsilon} E(\varepsilon, \theta)$$ + +with $E(0, \theta)$ computed thanks to +* $\partial_r \phi (0, \theta_1) = \left[\partial_r x \partial_x \phi + \partial_r y \partial_y \phi \right] (0, \theta_1)$, and +* $\partial_r \phi (0, \theta_2) = \left[\partial_r x \partial_x \phi + \partial_r y \partial_y \phi \right] (0, \theta_2)$, +* where $\theta_1$ and $\theta_2$ correspond to linearly independent directions. + + +(In the code, we chose $\theta_1 = \frac{\pi}{4}$ and $\theta_2 = - \frac{\pi}{4}$, and $\varepsilon = 10^{-12}$.) + + +## Unit tests + +The test are implemented in the `tests/geometryRTheta/polar_poisson/` folder +([polar_poisson](./../../../tests/geometryRTheta/polar_poisson/README.md)). + +The PolarSplineFEMPoissonSolver is tested on a circular mapping (CircularToCartesian) and on a Czarny mapping (CzarnyToCartesian) with +(the test cases are given in Emily Bourne's thesis [1]) + + * the Poisson coefficients: + * $\alpha(r) = \exp\left( - \tanh\left( \frac{r - r_p}{\delta_r} \right) \right)$, with $r_p = 0.7$ and $\delta_r = 0.05$, + * $\beta(r) = \frac{1}{\alpha(r)}$. + * for the following test functions: + * polar solution: $\phi(x, y) = C r(x,y)^6 (r(x,y) -1)^6 \cos(m\theta)$, + * with $C = 2^{12}1e-4$ and $m = 11$; + * cartesian solution: $\phi(x,y) = C (1+r(x,y))^6 (1 - r(x,y))^6 \cos(2\pi x) \sin(2\pi y)$, + * with $C = 2^{12}1e-4$. + + The VlasovPoissonSolver is tested on a circular mapping (CircularToCartesian) and on a Czarny mapping (CzarnyToCartesian) + with the same Poisson coeffiecients and for the cartesian solution. + + +## References + +[1] Emily Bourne, "Non-Uniform Numerical Schemes for the Modelling of Turbulence in the 5D GYSELA Code". December 2022. + +[2] 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 + + * ipoissonsolver.hpp : Define a base class for the Poisson solvers: IPoissonSolver. + * polarpoissonsolver.hpp : Define a Poisson solver using FEM on B-splines: PolarSplineFEMPoissonSolver. + * poisson_rhs_function.hpp : Define a rhs object (PoissonRHSFunction) for the Poisson equation (mainly used for vlasovpoissonsolver.hpp): PoissonRHSFunction. + * vlasovpoissonsolver.hpp : Define a class which solves the Poisson equation and computes the electric field: VlasovPoissonSolver. + diff --git a/src/geometryRTheta/poisson/ipoissonsolver.hpp b/src/geometryRTheta/poisson/ipoissonsolver.hpp index 77a977000..44d6367e4 100644 --- a/src/geometryRTheta/poisson/ipoissonsolver.hpp +++ b/src/geometryRTheta/poisson/ipoissonsolver.hpp @@ -6,13 +6,27 @@ #include +/** + * @brief Base class for Poisson solver. + */ class IPoissonSolver { public: virtual ~IPoissonSolver() = default; + /** + * @brief Compute the electrical potential and + * the electric field from the Poisson equation. + * + * @param[out] electrostatic_potential + * The solution of the Poisson equation. + * @param[out] electric_field + * The electric field @f$E = -\nabla \phi@f$. + * @param[in] allfdistribu + * The rhs of the Poisson equation. + */ virtual void operator()( DSpanRP electrostatic_potential, - DSpanRP electric_field, + VectorFieldSpan> electric_field, DViewRP allfdistribu) const = 0; }; diff --git a/src/geometryRTheta/poisson/nullpoissonsolver.cpp b/src/geometryRTheta/poisson/nullpoissonsolver.cpp deleted file mode 100644 index 99dc19bca..000000000 --- a/src/geometryRTheta/poisson/nullpoissonsolver.cpp +++ /dev/null @@ -1,5 +0,0 @@ -// SPDX-License-Identifier: MIT - -#include "nullpoissonsolver.hpp" - -void NullPoissonSolver::operator()(DSpanRP const, DSpanRP const, DViewRP const) const {} diff --git a/src/geometryRTheta/poisson/nullpoissonsolver.hpp b/src/geometryRTheta/poisson/nullpoissonsolver.hpp deleted file mode 100644 index fe9f5dc64..000000000 --- a/src/geometryRTheta/poisson/nullpoissonsolver.hpp +++ /dev/null @@ -1,18 +0,0 @@ -// SPDX-License-Identifier: MIT - -#pragma once - -#include - -#include "ipoissonsolver.hpp" - -class NullPoissonSolver : public IPoissonSolver -{ -public: - NullPoissonSolver() = default; - - ~NullPoissonSolver() override = default; - - void operator()(DSpanRP electrostatic_potential, DSpanRP electric_field, DViewRP allfdistribu) - const override; -}; diff --git a/src/geometryRTheta/poisson/poisson_rhs_function.hpp b/src/geometryRTheta/poisson/poisson_rhs_function.hpp new file mode 100644 index 000000000..d620a0b6f --- /dev/null +++ b/src/geometryRTheta/poisson/poisson_rhs_function.hpp @@ -0,0 +1,50 @@ +#pragma once + +#include + +#include + +#include "geometry.hpp" +#include "spline_interpolator_2d_rp.hpp" + + + +/** + * @brief Type of right-hand side (rhs) function of the Poisson equation. + */ +class PoissonRHSFunction +{ +private: + Spline2DView const m_coefs; + SplineEvaluator2D const& m_evaluator; + +public: + /** + * @brief Instantiate a PoissonRHSFunction. + * + * @param[in] coefs + * The bsplines coefficients of the right-hand side function. + * @param[in] evaluator + * Evaluator on bsplines. + */ + PoissonRHSFunction(Spline2DView coefs, SplineEvaluator2D const& evaluator) + : m_coefs(coefs) + , m_evaluator(evaluator) + { + } + + ~PoissonRHSFunction() {}; + + /** + * @brief Get the value of the function at a given coordinate. + * + * @param[in] coord_rp + * Polar coordinate where we want to evaluate the rhs function. + * + * @return A double with the value of the rhs at the given coordinate. + */ + double operator()(CoordRP const& coord_rp) const + { + return m_evaluator(coord_rp, m_coefs); + } +}; diff --git a/src/geometryRTheta/poisson/polarpoissonsolver.hpp b/src/geometryRTheta/poisson/polarpoissonsolver.hpp index 165e4983a..13dafa486 100644 --- a/src/geometryRTheta/poisson/polarpoissonsolver.hpp +++ b/src/geometryRTheta/poisson/polarpoissonsolver.hpp @@ -12,23 +12,24 @@ #include -//#include "chargedensitycalculator.hpp" -//#include "ipoissonsolver.hpp" - -template +#include "geometry.hpp" + +/** + * @brief Define a polar Poisson solver. + * + * Solve the following Poisson equation + * + * (1) @f$ L\phi = - \nabla \cdot (\alpha \nabla \phi) + \beta \phi = \rho @f$, in @f$ \Omega@f$, + * + * @f$ \phi = 0 @f$, on @f$ \partial \Omega@f$, + * + * (see in Emily Bourne's thesis "Non-Uniform Numerical Schemes for the Modelling of Turbulence + * in the 5D GYSELA Code". December 2022.) + * + */ class PolarSplineFEMPoissonSolver { public: - using BSplinesR = typename PolarBSplines::BSplinesR_tag; - using BSplinesP = typename PolarBSplines::BSplinesP_tag; - - using RDimR = typename BSplinesR::tag_type; - using RDimP = typename BSplinesP::tag_type; - -public: - using SplineRP = ddc::ChunkSpan>; - using SplinePolar = PolarSpline; - struct RBasisSubset { }; @@ -42,6 +43,12 @@ class PolarSplineFEMPoissonSolver { }; + /** + * @brief This struct is used to provide alternative tags. + * + * It is used to create new meshes. This is necessary because + * in DDC we can only create a single mesh for a tagged dimension. + */ template struct InternalTagGenerator { @@ -49,24 +56,68 @@ class PolarSplineFEMPoissonSolver }; private: + /** + * @brief Tag the first dimension for the quadrature. + */ using QDimR = InternalTagGenerator; + /** + * @brief Tag the second dimension for the quadrature. + */ using QDimP = InternalTagGenerator; + /** + * @brief Tag the first dimension for the quadrature mesh. + */ using QDimRMesh = ddc::NonUniformPointSampling; + /** + * @brief Tag the second dimension for the quadrature mesh. + */ using QDimPMesh = ddc::NonUniformPointSampling; + /** + * @brief Tag the quadrature domain in the first dimension. + */ using QuadratureDomainR = ddc::DiscreteDomain; + /** + * @brief Tag the quadrature domain in the second dimension. + */ using QuadratureDomainP = ddc::DiscreteDomain; + /** + * @brief Tag the quadrature domain. + */ using QuadratureDomainRP = ddc::DiscreteDomain; - using QuadratureMeshR = ddc::DiscreteElement; - using QuadratureMeshP = ddc::DiscreteElement; - using QuadratureMeshRP = ddc::DiscreteElement; + /** + * @brief Tag the elements (index) of the quadrature domain in the first dimension. + */ + using QuadratureIndexR = ddc::DiscreteElement; + /** + * @brief Tag the elements (index) of the quadrature domain in the second dimension. + */ + using QuadratureIndexP = ddc::DiscreteElement; + /** + * @brief Tag the elements (index) of the quadrature domain. + */ + using QuadratureIndexRP = ddc::DiscreteElement; + /** + * @brief Tag a vector on the first dimension of the quadrature mesh. + */ using QuadratureLengthR = ddc::DiscreteVector; + /** + * @brief Tag a vector on the second dimension of the quadrature mesh. + */ using QuadratureLengthP = ddc::DiscreteVector; + /** + * @brief Object storing a value and a value of the derivative + * of a 1D function. + */ struct EvalDeriv1DType { double value; double derivative; }; + /** + * @brief Object storing a value and a value of the derivatives + * in each direction of a 2D function. + */ struct EvalDeriv2DType { double value; @@ -74,18 +125,20 @@ class PolarSplineFEMPoissonSolver double poloidal_derivative; }; - using BSpline2DDomain = ddc::DiscreteDomain; - using IDimBSpline2D = ddc::DiscreteElement; - using IDimPolarBspl = ddc::DiscreteElement; - + /** + * @brief Tag an index of cell. + */ using CellIndex = ddc::DiscreteElement; + /** + * @brief Tag type of matrix element. + */ using MatrixElement = Eigen::Triplet; private: static constexpr int n_gauss_legendre_r = BSplinesR::degree() + 1; static constexpr int n_gauss_legendre_p = BSplinesP::degree() + 1; - static constexpr int n_overlap_cells = PolarBSplines::continuity + 1; + static constexpr int n_overlap_cells = PolarBSplinesRP::continuity + 1; static constexpr ddc::DiscreteVector n_non_zero_bases_r = ddc::DiscreteVector(BSplinesR::degree() + 1); @@ -101,9 +154,9 @@ class PolarSplineFEMPoissonSolver const int nbasis_p; // Domains - ddc::DiscreteDomain fem_non_singular_domain; - ddc::DiscreteDomain radial_bsplines; - ddc::DiscreteDomain polar_bsplines; + BSDomainPolar fem_non_singular_domain; + BSDomainR radial_bsplines; + BSDomainP polar_bsplines; QuadratureDomainR quadrature_domain_r; QuadratureDomainP quadrature_domain_p; @@ -116,7 +169,7 @@ class PolarSplineFEMPoissonSolver ddc::Chunk weights_p; // Basis Spline values and derivatives at Gauss-Legendre points - ddc::Chunk> + ddc::Chunk> singular_basis_vals_and_derivs; ddc::Chunk> r_basis_vals_and_derivs; @@ -125,21 +178,42 @@ class PolarSplineFEMPoissonSolver ddc::Chunk int_volume; - PolarSplineEvaluator m_polar_spline_evaluator; + PolarSplineEvaluator m_polar_spline_evaluator; Eigen::SimplicialLDLT> m_matrix; public: + /** + * @brief Instantiate a polar Poisson solver using FEM with B-splines. + * + * The equation we are studying: + * + * (1) @f$ L\phi = - \nabla \cdot (\alpha \nabla \phi) + \beta \phi = \rho @f$, in @f$ \Omega@f$, + * + * @f$ \phi = 0 @f$, on @f$ \partial \Omega@f$. + * + * @param[in] coeff_alpha + * The spline representation of the @f$ \alpha @f$ function in the + * definition of the Poisson equation. + * @param[in] coeff_beta + * The spline representation of the @f$ \beta @f$ function in the + * definition of the Poisson equation. + * @param[in] mapping + * The mapping from the logical domain to the physical domain where + * the equation is defined. + * + * @tparam Mapping A Curvilinear2DToCartesian class. + */ template PolarSplineFEMPoissonSolver( - ddc::ChunkSpan> coeff_alpha, - ddc::ChunkSpan> coeff_beta, + Spline2DView coeff_alpha, + Spline2DView coeff_beta, Mapping const& mapping) : nbasis_r(ddc::discrete_space().nbasis() - n_overlap_cells - 1) , nbasis_p(ddc::discrete_space().nbasis()) , fem_non_singular_domain( - ddc::discrete_space().non_singular_domain().remove_last( - ddc::DiscreteVector {nbasis_p})) + ddc::discrete_space().non_singular_domain().remove_last( + ddc::DiscreteVector {nbasis_p})) , radial_bsplines(ddc::discrete_space().full_domain().remove_first( ddc::DiscreteVector {n_overlap_cells})) , polar_bsplines(ddc::discrete_space().full_domain().take_first( @@ -160,8 +234,8 @@ class PolarSplineFEMPoissonSolver , points_p(quadrature_domain_p) , weights_r(quadrature_domain_r) , weights_p(quadrature_domain_p) - , singular_basis_vals_and_derivs(ddc::DiscreteDomain( - PolarBSplines::singular_domain(), + , singular_basis_vals_and_derivs(ddc::DiscreteDomain( + PolarBSplinesRP::singular_domain(), ddc::select(quadrature_domain_singular), ddc::select(quadrature_domain_singular))) , r_basis_vals_and_derivs(ddc::DiscreteDomain< @@ -171,7 +245,7 @@ class PolarSplineFEMPoissonSolver PBasisSubset, QDimPMesh>(non_zero_bases_p, quadrature_domain_p)) , int_volume(QuadratureDomainRP(quadrature_domain_r, quadrature_domain_p)) - , m_polar_spline_evaluator(g_polar_null_boundary_2d) + , m_polar_spline_evaluator(g_polar_null_boundary_2d) { const std::size_t ncells_r = ddc::discrete_space().ncells(); const std::size_t ncells_p = ddc::discrete_space().ncells(); @@ -221,7 +295,7 @@ class PolarSplineFEMPoissonSolver ddc::init_discrete_space(vect_points_p); // Find value and derivative of 1D bsplines in radial direction - ddc::for_each(quadrature_domain_r, [&](QuadratureMeshR const ir) { + ddc::for_each(quadrature_domain_r, [&](QuadratureIndexR const ir) { std::array data; DSpan2D vals(data.data(), n_non_zero_bases_r, 2); ddc::discrete_space().eval_basis_and_n_derivs(vals, get_coordinate(ir), 1); @@ -232,7 +306,7 @@ class PolarSplineFEMPoissonSolver }); // Find value and derivative of 1D bsplines in poloidal direction - ddc::for_each(quadrature_domain_p, [&](QuadratureMeshP const ip) { + ddc::for_each(quadrature_domain_p, [&](QuadratureIndexP const ip) { std::array data; DSpan2D vals(data.data(), n_non_zero_bases_p, 2); ddc::discrete_space().eval_basis_and_n_derivs(vals, get_coordinate(ip), 1); @@ -242,35 +316,35 @@ class PolarSplineFEMPoissonSolver } }); - auto singular_domain = PolarBSplines::singular_domain(); + auto singular_domain = PolarBSplinesRP::singular_domain(); // Find value and derivative of 2D bsplines covering the singular point - ddc::for_each(quadrature_domain_singular, [&](QuadratureMeshRP const irp) { - std::array singular_data; + ddc::for_each(quadrature_domain_singular, [&](QuadratureIndexRP const irp) { + std::array singular_data; std::array data; - DSpan1D singular_vals(singular_data.data(), PolarBSplines::n_singular_basis()); + DSpan1D singular_vals(singular_data.data(), PolarBSplinesRP::n_singular_basis()); DSpan2D vals(data.data(), n_non_zero_bases_r, n_non_zero_bases_p); - QuadratureMeshR ir = ddc::select(irp); - QuadratureMeshP ip = ddc::select(irp); + QuadratureIndexR ir = ddc::select(irp); + QuadratureIndexP ip = ddc::select(irp); - const ddc::Coordinate coord(get_coordinate(irp)); + const CoordRP coord(get_coordinate(irp)); // Calculate the value - ddc::discrete_space().eval_basis(singular_vals, vals, coord); + ddc::discrete_space().eval_basis(singular_vals, vals, coord); for (auto ib : singular_domain) { singular_basis_vals_and_derivs(ib, ir, ip).value = singular_vals[ib.uid()]; } // Calculate the radial derivative - ddc::discrete_space().eval_deriv_r(singular_vals, vals, coord); + ddc::discrete_space().eval_deriv_r(singular_vals, vals, coord); for (auto ib : singular_domain) { singular_basis_vals_and_derivs(ib, ir, ip).radial_derivative = singular_vals[ib.uid()]; } // Calculate the poloidal derivative - ddc::discrete_space().eval_deriv_p(singular_vals, vals, coord); + ddc::discrete_space().eval_deriv_p(singular_vals, vals, coord); for (auto ib : singular_domain) { singular_basis_vals_and_derivs(ib, ir, ip).poloidal_derivative = singular_vals[ib.uid()]; @@ -279,10 +353,10 @@ class PolarSplineFEMPoissonSolver // Find the integral volume associated with each point used in the quadrature scheme QuadratureDomainRP all_quad_points(quadrature_domain_r, quadrature_domain_p); - ddc::for_each(all_quad_points, [&](QuadratureMeshRP const irp) { - QuadratureMeshR const ir = ddc::select(irp); - QuadratureMeshP const ip = ddc::select(irp); - ddc::Coordinate coord(get_coordinate(ir), get_coordinate(ip)); + ddc::for_each(all_quad_points, [&](QuadratureIndexRP const irp) { + QuadratureIndexR const ir = ddc::select(irp); + QuadratureIndexP const ip = ddc::select(irp); + CoordRP coord(get_coordinate(ir), get_coordinate(ip)); int_volume(ir, ip) = abs(mapping.jacobian(coord)) * weights_r(ir) * weights_p(ip); }); @@ -293,16 +367,16 @@ class PolarSplineFEMPoissonSolver g_null_boundary_2d); constexpr int n_elements_singular - = PolarBSplines::n_singular_basis() * PolarBSplines::n_singular_basis(); + = PolarBSplinesRP::n_singular_basis() * PolarBSplinesRP::n_singular_basis(); const int n_elements_overlap - = 2 * (PolarBSplines::n_singular_basis() * BSplinesR::degree() * nbasis_p); + = 2 * (PolarBSplinesRP::n_singular_basis() * BSplinesR::degree() * nbasis_p); const int n_stencil_p = nbasis_p * min(int(1 + 2 * BSplinesP::degree()), nbasis_p); const int n_stencil_r = nbasis_r * (1 + 2 * BSplinesR::degree()) - (1 + BSplinesR::degree()) * BSplinesR::degree(); const int n_elements_stencil = n_stencil_r * n_stencil_p; // Matrix size is equal to the number Polar bspline - const int n_matrix_size = ddc::discrete_space().nbasis() - nbasis_p; + const int n_matrix_size = ddc::discrete_space().nbasis() - nbasis_p; Eigen::SparseMatrix matrix(n_matrix_size, n_matrix_size); const int n_matrix_elements = n_elements_singular + n_elements_overlap + n_elements_stencil; std::vector matrix_elements(n_matrix_elements); @@ -318,9 +392,9 @@ class PolarSplineFEMPoissonSolver quadrature_domain_singular, 0.0, ddc::reducer::sum(), - [&](QuadratureMeshRP const quad_idx) { - QuadratureMeshR const ir = ddc::select(quad_idx); - QuadratureMeshP const ip = ddc::select(quad_idx); + [&](QuadratureIndexRP const quad_idx) { + QuadratureIndexR const ir = ddc::select(quad_idx); + QuadratureIndexP const ip = ddc::select(quad_idx); std::optional> null(std::nullopt); return weak_integral_element( ir, @@ -337,11 +411,10 @@ class PolarSplineFEMPoissonSolver assert(matrix_idx == n_elements_singular); // Create domains associated with the 2D splines - ddc::DiscreteDomain central_radial_bspline_domain( + BSDomainR central_radial_bspline_domain( radial_bsplines.take_first(ddc::DiscreteVector {BSplinesR::degree()})); - BSpline2DDomain - non_singular_domain_near_centre(central_radial_bspline_domain, polar_bsplines); + BSDomainRP non_singular_domain_near_centre(central_radial_bspline_domain, polar_bsplines); const ddc::DiscreteDomain r_cells_near_centre( ddc::DiscreteElement {0}, @@ -350,7 +423,7 @@ class PolarSplineFEMPoissonSolver // Calculate the matrix elements where bspline products overlap the bsplines which cover the singular point ddc::for_each(singular_domain, [&](IDimPolarBspl const idx_test) { ddc::for_each(non_singular_domain_near_centre, [&](IDimBSpline2D const idx_trial) { - const IDimPolarBspl polar_idx_trial(PolarBSplines::get_polar_index(idx_trial)); + const IDimPolarBspl polar_idx_trial(PolarBSplinesRP::get_polar_index(idx_trial)); const ddc::DiscreteElement r_idx_trial( ddc::select(idx_trial)); const ddc::DiscreteElement p_idx_trial( @@ -391,9 +464,9 @@ class PolarSplineFEMPoissonSolver cell_quad_points, 0.0, ddc::reducer::sum(), - [&](QuadratureMeshRP const quad_idx) { - QuadratureMeshR const ir = ddc::select(quad_idx); - QuadratureMeshP const ip = ddc::select(quad_idx); + [&](QuadratureIndexRP const quad_idx) { + QuadratureIndexR const ir = ddc::select(quad_idx); + QuadratureIndexP const ip = ddc::select(quad_idx); return weak_integral_element( ir, ip, @@ -417,17 +490,17 @@ class PolarSplineFEMPoissonSolver // Calculate the matrix elements following a stencil ddc::for_each(fem_non_singular_domain, [&](IDimPolarBspl const polar_idx_test) { - const IDimBSpline2D idx_test(PolarBSplines::get_2d_index(polar_idx_test)); + const IDimBSpline2D idx_test(PolarBSplinesRP::get_2d_index(polar_idx_test)); const std::size_t r_idx_test(ddc::select(idx_test).uid()); const std::size_t p_idx_test(ddc::select(idx_test).uid()); // Calculate the index of the elements that are already filled - ddc::DiscreteDomain remaining_p( + BSDomainP remaining_p( ddc::DiscreteElement {p_idx_test}, ddc::DiscreteVector {BSplinesP::degree() + 1}); ddc::for_each(remaining_p, [&](auto const p_idx_trial) { IDimBSpline2D idx_trial(ddc::DiscreteElement(r_idx_test), p_idx_trial); - IDimPolarBspl polar_idx_trial(PolarBSplines::get_polar_index( + IDimPolarBspl polar_idx_trial(PolarBSplinesRP::get_polar_index( IDimBSpline2D(r_idx_test, pmod(p_idx_trial.uid())))); double element = get_matrix_stencil_element( idx_test, @@ -447,23 +520,23 @@ class PolarSplineFEMPoissonSolver = MatrixElement(polar_idx_trial.uid(), polar_idx_test.uid(), element); } }); - ddc::DiscreteDomain remaining_r( + BSDomainR remaining_r( ddc::select(idx_test) + 1, ddc::DiscreteVector { min(BSplinesR::degree(), ddc::discrete_space().nbasis() - 2 - r_idx_test)}); - ddc::DiscreteDomain relevant_p( + BSDomainP relevant_p( ddc::DiscreteElement { p_idx_test + ddc::discrete_space().nbasis() - BSplinesP::degree()}, ddc::DiscreteVector {2 * BSplinesP::degree() + 1}); - BSpline2DDomain trial_domain(remaining_r, relevant_p); + BSDomainRP trial_domain(remaining_r, relevant_p); ddc::for_each(trial_domain, [&](IDimBSpline2D const idx_trial) { const int r_idx_trial(ddc::select(idx_trial).uid()); const int p_idx_trial(ddc::select(idx_trial).uid()); - IDimPolarBspl polar_idx_trial(PolarBSplines::get_polar_index( + IDimPolarBspl polar_idx_trial(PolarBSplinesRP::get_polar_index( IDimBSpline2D(r_idx_trial, pmod(p_idx_trial)))); double element = get_matrix_stencil_element( idx_test, @@ -488,33 +561,44 @@ class PolarSplineFEMPoissonSolver m_matrix.compute(matrix); } - template - void operator()( - RHSFunction const& rhs, - ddc::ChunkSpan const, Domain> const coords_eval, - ddc::ChunkSpan result) const + + /** + * @brief Solve the Poisson equation. + * + * This operator returns the coefficients associated with the B-Splines + * of the solution @f$\phi@f$. + * + * @param[in] rhs + * The rhs @f$ \rho@f$ of the Poisson equation. + * The type is templated but we can use the PoissonRHSFunction + * class. + * @param[out] spline + * The spline representation of the solution @f$\phi@f$. + */ + template + void operator()(RHSFunction const& rhs, SplinePolar& spline) const { Eigen::VectorXd b( - ddc::discrete_space().nbasis() + ddc::discrete_space().nbasis() - ddc::discrete_space().nbasis()); // Fill b - ddc::for_each(PolarBSplines::singular_domain(), [&](IDimPolarBspl const idx) { + ddc::for_each(PolarBSplinesRP::singular_domain(), [&](IDimPolarBspl const idx) { b(idx.uid()) = ddc::transform_reduce( quadrature_domain_singular, 0.0, ddc::reducer::sum(), - [&](QuadratureMeshRP const quad_idx) { - QuadratureMeshR const ir = ddc::select(quad_idx); - QuadratureMeshP const ip = ddc::select(quad_idx); - ddc::Coordinate coord(get_coordinate(ir), get_coordinate(ip)); + [&](QuadratureIndexRP const quad_idx) { + QuadratureIndexR const ir = ddc::select(quad_idx); + QuadratureIndexP const ip = ddc::select(quad_idx); + CoordRP coord(get_coordinate(ir), get_coordinate(ip)); return rhs(coord) * singular_basis_vals_and_derivs(idx, ir, ip).value * int_volume(ir, ip); }); }); const std::size_t ncells_r = ddc::discrete_space().ncells(); ddc::for_each(fem_non_singular_domain, [&](IDimPolarBspl const idx) { - const IDimBSpline2D idx_2d(PolarBSplines::get_2d_index(idx)); + const IDimBSpline2D idx_2d(PolarBSplinesRP::get_2d_index(idx)); const std::size_t r_idx(ddc::select(idx_2d).uid()); const std::size_t p_idx(ddc::select(idx_2d).uid()); @@ -553,11 +637,10 @@ class PolarSplineFEMPoissonSolver cell_quad_points, 0.0, ddc::reducer::sum(), - [&](QuadratureMeshRP const quad_idx) { - QuadratureMeshR const ir = ddc::select(quad_idx); - QuadratureMeshP const ip = ddc::select(quad_idx); - ddc::Coordinate - coord(get_coordinate(ir), get_coordinate(ip)); + [&](QuadratureIndexRP const quad_idx) { + QuadratureIndexR const ir = ddc::select(quad_idx); + QuadratureIndexP const ip = ddc::select(quad_idx); + CoordRP coord(get_coordinate(ir), get_coordinate(ip)); double rb = r_basis_vals_and_derivs(ib_r, ir).value; double pb = p_basis_vals_and_derivs(ib_p, ip).value; return rhs(coord) * rb * pb * int_volume(ir, ip); @@ -569,23 +652,21 @@ class PolarSplineFEMPoissonSolver // Solve the matrix equation Eigen::VectorXd x = m_matrix.solve(b); - ddc::DiscreteDomain non_singular_2d_domain( + BSDomainRP non_singular_2d_domain( radial_bsplines.remove_last(ddc::DiscreteVector {1}), polar_bsplines); - ddc::DiscreteDomain dirichlet_boundary_domain( + BSDomainRP dirichlet_boundary_domain( radial_bsplines.take_last(ddc::DiscreteVector {1}), polar_bsplines); - ddc::DiscreteDomain polar_domain(ddc::discrete_space().full_domain()); - SplinePolar - spline(PolarBSplines::singular_domain(), - ddc::DiscreteDomain(radial_bsplines, polar_domain)); + BSDomainP polar_domain(ddc::discrete_space().full_domain()); + // Fill the spline - ddc::for_each(PolarBSplines::singular_domain(), [&](IDimPolarBspl const idx) { + ddc::for_each(PolarBSplinesRP::singular_domain(), [&](IDimPolarBspl const idx) { spline.singular_spline_coef(idx) = x(idx.uid()); }); ddc::for_each(fem_non_singular_domain, [&](IDimPolarBspl const idx) { - const IDimBSpline2D idx_2d(PolarBSplines::get_2d_index(idx)); + const IDimBSpline2D idx_2d(PolarBSplinesRP::get_2d_index(idx)); spline.spline_coef(idx_2d) = x(idx.uid()); }); ddc::for_each(dirichlet_boundary_domain, [&](IDimBSpline2D const idx) { @@ -593,7 +674,7 @@ class PolarSplineFEMPoissonSolver }); // Copy the periodic elements - ddc::DiscreteDomain copy_domain( + BSDomainRP copy_domain( radial_bsplines, polar_domain.remove_first( ddc::DiscreteVector(ddc::discrete_space().nbasis()))); @@ -604,15 +685,42 @@ class PolarSplineFEMPoissonSolver ddc::select(idx_2d) - ddc::discrete_space().nbasis()); }); + } + + /** + * @brief Solve the Poisson equation. + * + * This operator uses the other operator () and returns the values on + * the grid of the solution @f$\phi@f$. + * + * @param[in] rhs + * The rhs @f$ \rho@f$ of the Poisson equation. + * The type is templated but we can use the PoissonRHSFunction + * class. + * @param[in] coords_eval + * A Chunk of coordinates where we want to compute the solution. + * @param[out] result + * The values of the solution @f$\phi@f$ on the given coords_eval. + */ + template + void operator()(RHSFunction const& rhs, ViewRP const coords_eval, DSpanRP result) const + { + BSDomainP polar_domain(ddc::discrete_space().full_domain()); + SplinePolar + spline(PolarBSplinesRP::singular_domain(), + BSDomainRP(radial_bsplines, polar_domain)); + + (*this)(rhs, spline); m_polar_spline_evaluator(result, coords_eval, spline); } + private: static QuadratureDomainRP get_quadrature_points_in_cell(int cell_idx_r, int cell_idx_p) { - const QuadratureMeshR first_quad_point_r(cell_idx_r * n_gauss_legendre_r); - const QuadratureMeshP first_quad_point_p(cell_idx_p * n_gauss_legendre_p); + const QuadratureIndexR first_quad_point_r(cell_idx_r * n_gauss_legendre_r); + const QuadratureIndexP first_quad_point_p(cell_idx_p * n_gauss_legendre_p); constexpr QuadratureLengthR n_GL_r(n_gauss_legendre_r); constexpr QuadratureLengthP n_GL_p(n_gauss_legendre_p); const QuadratureDomainR quad_points_r(first_quad_point_r, n_GL_r); @@ -622,12 +730,12 @@ class PolarSplineFEMPoissonSolver template double weak_integral_element( - QuadratureMeshR ir, - QuadratureMeshP ip, + QuadratureIndexR ir, + QuadratureIndexP ip, EvalDeriv2DType const& test_bspline_val_and_deriv, EvalDeriv2DType const& trial_bspline_val_and_deriv, - ddc::ChunkSpan> coeff_alpha, - ddc::ChunkSpan> coeff_beta, + Spline2DView coeff_alpha, + Spline2DView coeff_beta, SplineEvaluator2D const& evaluator, Mapping const& mapping) { @@ -646,13 +754,13 @@ class PolarSplineFEMPoissonSolver template double weak_integral_element( - QuadratureMeshR ir, - QuadratureMeshP ip, + QuadratureIndexR ir, + QuadratureIndexP ip, EvalDeriv2DType const& test_bspline_val_and_deriv, EvalDeriv1DType const& trial_bspline_val_and_deriv_r, EvalDeriv1DType const& trial_bspline_val_and_deriv_p, - ddc::ChunkSpan> coeff_alpha, - ddc::ChunkSpan> coeff_beta, + Spline2DView coeff_alpha, + Spline2DView coeff_beta, SplineEvaluator2D const& evaluator, Mapping const& mapping) { @@ -671,13 +779,13 @@ class PolarSplineFEMPoissonSolver template double weak_integral_element( - QuadratureMeshR ir, - QuadratureMeshP ip, + QuadratureIndexR ir, + QuadratureIndexP ip, EvalDeriv1DType const& test_bspline_val_and_deriv_r, EvalDeriv2DType const& trial_bspline_val_and_deriv, EvalDeriv1DType const& test_bspline_val_and_deriv_p, - ddc::ChunkSpan> coeff_alpha, - ddc::ChunkSpan> coeff_beta, + Spline2DView coeff_alpha, + Spline2DView coeff_beta, SplineEvaluator2D const& evaluator, Mapping const& mapping) { @@ -696,14 +804,14 @@ class PolarSplineFEMPoissonSolver template double weak_integral_element( - QuadratureMeshR ir, - QuadratureMeshP ip, + QuadratureIndexR ir, + QuadratureIndexP ip, EvalDeriv1DType const& test_bspline_val_and_deriv_r, EvalDeriv1DType const& trial_bspline_val_and_deriv_r, EvalDeriv1DType const& test_bspline_val_and_deriv_p, EvalDeriv1DType const& trial_bspline_val_and_deriv_p, - ddc::ChunkSpan> coeff_alpha, - ddc::ChunkSpan> coeff_beta, + Spline2DView coeff_alpha, + Spline2DView coeff_beta, SplineEvaluator2D const& evaluator, Mapping const& mapping) { @@ -742,14 +850,14 @@ class PolarSplineFEMPoissonSolver template double templated_weak_integral_element( - QuadratureMeshR ir, - QuadratureMeshP ip, + QuadratureIndexR ir, + QuadratureIndexP ip, TestValDerivType const& test_bspline_val_and_deriv, TrialValDerivType const& trial_bspline_val_and_deriv, TestValDerivType const& test_bspline_val_and_deriv_p, TrialValDerivType const& trial_bspline_val_and_deriv_p, - ddc::ChunkSpan> coeff_alpha, - ddc::ChunkSpan> coeff_beta, + Spline2DView coeff_alpha, + Spline2DView coeff_beta, SplineEvaluator2D const& spline_evaluator, Mapping const& mapping) { @@ -763,7 +871,7 @@ class PolarSplineFEMPoissonSolver EvalDeriv1DType> || std::is_same_v); // Calculate coefficients at quadrature point - ddc::Coordinate coord(get_coordinate(ir), get_coordinate(ip)); + CoordRP coord(get_coordinate(ir), get_coordinate(ip)); const double alpha = spline_evaluator(coord, coeff_alpha); const double beta = spline_evaluator(coord, coeff_beta); @@ -796,8 +904,8 @@ class PolarSplineFEMPoissonSolver double get_matrix_stencil_element( IDimBSpline2D idx_test, IDimBSpline2D idx_trial, - ddc::ChunkSpan> coeff_alpha, - ddc::ChunkSpan> coeff_beta, + Spline2DView coeff_alpha, + Spline2DView coeff_beta, SplineEvaluator2D const& evaluator, Mapping const& mapping) { @@ -889,9 +997,9 @@ class PolarSplineFEMPoissonSolver cell_quad_points, 0.0, ddc::reducer::sum(), - [&](QuadratureMeshRP const quad_idx) { - QuadratureMeshR const r_idx = ddc::select(quad_idx); - QuadratureMeshP const p_idx = ddc::select(quad_idx); + [&](QuadratureIndexRP const quad_idx) { + QuadratureIndexR const r_idx = ddc::select(quad_idx); + QuadratureIndexP const p_idx = ddc::select(quad_idx); return weak_integral_element( r_idx, p_idx, diff --git a/src/geometryRTheta/poisson/vlasovpoissonsolver.hpp b/src/geometryRTheta/poisson/vlasovpoissonsolver.hpp new file mode 100644 index 000000000..a4f1cdd23 --- /dev/null +++ b/src/geometryRTheta/poisson/vlasovpoissonsolver.hpp @@ -0,0 +1,296 @@ +#pragma once + + +#include + +#include +#include + +#include +#include +#include +#include + +#include "ipoissonsolver.hpp" +#include "poisson_rhs_function.hpp" +#include "polarpoissonsolver.hpp" + + + +/** + * @brief Solve the Poisson equation and return the electric + * field for the coupled Vlasov equation. + * + * The Vlasov-Poisson equations are given by + * + * - (1) @f$ \partial_t \rho - E_y \partial_x \rho + E_x \partial_y\rho = 0 @f$, + * + * - (2) @f$ - \Delta \phi = \rho @f$, + * + * - (3) and @f$ E = -\nabla \phi @f$. + * + * The functions are defined on a logical domain, and the mapping from the logical + * domain to the physical domain is written @f$\mathcal{F}@f$. + * + * We here focus on equations (2) and (3). We first compute @f$ \phi @f$ + * on B-splines with the given Poisson solver. Then in the VlasovPoissonSolver::operator() + * we compute the electric field thanks to (3) using the B-splines coefficients. + * Depending on the given mapping, the computation at the center point is not + * always well-defined so we linearize around the center point as explained + * in Edoardo Zoni's article (https://doi.org/10.1016/j.jcp.2019.108889). + * + * With (3), we compute the electric field + * + * - @f$ E(r, \theta) = -\nabla_{x,y} \phi(r, \theta) @f$, + * + * with @f$\nabla_{x,y} \phi(r, \theta) @f$ from + * @f$ \nabla_{x,y} \phi(r, \theta) + * = (J_{\mathcal{F}}^{-1})^T \nabla_{r, \theta} \phi(r, \theta) @f$, + * + * and + * @f$ \nabla \phi = \sum_{ij} \partial_{x_i} \phi g^{i,j} e_j @f$ + * with @f$ g^{i,j} @f$ the coefficients of the inverse metric tensor + * and @f$ e_j @f$ the unnormalized local covariant base: + * + * + * - @f$ \nabla \phi(r,\theta) \cdot e_r = \partial_r \phi(r,\theta) g^{1,1} + \partial_\theta \phi(r,\theta)g^{2,1} @f$, + * + * - @f$ \nabla \phi(r,\theta) \cdot e_\theta = \partial_r \phi(r,\theta) g^{1,2} + \partial_\theta \phi(r,\theta)g^{2,2} @f$. + * + * + * + * For @f$ r < \varepsilon @f$, + * @f$ E(r, \theta) = \left( 1 - \frac{r}{\varepsilon} \right) E(0, \theta) + * + \frac{r}{\varepsilon} E(\varepsilon, \theta) @f$, + * + * with @f$ E(0, \theta) @f$ computed thanks to + * + * - @f$ \partial_r \phi (0, \theta_1) = \left[\partial_r x \partial_x \phi + * + \partial_r y \partial_y \phi \right](0, \theta_1) @f$, + * + * - @f$ \partial_r \phi (0, \theta_2) = \left[\partial_r x \partial_x \phi + * + \partial_r y \partial_y \phi \right] (0, \theta_2) @f$, + * + * where @f$ \theta_1 @f$ and @f$ \theta_2 @f$ correspond to + * linearly independent directions. + * + * + * The equation (1) is computed thanks to advection operator. + * + * + * @tparam Mapping A Curvilinear2DToCartesian class. + * + * + * @see PolarSplineFEMPoissonSolver + * + */ +template +class VlasovPoissonSolver : public IPoissonSolver +{ +public: + /** + * @brief Define a 2x2 matrix with an 2D array of an 2D array. + */ + using Matrix_2x2 = std::array, 2>; + +private: + Mapping const& m_mapping; + + SplineRPBuilder const& m_rhs_builder; + SplineRPEvaluator const& m_rhs_evaluator; + + PolarSplineEvaluator m_polar_spline_evaluator; + + PolarSplineFEMPoissonSolver const& m_poisson_solver; + + double const m_epsilon; + + + static constexpr int n_overlap_cells = PolarBSplinesRP::continuity + 1; + +public: + /** + * @brief Instantiate a VlasovPoissonSolver. + * + * @param[in] mapping + * The mapping @f$ \mathcal{F} @f$ from the logical domain to the physical domain. + * @param[in] builder + * The splines builder for the rhs. + * @param[in] evaluator + * The splines evaluator for the rhs. + * @param[in] poisson_solver + * The Poisson solver used to solve (2). + * @param[in] epsilon + * The parameter @f$ \varepsilon @f$ for the linearization of the + * electric field. + */ + VlasovPoissonSolver( + Mapping const& mapping, + SplineRPBuilder const& builder, + SplineRPEvaluator const& evaluator, + PolarSplineFEMPoissonSolver const& poisson_solver, + double const epsilon = 1e-12) + : m_mapping(mapping) + , m_rhs_builder(builder) + , m_rhs_evaluator(evaluator) + , m_polar_spline_evaluator(g_polar_null_boundary_2d) + , m_poisson_solver(poisson_solver) + , m_epsilon(epsilon) {}; + + ~VlasovPoissonSolver() {}; + + + /** + * @brief Compute the electric field from the Poisson equation. + * + * @param[out] electrostatic_potential + * The solution @f$\phi@f$ of the Poisson equation (2). + * @param[out] electric_field + * The electric field @f$E = -\nabla \phi@f$. + * @param[in] allfdistribu + * The rhs @f$\rho@f$ of the Poisson equation (2). + */ + void operator()( + DSpanRP electrostatic_potential, + VectorFieldSpan> electric_field, + DViewRP allfdistribu) const + { + auto const grid = ddc::get_domain(allfdistribu); + + FieldRP coords(grid); + ddc::for_each(grid, [&](IndexRP const irp) { coords(irp) = ddc::coordinate(irp); }); + + auto const dom_bsplinesRP = m_rhs_builder.spline_domain(); + + Spline2D allfdistribu_coef(dom_bsplinesRP); + + BSDomainR radial_bsplines(ddc::discrete_space().full_domain().remove_first( + ddc::DiscreteVector {n_overlap_cells})); + BSDomainP polar_domain(ddc::discrete_space().full_domain()); + SplinePolar electrostatic_potential_coef( + PolarBSplinesRP::singular_domain(), + BSDomainRP(radial_bsplines, polar_domain)); + + + m_rhs_builder(allfdistribu_coef, allfdistribu); + PoissonRHSFunction const charge_density_coord(allfdistribu_coef, m_rhs_evaluator); + + m_poisson_solver(charge_density_coord, electrostatic_potential_coef); + + + m_polar_spline_evaluator( + electrostatic_potential.span_view(), + coords.span_cview(), + electrostatic_potential_coef); + + + // Computation of the advection_field: + DFieldRP deriv_r_phi(grid); + DFieldRP deriv_p_phi(grid); + + + m_polar_spline_evaluator.deriv_dim_1( + deriv_r_phi.span_view(), + coords.span_cview(), + electrostatic_potential_coef); + m_polar_spline_evaluator.deriv_dim_2( + deriv_p_phi.span_view(), + coords.span_cview(), + electrostatic_potential_coef); + + ddc::for_each(grid, [&](IndexRP const irp) { + double const r = ddc::coordinate(ddc::select(irp)); + double const th = ddc::coordinate(ddc::select(irp)); + + if (r > m_epsilon) { + CoordRP const coord_rp(r, th); + + Matrix_2x2 J; // Jacobian matrix + m_mapping.jacobian_matrix(coord_rp, J); + Matrix_2x2 inv_G; // Inverse metric tensor + m_mapping.inverse_metric_tensor(coord_rp, inv_G); + + // E = -grad phi + double const electric_field_r + = -deriv_r_phi(irp) * inv_G[0][0] - deriv_p_phi(irp) * inv_G[1][0]; + double const electric_field_th + = -deriv_r_phi(irp) * inv_G[0][1] - deriv_p_phi(irp) * inv_G[1][1]; + + // compute the electric field in the physical domain (Cartesian domain) + ddcHelper::get(electric_field)(irp) + = electric_field_r * J[0][0] + electric_field_th * J[0][1]; + ddcHelper::get(electric_field)(irp) + = electric_field_r * J[1][0] + electric_field_th * J[1][1]; + + } else { + // Linearisation of the electric field + double const th1 = M_PI / 4.; + double const th2 = -M_PI / 4. + 2 * M_PI; + + // --- Value at r = 0: + CoordRP const coord_1_0(0, th1); + CoordRP const coord_2_0(0, th2); + + double const dr_x_1 = m_mapping.jacobian_11(coord_1_0); // dr_x (0, th1) + double const dr_y_1 = m_mapping.jacobian_21(coord_1_0); // dr_y (0, th1) + + double const dr_x_2 = m_mapping.jacobian_11(coord_2_0); // dr_x (0, th2) + double const dr_y_2 = m_mapping.jacobian_21(coord_2_0); // dr_y (0, th2) + + double deriv_r_phi_1 + = m_polar_spline_evaluator + .deriv_dim_1(coord_1_0, electrostatic_potential_coef); + double deriv_r_phi_2 + = m_polar_spline_evaluator + .deriv_dim_1(coord_2_0, electrostatic_potential_coef); + + double const deriv_y_phi_0 = 1 / (dr_y_2 - dr_y_1 * dr_x_2 / dr_x_1) + * (deriv_r_phi_2 - deriv_r_phi_1 * dr_x_2 / dr_x_1); + + double const deriv_x_phi_0 = 1 / dr_x_1 * (deriv_r_phi_1 - deriv_y_phi_0 * dr_y_1); + + // E = -grad phi + double const electric_field_x_0 = -deriv_x_phi_0; + double const electric_field_y_0 = -deriv_y_phi_0; + + + + // --- Value at r = m_epsilon: + CoordRP const coord_rp_epsilon(m_epsilon, th); + Matrix_2x2 J_eps; // Jacobian matrix + m_mapping.jacobian_matrix(coord_rp_epsilon, J_eps); + Matrix_2x2 inv_G_eps; // Inverse metric tensor + m_mapping.inverse_metric_tensor(coord_rp_epsilon, inv_G_eps); + + double deriv_r_phi_epsilon + = m_polar_spline_evaluator + .deriv_dim_1(coord_rp_epsilon, electrostatic_potential_coef); + double deriv_p_phi_epsilon + = m_polar_spline_evaluator + .deriv_dim_2(coord_rp_epsilon, electrostatic_potential_coef); + + + // E = -grad phi + double const electric_field_r = -deriv_r_phi_epsilon * inv_G_eps[0][0] + - deriv_p_phi_epsilon * inv_G_eps[1][0]; + double const electric_field_th = -deriv_r_phi_epsilon * inv_G_eps[0][1] + - deriv_p_phi_epsilon * inv_G_eps[1][1]; + + // compute the electric field in the physical domain (cartesian domain) + double const electric_field_x_epsilon + = electric_field_r * J_eps[0][0] + electric_field_th * J_eps[0][1]; + double const electric_field_y_epsilon + = electric_field_r * J_eps[1][0] + electric_field_th * J_eps[1][1]; + + + // --- Linearisation: + ddcHelper::get(electric_field)(irp) + = electric_field_x_0 * (1 - r / m_epsilon) + + electric_field_x_epsilon * r / m_epsilon; + ddcHelper::get(electric_field)(irp) + = electric_field_y_0 * (1 - r / m_epsilon) + + electric_field_y_epsilon * r / m_epsilon; + } + }); + } +}; diff --git a/tests/geometryRTheta/README.md b/tests/geometryRTheta/README.md index 3b1a29bd0..039e5de26 100644 --- a/tests/geometryRTheta/README.md +++ b/tests/geometryRTheta/README.md @@ -6,6 +6,6 @@ The `tests/geometryRTheta` folder contains all the code describing methods which - [2d_spline_interpolator](./2d_spline_interpolator/README.md) - Tests for the interpolator on 2D polar domain. - - Tests for the polar Poisson solver. +- [polar_poisson](./polar_poisson/README.md) - Tests for the Poisson solver on 2D polar domain. diff --git a/tests/geometryRTheta/polar_poisson/CMakeLists.txt b/tests/geometryRTheta/polar_poisson/CMakeLists.txt index 423dc4efe..83d0c39e6 100644 --- a/tests/geometryRTheta/polar_poisson/CMakeLists.txt +++ b/tests/geometryRTheta/polar_poisson/CMakeLists.txt @@ -29,3 +29,35 @@ foreach(MAPPING_TYPE "CIRCULAR_MAPPING" "CZARNY_MAPPING") set_property(TEST TestPoissonConvergence_${MAPPING_TYPE}_${SOLUTION} PROPERTY TIMEOUT 200) endforeach() endforeach() + + + + +foreach(MAPPING_TYPE "CIRCULAR_MAPPING" "CZARNY_MAPPING") + set(test_name "vlasov_polar_poisson_convergence_${MAPPING_TYPE}") + add_executable("${test_name}" + test_cases.cpp + vlasovpoissonsolver.cpp + ) + target_compile_features("${test_name}" PUBLIC cxx_std_17) + target_link_libraries("${test_name}" + PUBLIC + DDC::DDC + DDC::PDI_Wrapper + paraconf::paraconf + PDI::pdi + sll::splines + gslx::poisson_RTheta + gslx::paraconfpp + Eigen3::Eigen + gslx::geometry_RTheta + ) + target_compile_definitions("${test_name}" PUBLIC -D${MAPPING_TYPE}) + + find_package(Python3 REQUIRED COMPONENTS Interpreter) + + add_test(NAME TestPoissonConvergence_${MAPPING_TYPE}_${SOLUTION} + COMMAND "$" "${CMAKE_CURRENT_SOURCE_DIR}/test_vlasov_poisson.py" + "$") + set_property(TEST TestPoissonConvergence_${MAPPING_TYPE}_${SOLUTION} PROPERTY TIMEOUT 200) +endforeach() diff --git a/tests/geometryRTheta/polar_poisson/README.md b/tests/geometryRTheta/polar_poisson/README.md new file mode 100644 index 000000000..e8ae9d301 --- /dev/null +++ b/tests/geometryRTheta/polar_poisson/README.md @@ -0,0 +1,53 @@ +# Tests on the 2D polar poisson solver + + +The tests implemented in this folder test the 2D polar poisson solver implemented in the `src/geometryRTheta/poisson/` folder +( [poisson](./../../../src/geometryRTheta/poisson/README.md) ). + + +## Polar Poisson solver + +The PolarSplineFEMPoissonSolver is tested on a circular mapping (CircularToCartesian) and on a Czarny mapping (CzarnyToCartesian) with +(the test cases are given in Emily Bourne's thesis [1]). + +The studied equation is + +$$L\phi = - \nabla \cdot (\alpha \nabla \phi) + \beta \phi = \rho, $$ + +with + * $\alpha(r) = \exp\left( - \tanh\left( \frac{r - r_p}{\delta_r} \right) \right)$, with $r_p = 0.7$ and $\delta_r = 0.05$, + * $\beta(r) = \frac{1}{\alpha(r)}$. + + +The defined test functions are: + + * Polar solution: $\phi(x, y) = C r(x,y)^6 (r(x,y) -1)^6 \cos(m\theta)$, + * with $C = 2^{12}1e-4$ and $m = 11$; + * Cartesian solution: $\phi(x,y) = C (1+r(x,y))^6 (1 - r(x,y))^6 \cos(2\pi x) \sin(2\pi y)$, + * with $C = 2^{12}1e-4$. + + +The VlasovPoissonSolver is also tested on a circular mapping (CircularToCartesian) and on a Czarny mapping (CzarnyToCartesian) +for the same Poisson equation with the Cartesian solution. + + + + + + +## References + +[1] Emily Bourne, "Non-Uniform Numerical Schemes for the Modelling of Turbulence in the 5D GYSELA Code". December 2022. + + + +## Contents + + * polarpoissonfemsolver.cpp : it tests the PolarSplineFEMPoissonSolver. It solves the Poisson equation for a selected test case. + * test_cases.hpp : it defines RHS (ManufacturedPoissonTest) and exact solutions (PoissonSolution) of the Poisson equation. + * vlasovpoissonsolver.cpp : it tests the VlasovPoissonSolver. It solves the Poisson equation for a selected test case and computes the electric field in the physical domain. + * poisson.yaml : the parameters of the tests. + * test_poisson.py : it launches twice polarpoissonfemsolver.cpp test and checks the convergence order of the solution. + * test_vlasov_poisson.py : it launches twice vlasovpoissonsolver.cpp test and checks the convergence order or the solution and the electric field. + + diff --git a/tests/geometryRTheta/polar_poisson/poisson.yaml b/tests/geometryRTheta/polar_poisson/poisson.yaml index 3a2e8e399..655bdaac4 100644 --- a/tests/geometryRTheta/polar_poisson/poisson.yaml +++ b/tests/geometryRTheta/polar_poisson/poisson.yaml @@ -1,3 +1,3 @@ Mesh: r_size: 32 - p_size: 32 + p_size: 64 diff --git a/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp b/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp index 220ff55c8..dadcb8ee3 100644 --- a/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp +++ b/tests/geometryRTheta/polar_poisson/polarpoissonfemsolver.cpp @@ -20,7 +20,7 @@ #include "polarpoissonsolver.hpp" #include "test_cases.hpp" -using PoissonSolver = PolarSplineFEMPoissonSolver; +using PoissonSolver = PolarSplineFEMPoissonSolver; #if defined(CIRCULAR_MAPPING) using Mapping = CircularToCartesian; @@ -207,6 +207,5 @@ int main(int argc, char** argv) std::cout << "Max error : " << max_err << std::endl; PC_tree_destroy(&conf_voicexx); - return 0; } diff --git a/tests/geometryRTheta/polar_poisson/test_cases.hpp b/tests/geometryRTheta/polar_poisson/test_cases.hpp index 97bb0559e..99fe91330 100644 --- a/tests/geometryRTheta/polar_poisson/test_cases.hpp +++ b/tests/geometryRTheta/polar_poisson/test_cases.hpp @@ -10,10 +10,20 @@ #include "geometry.hpp" +/** + * @brief Base class for the exact solutions of the Poisson equation. + * + * @see PolarSplineFEMPoissonSolver + * @see VlasovPoissonSolver + */ template class PoissonSolution { public: + /** + * @brief Type the mapping function which converts the logical (polar) + * coordinates into the physical (Cartesian) coordinates. + */ using coordinate_converter_type = CurvilinearToCartesian; private: @@ -23,24 +33,61 @@ class PoissonSolution using RDimP = typename CurvilinearToCartesian::curvilinear_tag_p; protected: + /** + * @brief The mapping function which converts the logical (polar) + * coordinates into the physical (Cartesian) coordinates. + */ CurvilinearToCartesian const& m_coordinate_converter; public: + virtual ~PoissonSolution() = default; + + /** + * @brief Instantiate a PoissonSolution. + * + * @param[in] coordinate_converter + * The mapping function which converts the logical (polar) + * coordinates into the physical (Cartesian) coordinates. + */ PoissonSolution(CurvilinearToCartesian const& coordinate_converter) : m_coordinate_converter(coordinate_converter) { } + /** + * @brief Get the value of the function at a given coordinate. + * + * @param[in] coord + * The given polar coordinate. + * + * @return The value of the function at a given coordinate. + */ virtual double operator()(ddc::Coordinate const& coord) const = 0; }; +/** + * @brief Define a curvilinear solution of the Poisson equation. + * + * The solution is given by + * * @f$ \phi(x, y) = C r(x,y)^6 (r(x,y) -1)^6 \cos(m\theta) @f$, + * + * with @f$C = 2^{12}1e-4 @f$ and @f$ m = 11 @f$. + */ template class CurvilinearSolution : public PoissonSolution { public: + /** + * @brief Instantiate a CurvilinearSolution. + * + * @param[in] coordinate_converter + * The mapping function which converts the logical (polar) + * coordinates into the physical (Cartesian) coordinates. + */ CurvilinearSolution(CurvilinearToCartesian const& coordinate_converter) : PoissonSolution(coordinate_converter) { } + double operator()(ddc::Coordinate const& coord) const final { const double s = ddc::get(coord); @@ -49,14 +96,38 @@ class CurvilinearSolution : public PoissonSolution } }; +/** + * @brief Define a Cartesian solution of the Poisson equation. + * + * The solution is given by + * * @f$ \phi (x,y) = C (1+r(x,y))^6 (1 - r(x,y))^6 \cos(2\pi x) \sin(2\pi y) @f$, + * + * with @f$C = 2^{12}1e-4 @f$. + * + * Its x-derivative is + * * @f$ \partial_x \phi(x,y) = C \cdot 12 (r(x,y)^2 -1)^5 \cdot r(x,y) \cdot \partial_x r(x,y)\cdot \cos(2\pi x) \sin(2\pi y) + * - C (r(x,y)^2 -1)^6 \cdot 2\pi \sin(2\pi x) \sin(2\pi y), @f$ + * + * and its y-derivative, + * * @f$ \partial_y \phi(x,y) = C \cdot 12 (r(x,y)^2 -1)^5 \cdot r(x,y) \cdot \partial_y r(x,y) \cdot \cos(2\pi x) \sin(2\pi y) + * + C (r(x,y)^2 -1)^6 \cdot 2\pi \cos(2\pi x) \cos(2\pi y). @f$ + */ template class CartesianSolution : public PoissonSolution { public: + /** + * @brief Instantiate a CartesianSolution. + * + * @param[in] coordinate_converter + * The mapping function which converts the logical (polar) + * coordinates into the physical (Cartesian) coordinates. + */ CartesianSolution(CurvilinearToCartesian const& coordinate_converter) : PoissonSolution(coordinate_converter) { } + double operator()(ddc::Coordinate const& coord) const final { const double s = ddc::get(coord); @@ -67,24 +138,114 @@ class CartesianSolution : public PoissonSolution return 1e-4 * ipow(1 + s, 6) * ipow(1 - s, 6) * std::cos(2 * M_PI * x) * std::sin(2 * M_PI * y) / ipow(0.5, 12); } + + /** + * @brief Get the value of the x-derivative of the function at a given coordinate. + * + * @param[in] coord + * The given polar coordinate. + * + * @return The value of the x-derivative of the function at a given coordinate. + */ + double derivative_x(ddc::Coordinate const& coord) const + { + const ddc::Coordinate cart_coord + = PoissonSolution::m_coordinate_converter(coord); + const double x = ddc::get(cart_coord); + const double y = ddc::get(cart_coord); + + const double r = ddc::get(coord); + const double dx_r + = PoissonSolution::m_coordinate_converter.inv_jacobian_11( + coord); + + const double C = 1e-4 * ipow(2., 12); + return C * std::sin(2 * M_PI * y) + * (12. * r * dx_r * ipow(r * r - 1, 5) * std::cos(2 * M_PI * x) + - 2 * M_PI * ipow(r * r - 1, 6) * std::sin(2 * M_PI * x)); + } + + /** + * @brief Get the value of the y-derivative of the function at a given coordinate. + * + * @param[in] coord + * The given polar coordinate. + * + * @return The value of the y-derivative of the function at a given coordinate. + */ + double derivative_y(ddc::Coordinate const& coord) const + { + const ddc::Coordinate cart_coord + = PoissonSolution::m_coordinate_converter(coord); + const double x = ddc::get(cart_coord); + const double y = ddc::get(cart_coord); + + const double r = ddc::get(coord); + const double dy_r + = PoissonSolution::m_coordinate_converter.inv_jacobian_12( + coord); + + const double C = 1e-4 * ipow(2., 12); + return C * std::cos(2 * M_PI * x) + * (12. * r * dy_r * ipow(r * r - 1, 5) * std::sin(2 * M_PI * y) + + 2 * M_PI * ipow(r * r - 1, 6) * std::cos(2 * M_PI * y)); + } }; +/** + * @brief Defining the corresponding RHS of the Poisson equation + * for a given exact solution. + */ template class ManufacturedPoissonTest { public: + /** + * @brief Type the choosen solution of the Poisson equation. + */ using CurvilinearToCartesian = typename ChosenSolution::coordinate_converter_type; private: CurvilinearToCartesian const& m_coordinate_converter; public: + /** + * @brief Instantiate a ManufacturedPoissonTest. + * + * @param[in] coordinate_converter + * The mapping function which converts the logical (polar) + * coordinates into the physical (Cartesian) coordinates. + */ ManufacturedPoissonTest(CurvilinearToCartesian const& coordinate_converter) : m_coordinate_converter(coordinate_converter) { } + /** + * @brief Get the value of the RHS at the O-point. + * + * @param[in] coord + * The given polar coordinate. + * + * @return The value of the RHS at the O-point. + */ double solution_at_pole(ddc::Coordinate const& coord) const; + /** + * @brief Get the value of the RHS on the domain except the O-point. + * + * @param[in] coord + * The given polar coordinate. + * + * @return The value of the RHS on the domain except the O-point. + */ double non_singular_solution(ddc::Coordinate const& coord) const; + /** + * @brief Get the value of the RHS at any point of the domain. + * + * @param[in] coord + * The given polar coordinate. + * + * @return The value of the RHS at any point of the domain. + */ double operator()(ddc::Coordinate const& coord) const { if (ddc::get(coord) == 0.0) { diff --git a/tests/geometryRTheta/polar_poisson/test_vlasov_poisson.py b/tests/geometryRTheta/polar_poisson/test_vlasov_poisson.py new file mode 100644 index 000000000..7c2eda0e5 --- /dev/null +++ b/tests/geometryRTheta/polar_poisson/test_vlasov_poisson.py @@ -0,0 +1,70 @@ +#!/usr/bin/python3 +""" +A file used for testing polarpoissonfemsolver.cpp. +It provides 2 input files of different sizes and uses the results +to calculate the order of convergence. +Cubic splines should produce 4th order convergence. +As few points are used the convergence order will not +yet have converged so we simply verify that order > 3.5. + +Inputs: the executable associated with the file polarpoissonfemsolver.cpp. +""" +import subprocess +import sys + +import numpy as np + +executable = sys.argv[1] + +with open("poisson.yaml", "w", encoding="utf-8") as f: + print("Mesh:", file=f) + print(" r_size: 32", file=f) + print(" p_size: 64", file=f) + +with subprocess.Popen([executable, "poisson.yaml"], stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True) as p: + out, err = p.communicate() + assert p.returncode == 0 + +print(out) +if err: + print(err) + assert False + +out_lines = out.split('\n') +error_64 = [float(l.split(' ')[4]) for l in out_lines if "Max error function : " in l][0] +error_dx_64 = [float(l.split(' ')[4]) for l in out_lines if "Max error Ex : " in l][0] +error_dy_64 = [float(l.split(' ')[4]) for l in out_lines if "Max error Ey : " in l][0] + + +with open("poisson.yaml", "w", encoding="utf-8") as f: + print("Mesh:", file=f) + print(" r_size: 64", file=f) + print(" p_size: 128", file=f) + +with subprocess.Popen([executable, "poisson.yaml"], stdout=subprocess.PIPE, stderr=subprocess.PIPE, text=True) as p: + out, err = p.communicate() + assert p.returncode == 0 + +print(out) +if err: + print(err) + assert False + +out_lines = out.split('\n') +error_128 = [float(l.split(' ')[4]) for l in out_lines if "Max error function : " in l][0] +error_dx_128 = [float(l.split(' ')[4]) for l in out_lines if "Max error Ex : " in l][0] +error_dy_128 = [float(l.split(' ')[4]) for l in out_lines if "Max error Ey : " in l][0] + + +order = np.log(error_64/error_128) / np.log(2) +order_dx = np.log(error_dx_64/error_dx_128) / np.log(2) +order_dy = np.log(error_dy_64/error_dy_128) / np.log(2) + +print("Measured order of the electric potential : ", order) +print("Measured order of the electric field in x : ", order_dx) +print("Measured order of the electric field in y : ", order_dy) + + +assert 3.5 < order +assert 2.5 < order_dx +assert 2.5 < order_dy diff --git a/tests/geometryRTheta/polar_poisson/vlasovpoissonsolver.cpp b/tests/geometryRTheta/polar_poisson/vlasovpoissonsolver.cpp new file mode 100644 index 000000000..c69af1eac --- /dev/null +++ b/tests/geometryRTheta/polar_poisson/vlasovpoissonsolver.cpp @@ -0,0 +1,206 @@ +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include + +#include "geometry.hpp" +#include "paraconfpp.hpp" +#include "params.yaml.hpp" +#include "polarpoissonsolver.hpp" +#include "test_cases.hpp" +#include "vlasovpoissonsolver.hpp" + + + +using PoissonSolver = PolarSplineFEMPoissonSolver; + +#if defined(CIRCULAR_MAPPING) +using Mapping = CircularToCartesian; +#elif defined(CZARNY_MAPPING) +using Mapping = CzarnyToCartesian; +#endif +using DiscreteMapping = DiscreteToCartesian; + + + +using LHSFunction = CartesianSolution; +using RHSFunction = ManufacturedPoissonTest; + + +namespace fs = std::filesystem; + +int main(int argc, char** argv) +{ + ::ddc::ScopeGuard scope(argc, argv); + + PC_tree_t conf_voicexx; + if (argc == 2) { + conf_voicexx = PC_parse_path(fs::path(argv[1]).c_str()); + } else if (argc == 3) { + if (argv[1] == std::string_view("--dump-config")) { + std::fstream file(argv[2], std::fstream::out); + file << params_yaml; + return EXIT_SUCCESS; + } + } else { + std::cerr << "usage: " << argv[0] << " [--dump-config] " << std::endl; + return EXIT_FAILURE; + } + PC_errhandler(PC_NULL_HANDLER); + + std::chrono::time_point start_time + = std::chrono::system_clock::now(); + std::chrono::time_point end_time; + + CoordR const r_min(0.0); + CoordR const r_max(1.0); + IVectR const r_size(PCpp_int(conf_voicexx, ".Mesh.r_size")); + + CoordP const p_min(0.0); + CoordP const p_max(2.0 * M_PI); + IVectP const p_size(PCpp_int(conf_voicexx, ".Mesh.p_size")); + + std::vector r_knots(r_size + 1); + std::vector p_knots(p_size + 1); + + double const dr((r_max - r_min) / r_size); + double const dp((p_max - p_min) / p_size); + for (int i(0); i < r_size + 1; ++i) { + r_knots[i] = CoordR(r_min + i * dr); + } + for (int i(0); i < p_size + 1; ++i) { + p_knots[i] = CoordP(p_min + i * dp); + } + + // Creating mesh & supports + ddc::init_discrete_space(r_knots); + ddc::init_discrete_space(p_knots); + + ddc::init_discrete_space(SplineInterpPointsR::get_sampling()); + ddc::init_discrete_space(SplineInterpPointsP::get_sampling()); + + IDomainR const interpolation_domain_R(SplineInterpPointsR::get_domain()); + IDomainP const interpolation_domain_P(SplineInterpPointsP::get_domain()); + IDomainRP const grid(interpolation_domain_R, interpolation_domain_P); + + SplineRBuilder const r_builder(interpolation_domain_R); + SplinePBuilder const p_builder(interpolation_domain_P); + SplineRPBuilder const builder(grid); + +#if defined(CIRCULAR_MAPPING) + const Mapping mapping; +#elif defined(CZARNY_MAPPING) + const Mapping mapping(0.3, 1.4); +#endif + SplineEvaluator2D evaluator( + g_null_boundary_2d, + g_null_boundary_2d, + g_null_boundary_2d, + g_null_boundary_2d); + DiscreteMapping const discrete_mapping + = DiscreteMapping::analytical_to_discrete(mapping, builder, evaluator); + + ddc::init_discrete_space(discrete_mapping, r_builder, p_builder); + + auto const dom_bsplinesRP = builder.spline_domain(); + + DFieldRP coeff_alpha(grid); + DFieldRP coeff_beta(grid); + DFieldRP x(grid); + DFieldRP y(grid); + + ddc::for_each(grid, [&](IndexRP const irp) { + coeff_alpha(irp) + = std::exp(-std::tanh((ddc::coordinate(ddc::select(irp)) - 0.7) / 0.05)); + coeff_beta(irp) = 1.0 / coeff_alpha(irp); + const ddc::Coordinate coord(ddc::coordinate(irp)); + auto const cartesian_coord = mapping(coord); + x(irp) = ddc::get(cartesian_coord); + y(irp) = ddc::get(cartesian_coord); + }); + + Spline2D coeff_alpha_spline(dom_bsplinesRP); + Spline2D coeff_beta_spline(dom_bsplinesRP); + + builder(coeff_alpha_spline, coeff_alpha); + builder(coeff_beta_spline, coeff_beta); + + Spline2D x_spline_representation(dom_bsplinesRP); + Spline2D y_spline_representation(dom_bsplinesRP); + + builder(x_spline_representation, x); + builder(y_spline_representation, y); + + end_time = std::chrono::system_clock::now(); + std::cout << "Setup time : " + << std::chrono::duration_cast(end_time - start_time) + .count() + << "ms" << std::endl; + start_time = std::chrono::system_clock::now(); + + PoissonSolver const poisson_solver(coeff_alpha_spline, coeff_beta_spline, discrete_mapping); + VlasovPoissonSolver const + vlasov_poisson_solver(mapping, builder, evaluator, poisson_solver); + + + end_time = std::chrono::system_clock::now(); + std::cout << "Poisson initialisation time : " + << std::chrono::duration_cast(end_time - start_time) + .count() + << "ms" << std::endl; + + + LHSFunction lhs(mapping); + RHSFunction rhs(mapping); + FieldRP coords(grid); + DFieldRP result(grid); + VectorField> electric_field(grid); + + ddc::for_each(grid, [&](IndexRP const irp) { coords(irp) = CoordRP(ddc::coordinate(irp)); }); + + + DFieldRP rhs_vals(grid); + ddc::for_each(grid, [&](IndexRP const irp) { rhs_vals(irp) = rhs(coords(irp)); }); + + vlasov_poisson_solver(result.span_view(), electric_field, rhs_vals); + + + std::cout << "Solver time : " + << std::chrono::duration_cast(end_time - start_time) + .count() + << "ms" << std::endl; + + + double max_err_x = 0.0; + double max_err_y = 0.0; + ddc::for_each(grid, [&](IndexRP const irp) { + const double diff_x + = fabs(ddcHelper::get(electric_field)(irp) + lhs.derivative_x(coords(irp))); + const double diff_y + = fabs(ddcHelper::get(electric_field)(irp) + lhs.derivative_y(coords(irp))); + max_err_x = max_err_x > diff_x ? max_err_x : diff_x; + max_err_y = max_err_y > diff_y ? max_err_y : diff_y; + }); + std::cout << "Max error Ex : " << max_err_x << std::endl; + std::cout << "Max error Ey : " << max_err_y << std::endl; + + + double max_err = 0.0; + ddc::for_each(grid, [&](IndexRP const irp) { + const double err = fabs(result(irp) - lhs(coords(irp))); + max_err = max_err > err ? max_err : err; + }); + std::cout << "Max error function : " << max_err << std::endl << std::endl; + + PC_tree_destroy(&conf_voicexx); + return 0; +}