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<positions.size(); num += 2) { 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]

For a rigid pattern example, please take a look at how we do the pattern for a macroparticle per each N cells. The idea is that depending on the cell position we decide whether it should have a macroparticle in it (so here is the pattern applied which would be different in your case). Then we just return a matching density value to produce 0 or 1 (or can be more) macroparticles in this cell. This provides the right number of macroparticles following the pattern. Then we need to fix the position inside the cell to also follow the spatial pattern like here. (Normally this position would be random to provide a feel of continuous distribution). Btw your attached density.txt already uses it for probes (as our setups normally do).

From your attached picture I am not sure how the 2d pattern looks like for your setup. But, again, once you are able to express it via coordinates and have a fixed position inside a cell (perhaps choose the spatial steps accordingly), should be possible to implement.

Generally, please also note that in case the generic 2d / 3d density function somewhat works, it might be reasonable to just use that one (since it should be much quicker to do) and meanwhile construct the rest of the setup. Changing the initial distribution can be added at any point, from the code point of view it's just a change to a couple of template parameters, independent on the rest.

[cbontoiu]

One more thing here supposing I build a thin cylindrical shell of electrons. In the case of carbon nano-tubes it is known that the axial electrical/thermal conductivity is very high (j = 1e7-1e13 A/cm^2 and k = 3500 W/(mK)) but the same parameters are very small azimuthally. Is there a way to force the electrons to move only axially?

[sbastrakov]

Sorry, I am getting confused by the axial / azimutal coordinates. Could you please repharase using PIConGPU coordinates? Are you by any chance wondering about how particles move in a 2d simulation? In the latter case, the particle coordinates are x and y, each particle is an infinite string in z, but the movement is full 3d, so are the fields and currents.

[cbontoiu]

The simulation must be in 3D overall but the "free" carbon electrons should only move along the longitudinal direction say z and I take this axis as the tubes longitudinal axes. Any motion of the electrons along the picongpu y- and x-axes, which are transverse to the tubes longitudinal axes should be prevented, to account for the physical peculiar facts of carbon nano-tubes. Of course, in doing this the fields are still required to be rendered in 3D.

[n01r]

Hi @cbontoiu, sorry for not answering earlier since I'm occupied with thesis work atm. :)

I would like to address some of your questions:

1. Regarding your initialization: since the pattern is relatively complex, it should still be possible to describe it in code, you might find it helpful to have a look at https://github.com/hightower8083/PoGit. It is an input templater that uses Python to help you set up PIConGPU param and cfg files. @hightower8083 do you think you have some capacities to possibly help with questions if PoGit is used?
2. To restrict the allowed motion for electrons, a custom particle pusher has to be written. Please see the following details for directions:
   • Particle pushers live here. The most commonly used one is the so-called Boris pusher. You would have to create a copy and remove the motion in the unwanted directions.
   • Then you need to add your pusher to the pusher.param which you put inside your cloned and edited example.
   • Be aware that the standard pusher that will be sued for all particles by default is configured in species.param
   • Inside your speciesDefinition.param you finally have to configure your pusher for the electrons like so:

     ...
     particlePusher< particles::pusher::CNTDirectedElectrons >,
     ...

   • What is the central laser wavelength you are aiming to use? I assume you want to operate at intensities where the CNT still stays intact which puts the setup into a regime where a lot of physics happens that is inherently not described by the PIC method. Ionization, for instance, is added as a module on top of the PIC cycle. The models we have implemented so far work best with high power optical or near-infrared lasers to which the quasi-static approximation can be applied.
   • Unfortunately, electron-ion recombination has not been implemented, yet.
   • Since the Carbon ions are very heavy, compared to the electrons, they will move much less, anyway. If you would like to further restrict their motion you might have to write a custom pusher for them as well. Any motion of particles stems from an electromagnetic field interaction and the force is inherently connected to the ion's charge state.

Hope that helps already and you have a better feeling about how to do realize your setup!