Parabolic profile

[cbontoiu]

Hi,

I want to simulate a target which looks like a hollow thick tube, but it needs to have a parabolic profile for the atoms density starting from its core towards the outer edge. The curve I want is shown in this image, where the two vertical blue dashed line indicate the inner (min radius) and outer (max radius) edge of the hollow tube. The expression of the parabolic line is

density = Rbase + Rfact*pow(r, 2);
Rfact = 0.52e39;  Rbase = 0.9e27; 

and the density is in atoms/m^3 because I want to start with unionized atoms

#define _USE_MATH_DEFINES
#pragma once
#include <cmath>
#include "picongpu/particles/densityProfiles/profiles.def"
#include <pmacc/preprocessor/struct.hpp>

namespace picongpu {
    namespace SI {
        constexpr float_64 BASE_DENSITY_SI = 1;
    } // end of SI namespace

namespace densityProfiles {
    constexpr double SIDEY         = 2.000e-05; // [m] longitudinal y-size
    constexpr double SHIFTZ        = 3.500e-06; // [m] horizontal z-middle
    constexpr double SHIFTX        = 3.500e-06; // [m] vertical x-middle
    constexpr double MARGINY_LEFT  = 2.500e-06; // [m] left margin
    constexpr double MARGINY_RIGHT = 2.500e-06; // [m] right margin
    constexpr double INRAD         = 3.000e-07; // [m] inner target radius
    constexpr double OUTRAD        = 2.350e-06; // [m] outer target radius
    constexpr double Rfact         = 5.200e+38; // [m^-3] atoms density scaling factor
    constexpr double Rbase         = 9.000e+26; // [m^-3] atoms density baseline

    struct Functor {
        HDINLINE float_X operator()( const floatD_64& position_SI, const float3_64& cellSize_SI ) {
              const float_64 x( position_SI.x() );
              const float_64 y( position_SI.y() );
              const float_64 z( position_SI.z() );
             float_64 dens = 0.0;
             double r = sqrt(pow(z - SHIFTZ, 2) + pow(x - SHIFTX, 2));
             if( y >= MARGINY_LEFT && y <= SIDEY - MARGINY_RIGHT && r <= OUTRAD && r >= INRAD )
                   dens = Rbase + Rfact*pow(r, 2);
              // safety check: all parts of the function MUST be > 0
              dens *= float_64( dens >= 0.0 );
              return dens;
        } // end density operator
        }; // end of Functor

        // definition of free formula profile
    using Geometry = FreeFormulaImpl< Functor >;
} // end of densityProfiles namespace
} // end of picongpu namespace

Please not that I set BASE_DENSITY_SI = 1 but dens = Rbase + Rfact*pow(r, 2); should account for the correct scaling if, internally PIConGPU take the product between BASE_DENSITY_SI and dens.

It seems that my approach is not fine, probably due to BASE_DENSITY_SI = 1 and here are some strange behaviors I have seen:


Do you see a better solution to my problem, please. Thank you

[psychocoderHPC]

Could you give a short explanation of what your issue is, I see the images but do not know what the result should be. If you miss the gradience place have a look into particles.param and the value MIN_WEIGHTING it could be because of the normalization to the BASE_DENSITY_SI your typical weighting is less than MIN_WEIGTING and therefore only one particle with an increased will be added. Could you please post the beginning of stdout that we can see your UNITS and the typical macro particle weighting?

[cbontoiu]

It seems that regions of the simulations domain are completely removed. Not only in the images above but also in the particle data written. For those regions the electric field is blank as well. My MIN_WEIGHTING = 1.000e-06; but the average weight is 883

WEIGHT = CELL_WIDTH_SI * CELL_HEIGHT_SI * CELL_DEPTH_SI * averageAtomsDensity  / PPC
averageAtomsDensity = AllAtoms / (volAll - volBore);  # atoms/m^3

# SIMULATION DOMAIN -------------------------------- 
uXspace            = 2.034884e-08 # [m] mesh unit cell in x direction
uYspace            = 2.016129e-08 # [m] mesh unit cell in y direction
uZspace            = 2.034884e-08 # [m] mesh unit cell in z direction
CELL_VOLUME        = 8.348290e-24 # [m^3] cell volume
WEIGHT             = 883.62 

I attach the project: input.zip

[psychocoderHPC]

I tested your example and this shows up

   Estimates are based on DensityRatio to BASE_DENSITY of each species
   (see: density.param, speciesDefinition.param).
   It and does not cover other forms of initialization
PIConGPUVerbose PHYSICS(1) | species e: omega_p * dt <= 0.1 ? (omega_p * dt = 2.20373e-15)
PIConGPUVerbose PHYSICS(1) | species C: omega_p * dt <= 0.1 ? (omega_p * dt = 8.90803e-17)
PIConGPUVerbose PHYSICS(1) | y-cells per wavelength: 14.88
PIConGPUVerbose PHYSICS(1) | macro particles per device: 181764096
PIConGPUVerbose PHYSICS(1) | typical macro particle weighting: 2.78276e-24
PIConGPUVerbose PHYSICS(1) | UNIT_SPEED 2.99792e+08
PIConGPUVerbose PHYSICS(1) | UNIT_TIME 3.9063e-17
PIConGPUVerbose PHYSICS(1) | UNIT_LENGTH 1.17108e-08
PIConGPUVerbose PHYSICS(1) | UNIT_MASS 2.53493e-54
PIConGPUVerbose PHYSICS(1) | UNIT_CHARGE 4.45848e-43
PIConGPUVerbose PHYSICS(1) | UNIT_EFIELD 4.36348e+13
PIConGPUVerbose PHYSICS(1) | UNIT_BFIELD 145550
PIConGPUVerbose PHYSICS(1) | UNIT_ENERGY 2.27828e-37
PIConGPUVerbose PHYSICS(1) | Resolving Debye length for species "e"?
PIConGPUVerbose PHYSICS(1) | Estimate used momentum variance in 11900 supercells with at least 10 macroparticles each
PIConGPUVerbose PHYSICS(1) | 11900 (100 %) supercells had local Debye length estimate not resolved by a single cell
PIConGPUVerbose PHYSICS(1) | Estimated weighted average temperature 0 keV and corresponding Debye length 0 m.

Because of the normalization of PIConGPU in unit.param the typical macro weighting is 2.78276e-24. IMO this is the reason for the issue. I tried to increase the BASE_DENSITY by 1e24 and divide the density in density.param by 1e24. This brought back the typical weighting and the units into a reasonable range.

Could you please try

namespace picongpu
{
    namespace SI
    {
        constexpr float_64 BASE_DENSITY_SI = 1e24;
    } // namespace SI

    namespace densityProfiles
    {
        constexpr double SIDEY = 2.000e-05; // [m] longitudinal y-size
        constexpr double SHIFTZ = 3.500e-06; // [m] horizontal z-middle
        constexpr double SHIFTX = 3.500e-06; // [m] vertical x-middle
        constexpr double MARGINY_LEFT = 2.500e-06; // [m] left margin
        constexpr double MARGINY_RIGHT = 2.500e-06; // [m] right margin
        constexpr double INRAD = 3.000e-07; // [m] inner target radius
        constexpr double OUTRAD = 2.350e-06; // [m] outer target radius
        constexpr double Rfact = 5.200e+38; // [m^-3] atoms density scaling factor
        constexpr double Rbase = 9.000e+26; // [m^-3] atoms density baseline

        struct Functor
        {
            HDINLINE float_X operator()(const floatD_64& position_SI, const float3_64& cellSize_SI)
            {
                const float_64 x(position_SI.x());
                const float_64 y(position_SI.y());
                const float_64 z(position_SI.z());
                float_64 dens = 0.0;
                double r = sqrt(pow(z - SHIFTZ, 2) + pow(x - SHIFTX, 2));
                if(y >= MARGINY_LEFT && y <= SIDEY - MARGINY_RIGHT && r <= OUTRAD && r >= INRAD)
                    dens = Rbase + Rfact * pow(r, 2);
                // safety check: all parts of the function MUST be > 0
                dens *= float_64(dens >= 0.0);
                return dens / SI::BASE_DENSITY_SI;
            } // end density operator

        }; // end of Functor

        // definition of free formula profile
        using Geometry = FreeFormulaImpl<Functor>;

    } // namespace densityProfiles

} // namespace picongpu

In general this should not change your setup but avoid that particles get not initialized because of a weighting less than MIN_WEIGHTING

[cbontoiu]

Hi @psychocoderHPC

Yes, indeed this was the solution. I will close the ticket. Thank you for your help!