add the possibility to have a parabolic plasma profile

doc/source/run/parameters.rst

@@ -92,6 +92,10 @@ Plasma parameters
 * ``plasma.radius`` (`float`) optional (default `infinity`)
     Radius of the plasma. Set a value to run simulations in a plasma column.
 
+* ``plasma.parabolic_curvature`` (`float`) optional (default `0.`)
+    Curvature of a parabolic plasma profile. The plasma density is set to
+    :math:`\mathrm{plasma.parabolic\_curvature}*r^2 + \mathrm{plasma.density}`.
+
 * ``plasma.max_qsa_weighting_factor`` (`float`) optional (default `35.`)
     The maximum allowed weighting factor :math:`\gamma /(\psi+1)` before particles are considered
     as violating the quasi-static approximation and are removed from the simulation.

src/particles/PlasmaParticleContainer.H

@@ -45,19 +45,24 @@ public:
      * \param[in] a_thermal_momentum_mean average momentum (3d)
      * \param[in] a_density plasma density (SI)
      * \param[in] a_radius plasma radius. Only particles with x**2+y**2<a_radius**2 are injected
+     * \param[in] a_parabolic_curvature defines the curvature of a transerse plasma profile with
+     *            a_parabolic_curvature*r^2 + a_density
      */
     void InitParticles (
         const amrex::IntVect& a_num_particles_per_cell,
         const amrex::RealVect& a_thermal_momentum_std,
         const amrex::RealVect& a_thermal_momentum_mean,
         const amrex::Real a_density,
-        const amrex::Real a_radius);
+        const amrex::Real a_radius,
+        const amrex::Real a_parabolic_curvature);
 
     amrex::Real m_density {0}; /**< Density of the plasma */
     /** maximum weighting factor gamma/(Psi +1) before particle is regarded as violating
      *  the quasi-static approximation and is removed */
     amrex::Real m_max_qsa_weighting_factor {35.};
     amrex::Real m_radius {std::numeric_limits<amrex::Real>::infinity()}; /**< radius of the plasma */
+    /** defines the parabolic curvature of a plasma profile m_parabolic_curvature*r^2 + m_density */
+    amrex::Real m_parabolic_curvature {0.};
     amrex::IntVect m_ppc {0,0,1}; /**< Number of particles per cell in each direction */
     amrex::RealVect m_u_mean {0,0,0}; /**< Avg momentum in each direction normalized by m*c */
     amrex::RealVect m_u_std {0,0,0}; /**< Thermal momentum in each direction normalized by m*c */

src/particles/PlasmaParticleContainer.cpp

@@ -8,6 +8,7 @@ PlasmaParticleContainer::PlasmaParticleContainer (amrex::AmrCore* amr_core)
     amrex::ParmParse pp("plasma");
     pp.query("density", m_density);
     pp.query("radius", m_radius);
+    pp.query("parabolic_curvature", m_parabolic_curvature);
     pp.query("max_qsa_weighting_factor", m_max_qsa_weighting_factor);
     amrex::Vector<amrex::Real> tmp_vector;
     if (pp.queryarr("ppc", tmp_vector)){
@@ -34,7 +35,7 @@ PlasmaParticleContainer::InitData ()
     reserveData();
     resizeData();
 
-    InitParticles(m_ppc,m_u_std, m_u_mean, m_density, m_radius);
+    InitParticles(m_ppc,m_u_std, m_u_mean, m_density, m_radius, m_parabolic_curvature);
 
     m_num_exchange = TotalNumberOfParticles();
 }

src/particles/PlasmaParticleContainerInit.cpp

@@ -10,7 +10,8 @@ InitParticles (const amrex::IntVect& a_num_particles_per_cell,
                const amrex::RealVect& a_u_std,
                const amrex::RealVect& a_u_mean,
                const amrex::Real a_density,
-               const amrex::Real a_radius)
+               const amrex::Real a_radius,
+               const amrex::Real a_parabolic_curvature)
 {
     HIPACE_PROFILE("PlasmaParticleContainer::InitParticles");
 
@@ -120,9 +121,11 @@ InitParticles (const amrex::IntVect& a_num_particles_per_cell,
                     y >= a_bounds.hi(1) || y < a_bounds.lo(1) ||
                     x*x + y*y > a_radius*a_radius ) continue;
 
+                const amrex::Real radius = std::sqrt(x*x + y*y);
+                const amrex::Real parabolic_factor = a_parabolic_curvature * radius * radius;
+
                 amrex::Real u[3] = {0.,0.,0.};
-                ParticleUtil::get_gaussian_random_momentum(u, a_u_mean,
-                                                           a_u_std);
+                ParticleUtil::get_gaussian_random_momentum(u, a_u_mean, a_u_std);
 
                 ParticleType& p = pstruct[pidx];
                 p.id()   = pid + pidx;
@@ -131,8 +134,8 @@ InitParticles (const amrex::IntVect& a_num_particles_per_cell,
                 p.pos(1) = y;
                 p.pos(2) = z;
 
-                arrdata[PlasmaIdx::w        ][pidx] = a_density * scale_fac;
-                arrdata[PlasmaIdx::w0       ][pidx] = a_density * scale_fac;;
+                arrdata[PlasmaIdx::w        ][pidx] = (a_density + parabolic_factor) * scale_fac;
+                arrdata[PlasmaIdx::w0       ][pidx] = (a_density + parabolic_factor) * scale_fac;
                 arrdata[PlasmaIdx::ux       ][pidx] = u[0] * phys_const.c;
                 arrdata[PlasmaIdx::uy       ][pidx] = u[1] * phys_const.c;
                 arrdata[PlasmaIdx::psi      ][pidx] = 0.;