Picongpu: Density generation

[cbontoiu]

I need to generate ribbons of carbon atoms density and got into a problem as you can see in the file attached. I would ask help with the following issues:

// 1 ...............................
It seems that my new constants declared within namespace densityProfiles {} are accepted by pic-build and I can I use them to generate variables y0(), y1(), ... which enter in the if statements.

However I want to automate the process and use to y(position_SI.y()) function to fill in an array double positions[4*NT] such that density can be assigned in a for loop like
for(int num = 0; num {
if(y >= positions[num] && y <= positions[num+1])
dens = 1.0;
}

However trying to use the array got me the error "picongpu::densityProfiles::positions" cannot be directly written in a device function.
What are my options in this case?

EDIT: THIS WORKS MOVING THE ARRAY DECLARATION INSIDE HDINLINE float_X
operator()( const floatD_64& position_SI, const float3_64& cellSize_SI )
{}

// 2...................
I am curious what is
HDINLINE float_X
operator()( const floatD_64& position_SI, const float3_64& cellSize_SI )
{}
since we declare the y-distribution inside it? Does it mean that the x-distribution cannot be controlled? How would I build regions of discrete density in x- and y- directions independently?

Thank you!
Cristian

density.txt

[sbastrakov]

Hello @cbontoiu .

For the first question, the issue is that global double positions[4*NT]; will stay on the host side, which is not compatible with the density computation, in your case FlatFoilWithRampFunctor::operator(), being on the device side. What you can do is indeed moving the array declaration and definition inside the operator(). Optionally, the array initialization can be separated to a new function, in this case please do not forget to make it HDINLINE. Then the array will be created each time the function is called. This is some computational redundancy but should be affortable generally, and especially for setting things up. Alternatively, you can make the array be a member of the struct and initialize it in a constructor.

Not sure if I got the second question right. That operator() returns normalized density value based on spatial position. The position is of type floatD_64: 2d or 3d vector depending on the simulation dimensionality. So via .x() and .y() of position_SI one can fully control the XY distribution (and full 3d additionally with .z())

[cbontoiu]

Great! Thanks.
Now, given
const float_64 y( position_SI.y() * 1.e9 );
how would I replace
constexpr float_64 y0(margin);
positions[0] = y0;

by something more direct like this
positions[0] = y0(margin);
(which does not work)?

Now the error is "error: cannot determine which instance of overloaded function "y0" is intended".

[sbastrakov]

As far as I understand, these y0, y1, ... y5 variables are merely used for initialization of the positions array. In this case you can simply do

positions[0] = margin;
positions[1] = positions[0] + wall_thck; // here use positions[0] instead of y0, etc.
...

Independent of this change, to hide away this initialization, you can make a function like HDINLINE void initPositions( double positions[] ) { /* init the positions here */ } and in the operator() just have smth like

double positions[4*NT];
initPositions( positions );

[cbontoiu]

OK, it works but in this case what about the
const float_64 y( position_SI.y() * 1.e9 );
and similarly for x and z?
Aren't they essential in generating the position?
Doing this assignment just skips calling these functions.

[sbastrakov]

My previous post was just a technical rewriting of your code attached earlier. Which was expressed in terms of y only in both your original version and my suggested one. Do you want the ribbon boundaries to be defined not just in terms of y, but somehow configured in the 3d space? That is totally possible, perhaps with some more lines of code. We could help with that, just need more information on the configuration then.

[sbastrakov]

Please note that this function we are talking about does not generate any positions at all. It just computes the distribution density values at the given position and returns it. It is a particle generator part of the code that will call this function at various positions and sample macroparticles accordingly. As a user you do not directly control this logic, just provide the density computation functor.

[cbontoiu]

Yes, I ultimately want to achieve a distribution of single carbon atoms such as in the picture attached with

• inner radius: 1.35 nm
• carbon-carbon bond length: 0.1421 nm

but because I don't know how to do it and I am learning to use the code and converting myself from Java to C++, I thought I would start with a 2D approximation.

To the Physics part, these things are necessary for setting a correct target:

• simultaneously implement ionization thresholds: 11.26 eV, 24.38 eV, 47.89 eV, 64.49 eV, 392.1 eV and 490 eV for the six electrons;
• setting the carbon atoms as fixed or greatly limit their mobility, maybe as a function of the local ionization degree, since their bonds are initially very strong;
• allowing ion-electron recombination;
• setting only one electron per atom as a free electron, confined to the surface of the tube or within some shell region along it.

[sbastrakov]

Thanks for clarifications.

With the density one has two options. Either set it up as a distribution function, as you started doing, then you do not have explicit control over where exactly macroparticles are placed, they will be just sampled according to this function. For simple regular patterns, one could also implement it explicitly, e.g. what EveryNthCellImpl does. Then one has a full control over placement and other particle initial conditions, but with more programming effort (since we do not have an existing pattern for this). Sorry, I am still not sure which of the two options you would prefer.

With the physics part perhaps @PrometheusPi or @n01r could comment.

[cbontoiu]

Thanks for the help and availability to look into this problem. I prefer the explicit macroparticle pattern but I didn't know that this is possible. So I started by adapting an example from those coming with the code. Which other example would be closer to the macroparticle pattern approach?

In fact, I think the most physically correct approach would contain (1) the ions (single ionization) distributed in a pattern as shown by the picture, using the macroparticle pattern and (2) the free electrons distributed in a thin shell along the tube shape with uniform density. I am not sure if the two approaches can be reconciliated but I am eager to explore.

[sbastrakov]

[cbontoiu]

[sbastrakov]

[cbontoiu]

[n01r]