Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Thomas-Fermi impact ionization model #1754

Merged

Conversation

n01r
Copy link
Member

@n01r n01r commented Jan 9, 2017

This pull request adds the Thomas-Fermi impact/collisional ionization model to PIConGPU.

The model is implemented using the field ionization framework introducing grid-mapped and particle-position-interpolated density and "temperature" fields as input parameters. For more information see this issue.

The implementation also makes use of multiple FieldTmp slots (#1703).

Preliminary compile- and runtime tests have been successful and show good agreement of HDF5 output and analytical expectation of charge states.

I propose to do a follow-up PR with a new ThinFoil example to free up the LaserWakefield example.

ToDo
  • Clean up the commit history and regroup into meaningful units

Please refrain from commenting on small inconsistencies until the above box has been checked.

  • Replace round2Nearest(avgChargeState) by RNG generated distribution between floor and ceil (Update 2017-01-10)

@n01r
Copy link
Member Author

n01r commented Jan 9, 2017

@ax3l Does it make sense to add this PR to the Collisional Ionization Methods project?

@ax3l ax3l self-requested a review January 9, 2017 15:51
@ax3l ax3l added component: core in PIConGPU (core application) feature labels Jan 9, 2017
@ax3l ax3l self-assigned this Jan 9, 2017
@ax3l ax3l added this to the Next Stable: 0.3.0 / 1.0.0 milestone Jan 9, 2017
@ax3l
Copy link
Member

ax3l commented Jan 9, 2017

yes sure, done!

I propose to do a follow-up PR with a new ThinFoil example to free up the LaserWakefield example.

yes, sounds good!

@n01r n01r force-pushed the topic-ThomasFermiNewFieldTmp branch 3 times, most recently from ef5e08a to ed83deb Compare January 10, 2017 10:58
/** @TODO replace the static_cast<float_64> by casts to float_X
* or leave out entirely and compute everything in PIConGPU scaled units
*/
float_64 const densityUnit = static_cast<float_64>(particleToGrid::derivedAttributes::Density().getUnit()[0]);
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do calculations in numerically stable scaled units if possible.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems not really possible in general right now. Maybe at a later point.

float_64 x = TFAlpha * math::pow(Q,TFBeta_temp);

/* Thomas-Fermi average ionization state */
float_X ZStar = static_cast<float_X>(protonNumber * x / (float_64(1.) + x + math::sqrt(float_64(1.) + float_64(2.)*x)));
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thomas-fermi_to_not_infringe_copyright

Let's see if we can go from float_64 to float_X.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ax3l could you suggest a way to calculate this in more numerically stable units?
I don't really see it right now. Normally this is fine because you need to arrive at the correct units in the end and can just multiply after all calculations back into SI. But since this is a fitted formula what do I do here?

E.g. with a terms like math::pow(R,B) and math::pow(R,C)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yep, it's nasty but possible.
first start by checking the dimensions of the constants/fit parameters (and add them). We can then normalize those quantities in formula accordingly. If the formula is SI already, then it's just simple scaling of the constants before doing the calc.

Copy link
Member Author

@n01r n01r Jan 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

after offline discussion: So do we agree on keeping the calculation in double precision for now since there is no easy way of retracing how and with which units this was fitted?

/** determine the new number of bound electrons from the TF ionization state by rounding to nearest charge state
* @TODO introduce partial macroparticle ionization / ionization distribution at some point
*/
float_X newBoundElectrons = protonNumber - float_X(math::float2int_rn(ZStar));
Copy link
Member Author

@n01r n01r Jan 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Follow this up with a commit that calculates

/* works only for positive charge states */
float_X chargeState  = math::floor(ZStar) + static_cast<float_X>(randNr < math::modf(ZStar))
float_X newBoundElectrons = protonNumber - chargeState;

where randNr is a random number.
If for example ZStar = 3.6 this would yield a charge state of 4+ in 60% and a charge state of 3+ in 40% of the cases.

/** convert kinetic energy in J to "temperature" in eV by assuming an ideal electron gas
* E_kin = 3/2 k*T
*/
float_64 temperature = kinEnergy * UNITCONV_Joule_to_keV * float_64(1.e3) * float_64(2./3.);
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What do we actually want to call an "electron temperature"?

@ax3l , @HighIander : ideas?

ToDo: Implement cut-off at 50 keV (?) / Make the cut-off a user-definable parameter.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, the cut-off is an important problem and there is a very recent scientific discussion about exactly this topic here: Phys. Rev. Lett. 113, 255001 (2014) and in the corresponding comment and reply.
Summarizing, I think a 30 keV cutoff as default is OK for now, but one should allow for a user-defined value and (later) even better would allow for a real fit with a Maxwell to the low-energy part or getting the thermal velocity by any other means "on-the-fly".

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, we should keep the assumption that one can derive a "temperature" in the simple model from the mean with user-defined cut-off and should then rather check against one of the many extensions of TF, if one is better suited non-thermal distributions for our use cases (see the references from last week in David Salzmann, 1998, Ch. 2.4).

@n01r n01r force-pushed the topic-ThomasFermiNewFieldTmp branch from ed83deb to cb6c60e Compare January 24, 2017 10:14
@n01r
Copy link
Member Author

n01r commented Jan 24, 2017

Go ahead @ax3l , @psychocoderHPC you can review from here. 👍 .
No need to introduce new feature parts all in this PR.

@ax3l ax3l changed the title [WIP] Add Thomas-Fermi impact ionization model Add Thomas-Fermi impact ionization model Jan 24, 2017
@ax3l ax3l changed the title Add Thomas-Fermi impact ionization model Thomas-Fermi impact ionization model Jan 24, 2017
@ax3l
Copy link
Member

ax3l commented Jan 27, 2017

@n01r compile error, e.g.

src/picongpu/include/particles/ionization/byCollision/ThomasFermi/AlgorithmThomasFermi.hpp(136):
  error: no instance of function template "PMacc::algorithms::math::fmod" matches the argument list
    argument types are: (picongpu::float_X)

@psychocoderHPC
Copy link
Member

psychocoderHPC commented Jan 27, 2017 via email

@n01r
Copy link
Member Author

n01r commented Feb 3, 2017

Yeah, corrected that offline already and had it in a git stash which I then forgot 🤦‍♂️

@n01r n01r force-pushed the topic-ThomasFermiNewFieldTmp branch from cb6c60e to c2d52bb Compare February 3, 2017 13:01
@n01r
Copy link
Member Author

n01r commented Feb 6, 2017

@ax3l checks are passing, could you review this week, please?

/** \file AlgorithmThomasFermi.hpp
*
* IONIZATION ALGORITHM for the Thomas-Fermi model
*
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you add the original paper and maybe other materials you found relevant for your implementation at this point?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

-> done so in the .def file, very good!

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please remove the empty line and keep one block. had problems with adding newlines in the \file description in the last days.

Write like this please:

/** \file AlgorithmThomasFermi.hpp
 *
 * IONIZATION ALGORITHM for the Thomas-Fermi model
 * - implements the calculation of average ionization degree and changes charge states
 *   by decreasing the number of bound electrons
 * - is called with the IONIZATION MODEL, specifically by setting the flag in @see speciesDefinition.param
 */

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

Copy link
Member

@ax3l ax3l left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Implementation and test results look good to me, it would be nice if you could upload a test result and describe what you did offline to the PR, too.

Sadly it's a bit hard to find out the correct scaling without crawling into details of the fit that was used and badly documented in the original paper. But this will be ok for now.

A few inline comments on wording, modern style for new code and a few questions. Please check your include orders, I already changed a few in #1800 so you can take a look (described there) and adjust here in the new files accordingly.

Thanks for the clean implementation and testing! ✨

* \brief calculation for the Thomas-Fermi pressure ionization model
*
* This model uses local density and "temperature" values as input
* parameters. To be able to speak of a "temperature" an equilibrium state
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Phrasing suggestion: To be able to speak of -> A physical temperature requires a defined equilibrium state ...

* would be required. Typical high power laser-plasma interaction is highly
* non-equilibrated, though. The name "temperature" is kept to illustrate
* the origination from the Thomas-Fermi model. It is nevertheless
* more accurate to think of it as an averaged kinetic energy.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

... to think of it as an averaged kinetic energy, which is not backed by the Model and should therefore be only used with a certain suspicion in such Non-LTE scenarios.

{

/* @TODO replace the float_64 with float_X and make sure the values are scaled to PIConGPU units */
const float_64 protonNumber = GetAtomicNumbers<ParticleType>::type::numberOfProtons;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think both lines can be constexpr

if (chargeState < protonNumber)
{
/* atomic mass number (usually A) A = N + Z */
float_64 massNumber = neutronNumber + protonNumber;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should also be a constexpr

float_64 const densityUnit = static_cast<float_64>(particleToGrid::derivedAttributes::Density().getUnit()[0]);
float_64 const kinEnergyDensityUnit = static_cast<float_64>(particleToGrid::derivedAttributes::EnergyDensity().getUnit()[0]);
/* convert from kinetic energy density to average kinetic energy per particle */
float_64 kinEnergy = (kinEnergyDensity / density) * (kinEnergyDensityUnit / densityUnit);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

const float_64 kinEnergy


#pragma once

#include <boost/type_traits/integral_constant.hpp>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please move this to the bottom, just before stdlib includes

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

there are no stdlib includes here, right? ^^; so this will be placed at the last place in the include list

/* alias for the single macro-particle */
auto particle = ionFrame[localIdx];
/* particle position, used for field-to-particle interpolation */
floatD_X pos = particle[position_];
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

const

typedef typename RNGFactory::GetRandomType<Distribution>::type RandomGen;
RandomGen randomGen;

typedef MappingDesc::SuperCellSize TVec;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

that is a very nasty renaming, because you use it below as actual mapping to get a 3D index inside the supercell.
not sure if the rename is needed at all, otherwise just make it shorter, e.g., using SuperCellSize = MappingDesc::SuperCellSize;

struct ThomasFermi_Impl
{

typedef T_DestSpecies DestSpecies;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if you don't mind, please use the C++11 using NewType = OtherType; since it's more readable for new code.

namespace ionization
{

struct AlgorithmNone;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

do you really need that?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is the fallback like in fieldIonizationCalc.def. There should be other ways to do a fallback, though.

Copy link
Member

@ax3l ax3l Feb 13, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yep, both should be removed. if no algo is used, just don't add the flag to the Particle class in the species. can be done in a follow-up PR for the field ioniz.

@n01r n01r force-pushed the topic-ThomasFermiNewFieldTmp branch from 7933c3b to 7b31142 Compare February 13, 2017 14:09
@n01r
Copy link
Member Author

n01r commented Feb 13, 2017

Merge conflicts in ionizerConfig.param - also still recent fix to ionizationMethods.hpp in #1824.
Need to rebase this against dev before finally uploading some test results and the eventual merge.

@n01r n01r force-pushed the topic-ThomasFermiNewFieldTmp branch from 7b31142 to 35099ff Compare February 13, 2017 14:37
@n01r n01r force-pushed the topic-ThomasFermiNewFieldTmp branch 2 times, most recently from 1e4043b to 89691fa Compare February 14, 2017 13:32
@n01r
Copy link
Member Author

n01r commented Feb 14, 2017

I force-pushed the rebased and amended commit history. The branch was compile- and runtested for a couple of simulation steps. Before I will test the more complicated setup of a spatial temperature gradient I created a little steady-state scenario.

It is in an unphysical regime but that doesn't matter for the model functionality.

We have:

  • neutral nitrogen gas at a homogenous density of 1e25 m^-3
  • electrons at the same density with an initial temperature of 17eV
  • no pushers
  • no current
  • ThomasFermi ionization model (of course)
  • 3 time steps (already 1 is enough)

We expect an average charge state of roughly 4.5.
nitrogen_z_depending_on_t_thomasfermi_theoretical

That leaves us with an expected average number of boundElectrons equal to 2.5.
This is the statistics overview of the boundElectrons attribute of the 4,194,291 nitrogen ions after the first time step. ⬇️
steady_state_nitrogen_ionization_boundelectrons

Looks successful to me.

Copy link
Member

@ax3l ax3l left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Slight changes and rebase pls :)

@@ -121,6 +122,10 @@ struct Hydrogen
* -> circularly polarized lasers
* - Keldysh : Keldysh ionization model
*
* - ThomasFermi : statistical impact ionization based on Thomas-Fermi
* atomic model
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

after rebase, please add this section to ionizerConfig.param now

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

float_64 constexpr nAvogadro = SI::N_AVOGADRO;
float_64 constexpr convM3ToCM3 = 1.e6;

float_64 const convToMassDensity = densityUnit * massNumber / nAvogadro / convM3ToCM3;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

massNumber / nAvogadro / convM3ToCM3 can also be expressed as constexpr :-)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

constexpr float_64 cleverNameYouChoose = massNumber / nAvogadro / convM3ToCM3;
float_64 const convToMassDensity = densityUnit * cleverNameYouChoose;

float_64 const convToMassDensity = densityUnit * massNumber / nAvogadro / convM3ToCM3;
float_64 const massDensity = density * convToMassDensity;

float_64 constexpr atomicTimesMassNumber = protonNumber * massNumber;
Copy link
Member

@ax3l ax3l Feb 17, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

review comment not applied fully. Please write

float_64 constexpr invAtomicTimesMassNumber = float_64(1.) / (protonNumber * massNumber);
float_64 const R = massDensity * invAtomicTimesMassNumber;

instead.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, you're right. Done.

float_64 const C = TFC1 * T_F + TFC2;

/* requires mass density in g/cm^3 */
float_64 constexpr nAvogadro = SI::N_AVOGADRO;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please keep constexpr as the first argument. constexpr is not a type specifier which can be read right to left but is a specifier for the expression you are writing currently. discussed again today with @psychocoderHPC

example difference to const: not all arguments on the right hand side need to be const to define a const.
but: to define a constexpr variable, the expression itself needs to be a constant expression.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the clarification!
Done!

@@ -0,0 +1,77 @@
/**
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

/** -> /*

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

@@ -0,0 +1,30 @@
/**
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

/** -> /*

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

@@ -138,6 +139,10 @@ struct Hydrogen
* -> circularly polarized lasers
* - Keldysh : Keldysh ionization model
*
* - ThomasFermi : statistical impact ionization based on Thomas-Fermi
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

after rebase, please add this block in ionizerConfig.param

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done


#pragma once

/** \file ionizers.hpp
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

above #pragma once pls

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

/** \file AlgorithmThomasFermi.hpp
*
* IONIZATION ALGORITHM for the Thomas-Fermi model
*
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please remove the empty line and keep one block. had problems with adding newlines in the \file description in the last days.

Write like this please:

/** \file AlgorithmThomasFermi.hpp
 *
 * IONIZATION ALGORITHM for the Thomas-Fermi model
 * - implements the calculation of average ionization degree and changes charge states
 *   by decreasing the number of bound electrons
 * - is called with the IONIZATION MODEL, specifically by setting the flag in @see speciesDefinition.param
 */


#pragma once

/* \file ionizers.def
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please above #pragma once

Copy link
Member Author

@n01r n01r Feb 22, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done, also changed to /** \file ...

{
namespace ionization
{

Copy link
Member

@ax3l ax3l Feb 17, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

plase add a thomasFermi namespace, after rebase you see that you can skip the particles namespace here

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added a particles namespace as well - but just for the TF fitting parameters for now to not break the methods which use the atomic and effective atomic numbers.

I still had to add a lot of thomasFermi:: though. (Unless every model gets its own namespace.)

* \url http://www.sciencedirect.com/science/article/pii/S0065219908601451
* doi:10.1016/S0065-2199(08)60145-1
*/
float_X constexpr TFAlpha = 14.3139;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

keep constexpr left pls (see above)

Marco Garten added 6 commits February 21, 2017 17:56
Adds a statistical collisional/impact ionization model based on the
Thomas-Fermi atomic model.

Files: - AlgorithmThomasFermi.hpp: calculates an average charge state
       - TFFittingParameters.def: provides the numerical fitting parameters
       - ThomasFermi.def: defines the model for use in speciesDefinition.param
       - ThomasFermi_Impl.hpp: prepares necessary fields, caches data,
                               calls the calculation
Adds forward declarations and hints to implementation to prepare the model
for being chosen in the params.
Depending on the effective charge state `ZStar` which is usually a
floating point value the new charge state for the particle will be determined
by `floor(ZStar)` plus the result of a Monte-Carlo process. (0 or 1)
@n01r
Copy link
Member Author

n01r commented Feb 22, 2017

offline talk: refactor all (new) typedefs to using

- Replaced `/** Copyright` with `/* Copyright` in license headers to fix doxygen
- moved all `constexpr` to the very left again
- put the TF fitting parameters into `namespace thomasFermi`
- replaced all `typedef`s by C++11 `using`
@n01r n01r force-pushed the topic-ThomasFermiNewFieldTmp branch from 89691fa to bf1999e Compare February 22, 2017 12:59
@n01r
Copy link
Member Author

n01r commented Feb 22, 2017

offline talk: dc.releaseData the RNG used by ThomasFermi_Impl.hpp still produces a seg fault after the simulation has finished but this will be changed later in DataConnector.hpp
update #1886

Apart from that the runtime checks seem okay after rebase and changes.

@n01r
Copy link
Member Author

n01r commented Feb 22, 2017

Merrrrrr-ge me .... harrrrrr-der!

Copy link
Member

@ax3l ax3l left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wonderful, well done!

@ax3l ax3l removed the request for review from psychocoderHPC February 22, 2017 15:26
@ax3l ax3l merged commit c9fa051 into ComputationalRadiationPhysics:dev Feb 22, 2017
@n01r n01r deleted the topic-ThomasFermiNewFieldTmp branch February 22, 2017 15:28
ax3l added a commit to ax3l/picongpu that referenced this pull request Apr 27, 2017
little follow-up to ComputationalRadiationPhysics#1754 adding a missing include that caused
compile issues
ax3l added a commit to ax3l/picongpu that referenced this pull request Apr 28, 2017
little follow-up to ComputationalRadiationPhysics#1754 adding a missing include that caused
compile issues
ax3l added a commit to ax3l/picongpu that referenced this pull request May 11, 2017
little follow-up to ComputationalRadiationPhysics#1754 adding a missing include that caused
compile issues
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
component: core in PIConGPU (core application)
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants