ComputationalRadiationPhysics / picongpu

Performance-Portable Particle-in-Cell Simulations for the Exascale Era :sparkles:
https://picongpu.readthedocs.io
Other
694 stars 218 forks source link

total charge #4518

Closed cbontoiu closed 1 year ago

cbontoiu commented 1 year ago

Hi

When I want to work with unionized Carbon atoms my speciesInitialization.param file is:

#pragma once
#include "picongpu/particles/InitFunctors.hpp"
namespace picongpu {
   namespace particles {
      using InitPipeline = bmpl::vector<
         CreateDensity< densityProfiles::Geometry, startPosition::Random2ppc, Carbon >,
    Manipulate< manipulators::unary::RandomPosition, Carbon>,
    ManipulateDerive< manipulators::binary::UnboundElectronsTimesWeighting, Carbon, Electrons >
      >; // end InitPipeline
      } // end of particles namespace
   } // end of picongpu namespace

If I want to work with completely ionized Carbon atoms my speciesInitialization.param file is:

#pragma once
#include "picongpu/particles/InitFunctors.hpp"
namespace picongpu {
   namespace particles {
      using InitPipeline = bmpl::vector<
      CreateDensity< densityProfiles::Geometry, startPosition::Random2ppc, Carbon >,
      Manipulate< manipulators::unary::RandomPosition, Carbon>,
      Manipulate< manipulators::SixTimesIonized, Carbon >,
      ManipulateDerive< manipulators::binary::UnboundElectronsTimesWeighting, Carbon, Electrons >
      >; // end InitPipeline
   } // end of particles namespace
} // end of picongpu namespace

where the particle.param file contains

#pragma once
#include "picongpu/particles/manipulators/manipulators.def"
#include "picongpu/particles/startPosition/functors.def"
#include "picongpu/particles/traits/GetAtomicNumbers.hpp"
#include <pmacc/math/operation.hpp>
namespace picongpu {
   namespace particles {
      namespace manipulators {  
         struct SixTimesIonizedImpl {
        template< typename T_Particle > DINLINE void operator()( T_Particle& particle ) {
           constexpr float_X protonNumber = GetAtomicNumbers< T_Particle >::type::numberOfProtons;
           particle[ boundElectrons_ ] = protonNumber - 6.0_X;
       }
    };// end of SixTimesIonizedImpl
    using SixTimesIonized = generic::Free< SixTimesIonizedImpl >;
    } // end of manipulators namespace
   } // end of particles namespace
} // end of picongpu namespace

I calculate the total "free electron" charge in the target as

totCharge = 0.0;
Ecd, info_Ecd = ts_densityQ.get_field(iteration = ITER[snap], field = 'e_all_chargeDensity', plot=False);
for i in range(CELLSZ):
   for j in range(CELLSY):
      for k in range(CELLSX):
         totCharge += Ecd[i][j][k]*uZspace*uYspace*uXspace

where uZspace, uYspace and uXspace are the mesh sizes in the 3D volume.

I would expect no electron charge in the beginning of the simulation for which I used unionized Carbon atoms, but as shown in the image below, the value is identical to the simulation with completely ionized Carbon atoms. Is this the expected behaviour?

ionization

psychocoderHPC commented 1 year ago

@pordyna @wmles Could you please have a look at this issue?

psychocoderHPC commented 1 year ago

IMO both of your cases are equal if you have not defined a chargeRatio. Bound electrons is by default zero https://github.com/ComputationalRadiationPhysics/picongpu/blob/26cbd45b2c978959301d08d7097d687ecdf9b8a1/include/picongpu/param/speciesAttributes.param#L123 so you ions are fully ionized. ManipulateDerive< manipulators::binary::UnboundElectronsTimesWeighting, Carbon, Electrons > in the first case will only adjust the weighting in case you have defined a charge ratio for carbons. In the second case you call Manipulate< manipulators::SixTimesIonized, Carbon >, but this will do nothing because the attribute boundElectrons is already zero. To start with a setup with ions only that are not already fully ionized you should use a manipulator which set boundElectrons to the number of protons.

// particle.param
// sry for the bad name, I am not a physicist ^^
struct NotIonized {
        template< typename T_Particle > DINLINE void operator()( T_Particle& particle ) {
           constexpr float_X protonNumber = GetAtomicNumbers< T_Particle >::type::numberOfProtons;
           particle[ boundElectrons_ ] = protonNumber ;
       }
    };// end of SixTimesIonizedImpl
    using NotIonized = generic::Free< NotIonized >;

// speciesInitializtion.param
// ...
Manipulate< manipulators::NotIonized, Carbon >
// ..

I hope this solves the question. I am wondering why you call ManipulateDerive< manipulators::binary::UnboundElectronsTimesWeighting, Carbon, Electrons > in your first setup for unionized ions. If you use the functor I suggested this will create electrons with zero weighting but maybe I misunderstand something.

PrometheusPi commented 1 year ago

Exactly, as @psychocoderHPC pointed out, both cases create free electrons and you need to use NotIonized. In PIConGPU, the default is fully ionized.

Why are using Manipulate< manipulators::unary::RandomPosition, Carbon>, is the second step of your init pipeline?

psychocoderHPC commented 1 year ago

@cbontoiu Manipulate< manipulators::unary::RandomPosition, Carbon>, should not be required because you use in CreateDensity already a random in cell position. If you skip this manipulator you should have slide faster initialization.

BrianMarre commented 1 year ago

@PrometheusPi @psychocoderHPC @cbontoiu Due to future features the NotIonized Functor has changed, we now have a dedicated functor for charge state initialization that will become important in the future in certain situations.

https://github.com/ComputationalRadiationPhysics/picongpu/blob/26cbd45b2c978959301d08d7097d687ecdf9b8a1/include/picongpu/particles/atomicPhysics/SetToAtomicGroundStateForChargeState.hpp

It is already in use in the atomicPhysics example [picongpu/share/picongpu/examples/atomicPhysics], including in it in the other examples is on my to-do list but I have not found time yet.

the changed SetIonization Functor is exemplary defined here https://github.com/ComputationalRadiationPhysics/picongpu/blob/26cbd45b2c978959301d08d7097d687ecdf9b8a1/share/picongpu/examples/atomicPhysics/include/picongpu/param/particle.param

and used here in the init pipeline https://github.com/ComputationalRadiationPhysics/picongpu/blob/26cbd45b2c978959301d08d7097d687ecdf9b8a1/share/picongpu/examples/atomicPhysics/include/picongpu/param/speciesInitialization.param#L49

cbontoiu commented 1 year ago

Hi, thanks for your help. Here are some questions for my Carbon target:

  1. If I am not interested about Carbon ions motion can I remove DefaultParticleAttributes from using CarbonAttributes = MakeSeq_t<DefaultParticleAttributes, atomicConfigNumber<configNumber_Carbon>, boundElectrons>; if there is any gain in performance of memory?

  2. Where would I define levelVector_Carbon, // name used in species Definition which appears in atomicPhysics.param? In your case this is levelVector_Argon.

2023-03-28_16-17

  1. I couldn't find documentation for how to define the stateDataFileName = "../CarbonLevels.txt"and transitionDataFileName = "../CarbonTransitions.txt"; files. Where are these files in your case?

2023-03-28_16-19

  1. If I want to work with neutral Carbon atoms, how would I change the struct SetIonTo....Ionized struct? Something like this?

      {
         template<typename T_Particle>
         DINLINE void operator()(T_Particle& particle)
            {
                constexpr float_X targetChargeState = 0;
                constexpr float_X numberBoundElectrons
                       = picongpu::traits::GetAtomicNumbers<T_Particle>::type::numberOfProtons;
    
                PMACC_ASSERT_MSG(numberBoundElectrons >= targetChargeState, "too few protons...");
                picongpu::particles::atomicPhysics::SetToAtomicGroundStateForChargeState{}(
                                particle, numberBoundElectrons - targetChargeState);
            } // end DINLINE operator
       }; // end SetIonTo0Ionized

Finally, these are my relevant files for working with completely ionized Carbon ions. Please, could you check to make sure I am not missing anything, and I am aligned with the latest dev version of PIConGPU, as well as being at an optimum for memory and performance? The compilation error is now /home/quasar/LASER3DEFFDENSCNT/input/include/picongpu/param/speciesInitialization.param(12): error: expected a type specifier /usr/include/boost/mpl/transform.hpp(42): error: class "mpl_::na" has no member "state" probably related to state and transition missing files.

speciesInitialization.txt speciesDefinition.txt particle.txt atomicPhysics.txt

Thank you.

BrianMarre commented 1 year ago

regarding 1.) Short answer: No

For ionization to work a particle species needs the following attributes:

and the following flags

(if I remember correctly, the compiler will tell you if I forgot something) You may therefore only reduce the execution time a little by using the Free pusher initialize the ions as non-moving and distribute the electron absolutely uniformly on the GPUs via a custom GPU to cell mapping.

In general yes non moving particles are possible in PIConGPU, we call them probe particles. You need to keep at least the position macro particle attribute, otherwise the position of a particle is not defined. You may remove the momentum attribute to save memory and some compute time, this will create non-moving(probe) macro particles. If you do that you also have to either remove the pusher flag or use the probe pusher(used for probing local field conditions instead of moving the particle), remove the current flag(current deposition relies on the momentum attribute) and may remove the chargeRatio(uncharged particle, probably not what you intended)/massRatio flag (removing flags saves a lot less in memory than removing attributes since they are species constants, but might make a difference if you need to squeeze out the last bit). Unfortunately this is not possible if you want to keep ionization since we need a macro-ion momentum for momentum conservation in the ionization, not currently configurable but can be deactivated by hand, ask me if you really want to do that.

regarding 3.) and 4.) ignore it and remove it from your setup. These options are for only for FLYonPIC, atomicPhysics excited states modelling, which requires considerable additional compute and memory resources, the current dev has a prototype version, new version upcoming upcoming soon.

regarding 5.) looks good :+1:

regarding initializing an ion as charge state 0: Would not recommend, since PIC assumes that all macro-ions are free, and in reality neutral carbon atoms are bound. Therefore we usually start with at least charge State 1, to make sure that

may happen by accident if you did not co-locate them in the init-pipeline, or have a mismatched initial ionization state of the ion species and densityRatio of the electron species,

These shadow background fields are constant and will influence the particle movement.

If you really want to start in a almost charge neutral setup, I would recommend to at least create some once ionized ions and the corresponding co-located free electrons, in addition to the bulk non-ionized electrons, with density ratio of once ionized ions to non-ionized ions according to saha-equation for the initial temperature you want to model

regarding your setup:

BrianMarre commented 1 year ago

a few further points,

And you do no need the atomicPhysics.param, since you do not use FLYonPIC.

EDIT: I was wrong, we actually have a fallback, if a species lacks the densityRatio flag and is put into a DensityWeighting functor a default densityratio of 1 will be used.

cbontoiu commented 1 year ago

@BrianMarre

I am grateful for this extensive reply! It is much appreciated.

I eliminated the atomicPhysics.param and I will stick to ionization state 1 and I want to have the ionizeable carbon ions in the system.

  1. Are the following settings fine for speciesDefinition.param?
/*--------------------------- electrons --------------------------------------*/
/* ratio relative to BASE_CHARGE and BASE_MASS */
value_identifier( float_X, MassRatioElectrons, 1.0 );
value_identifier( float_X, ChargeRatioElectrons, 1.0 );
/*------------------------------- C ------------------------------------------*/
/* ratio relative to BASE_CHARGE and BASE_MASS */
value_identifier( float_X, MassRatioCarbon, 22032.0 );
value_identifier( float_X, ChargeRatioCarbon, -6.0 );
/* ratio relative to BASE_DENSITY (n_e) */
value_identifier( float_X, DensityRatioCarbon, 1. / 1. );
  1. For adding some initial temperature, are the following lines:
    struct TemperatureParam {
    // Initial temperature unit: keV
    static constexpr float_64 temperature = 1.;
    };
    using AddTemperature = unary::Temperature<TemperatureParam>;

    sufficient in namespace manipulators of particle.param? I added one more line in the InitPipeline vector of speciesInitialization.param:

Manipulate<manipulators::AddTemperature, Electrons>

cbontoiu commented 1 year ago

@psychocoderHPC, @PrometheusPi, @BrianMarre

Thank you for your help! I will not mess up again with the ionization state. I'll set it to 1 and then let the laser ionize as much as it can.

Side note: I am working with Carbon in solid-state materials such as nanotubes or graphene. Each atom shares 3 of the 6 electrons in sigma and pi bonds with neighbouring Carbon atoms. I assume that the first ionization potential of such structures is low enough to be equivalent to considering all atoms ionized once, but the choice is arbitrary. When showing LWFA results starting with ionization state 3, one of our reviewers had doubts that the laser pulse (though very intense) can ionize the Carbon atoms in reality, so this is why I started working with neutral atoms in the past.

My model compiles now and this is how the files look like. I have concerns about the boundElectrons variable in using CarbonAttributes = MakeSeq_t<DefaultParticleAttributes, boundElectrons>; where does it come from and do I need to specify it in speciesInitialization.param?

1 4 3

BrianMarre commented 1 year ago

~~your charge sign convention is still wrong, the PIConGPU internal unit of charge is by default the elementary charge. Therefore ChargeRatioElectron = -1.0 and ChargeRatioCarbon = 6.0, otherwise seems fine to me.~~ (EDIT: I was wrong, mea culpa, ChargeRatios values are relative to BASE_CHARGE=-e not UNIT_CHARGE=e)

regarding boundElectron, it is defined in speciesAttributes.param that all works out of the box and you do not need to do anything.

Short side note:

psychocoderHPC commented 1 year ago

regarding boundElectron, it is defined in speciesAttributes.param that all works out of the box and you do not need to do anything.

Do clarify this a little bit more. Each attribute is a value_identifier and the last parameter is the default value if it is not set somehow e.g. via a manipulator in speciesInitialization.param. All pre-defined attributes will automatically be visible in your IO (openPMD) too. using CarbonAttributes = MakeSeq_t<DefaultParticleAttributes, boundElectrons>; is simply extending the list DefaultParticleAttributes by one attribute boundElectrons

BrianMarre commented 1 year ago

@cbontoiu I was wrong, your original ChargeRatio flag values were already correct, sorry. I got Base_Charge=-e and unitCharge=e mixed up.