Skip to content

Commit 1449b27

Browse files
BrianMarreBrianMarre
committed
add second branch from SCFLY
1 parent bebdda7 commit 1449b27

File tree

8 files changed

+171
-32
lines changed

8 files changed

+171
-32
lines changed

‎include/picongpu/particles/atomicPhysics/IPDModel.param‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,15 @@
2626

2727
namespace picongpu::atomicPhysics
2828
{
29+
/** switch for Stewart-Pyatt ionization potential depression model
30+
*
31+
* true =^= use minimum of additional scfly term and Stewart-Pyatt ionization potential depression,
32+
* false =^= use Stewart-Pyatt ionization potential depression directly
33+
*/
34+
constexpr bool useSCFLYextendedStewartPyattIPD = true;
35+
2936
//! configuration of used ionization potential depression Model
3037
using IPDModel = picongpu::particles::atomicPhysics::ionizationPotentialDepression::StewartPyattIPD<
31-
picongpu::particles::atomicPhysics::ionizationPotentialDepression::RelativisticTemperatureFunctor>;
38+
picongpu::particles::atomicPhysics::ionizationPotentialDepression::RelativisticTemperatureFunctor,
39+
useSCFLYextendedStewartPyattIPD>;
3240
} // namespace picongpu::atomicPhysics

‎include/picongpu/particles/atomicPhysics/debug/TestIonizationPotentialDepression.hpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -66,12 +66,14 @@ namespace picongpu::particles::atomicPhysics::debug
6666
float_64 const correctIPDValue = 6.306390823271927;
6767

6868
using StewartPyattIPD = particles::atomicPhysics::ionizationPotentialDepression::template StewartPyattIPD<
69-
particles::atomicPhysics::ionizationPotentialDepression::RelativisticTemperatureFunctor>;
69+
particles::atomicPhysics::ionizationPotentialDepression::RelativisticTemperatureFunctor,
70+
false>;
7071

7172
auto const superCellConstantInput = StewartPyattIPD::SuperCellConstantInput{
7273
static_cast<float_X>(temperatureTimesk_Boltzman),
7374
static_cast<float_X>(debyeLength),
74-
static_cast<float_X>(zStar)};
75+
static_cast<float_X>(zStar),
76+
static_cast<float_X>(electronDensity / (sim.unit.length() * sim.unit.length() * sim.unit.length()))};
7577

7678
// eV
7779
float_64 const ipd = static_cast<float_X>(StewartPyattIPD::ipd(superCellConstantInput, chargeState));

‎include/picongpu/particles/atomicPhysics/ionizationPotentialDepression/LocalIPDInputFields.hpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,5 +28,6 @@
2828
#pragma once
2929

3030
#include "picongpu/particles/atomicPhysics/ionizationPotentialDepression/localHelperFields/DebyeLengthField.hpp"
31+
#include "picongpu/particles/atomicPhysics/ionizationPotentialDepression/localHelperFields/FreeElectronDensityField.hpp"
3132
#include "picongpu/particles/atomicPhysics/ionizationPotentialDepression/localHelperFields/TemperatureEnergyField.hpp"
3233
#include "picongpu/particles/atomicPhysics/ionizationPotentialDepression/localHelperFields/ZStarField.hpp"

‎include/picongpu/particles/atomicPhysics/ionizationPotentialDepression/StewartPyattIPD.hpp‎

Lines changed: 75 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -48,14 +48,18 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression
4848
float_X debyeLength;
4949
// unitless
5050
float_X zStar;
51+
// 1/unit_length^3
52+
float_X freeElectronDensity;
5153

5254
constexpr StewartPyattSuperCellConstantInput(
53-
float_X temperatureTimesk_Boltzman_i,
54-
float_X debyeLength_i,
55-
float_X zStar_i)
55+
float_X const temperatureTimesk_Boltzman_i,
56+
float_X const debyeLength_i,
57+
float_X const zStar_i,
58+
float_X const freeElectronDensity_i)
5659
: temperatureTimesk_Boltzman(temperatureTimesk_Boltzman_i)
5760
, debyeLength(debyeLength_i)
5861
, zStar(zStar_i)
62+
, freeElectronDensity(freeElectronDensity_i)
5963
{
6064
}
6165
};
@@ -69,7 +73,7 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression
6973
* @tparam T_TemperatureFunctor term A to average over for all macro particles according to equi-partition theorem,
7074
* average(A) = k_B * T, must follow
7175
*/
72-
template<typename T_TemperatureFunctor>
76+
template<typename T_TemperatureFunctor, bool T_useSCFLYextendedStewartPyattIPD>
7377
struct StewartPyattIPD : IPDModel
7478
{
7579
using SuperCellConstantInput = detail::StewartPyattSuperCellConstantInput;
@@ -143,17 +147,23 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression
143147
= std::make_unique<s_IPD::localHelperFields::DebyeLengthField<picongpu::MappingDesc>>(mappingDesc);
144148
dataConnector.consume(std::move(debyeLengthField));
145149

150+
// local k_Boltzman * Temperature field, in eV
151+
auto temperatureEnergyField
152+
= std::make_unique<s_IPD::localHelperFields::TemperatureEnergyField<picongpu::MappingDesc>>(
153+
mappingDesc);
154+
dataConnector.consume(std::move(temperatureEnergyField));
155+
146156
// z^star IPD input field, z^star = = average(q^2) / average(q) ;for q charge number of ion, unitless,
147157
// not weighted
148158
auto zStarField
149159
= std::make_unique<s_IPD::localHelperFields::ZStarField<picongpu::MappingDesc>>(mappingDesc);
150160
dataConnector.consume(std::move(zStarField));
151161

152-
// local k_Boltzman * Temperature field, in eV
153-
auto temperatureEnergyField
154-
= std::make_unique<s_IPD::localHelperFields::TemperatureEnergyField<picongpu::MappingDesc>>(
162+
// unit: 1/unit_length^3
163+
auto freeElectronDensityField
164+
= std::make_unique<s_IPD::localHelperFields::FreeElectronDensityField<picongpu::MappingDesc>>(
155165
mappingDesc);
156-
dataConnector.consume(std::move(temperatureEnergyField));
166+
dataConnector.consume(std::move(freeElectronDensityField));
157167
//@}
158168
}
159169

@@ -205,7 +215,7 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression
205215
T_AtomicPhysicsIonSpeciesList,
206216
s_IPD::stage::ApplyIPDIonization<
207217
boost::mpl::_1,
208-
StewartPyattIPD<T_TemperatureFunctor>,
218+
StewartPyattIPD<T_TemperatureFunctor, T_useSCFLYextendedStewartPyattIPD>,
209219
std::integral_constant<bool, T_SkipFinishedSuperCell>>>;
210220
ForEachIonSpeciesApplyIPDIonization{}(mappingDesc);
211221
};
@@ -220,21 +230,28 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression
220230
* @param zStarBox deviceDataBox giving access to the local z^Star value, = average(q^2) / average(q),
221231
* for all local superCells, unitless, not weighted
222232
*/
223-
template<typename T_DebyeLengthBox, typename T_TemperatureEnergyBox, typename T_ZStarBox>
233+
template<
234+
typename T_DebyeLengthBox,
235+
typename T_TemperatureEnergyBox,
236+
typename T_ZStarBox,
237+
typename T_FreeElectronDensityBox>
224238
HDINLINE static SuperCellConstantInput getSuperCellConstantInput(
225239
pmacc::DataSpace<simDim> const superCellFieldIdx,
226240
T_DebyeLengthBox const debyeLengthBox,
227241
T_TemperatureEnergyBox const temperatureEnergyBox,
228-
T_ZStarBox const zStarBox)
242+
T_ZStarBox const zStarBox,
243+
T_FreeElectronDensityBox const freeElectronDensityBox)
229244
{
230245
// eV, not weighted
231246
float_X temperatureTimesk_Boltzman = temperatureEnergyBox(superCellFieldIdx);
232247
// UNIT_LENGTH, not weighted
233248
float_X debyeLength = debyeLengthBox(superCellFieldIdx);
234249
// unitless, not weighted
235250
float_X zStar = zStarBox(superCellFieldIdx);
251+
// 1/unit_length^3
252+
float_X freeElectronDensity = freeElectronDensityBox(superCellFieldIdx);
236253

237-
return SuperCellConstantInput(temperatureTimesk_Boltzman, debyeLength, zStar);
254+
return SuperCellConstantInput(temperatureTimesk_Boltzman, debyeLength, zStar, freeElectronDensity);
238255
}
239256

240257
/** calculate ionization potential depression
@@ -255,20 +272,57 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression
255272
// UNIT_ENERGY/eV
256273
constexpr float_X eV = sim.pic.get_eV();
257274

275+
float_X const chargeState_X = static_cast<float_X>(chargeState);
276+
258277
// UNIT_MASS * UNIT_LENGTH^3 / UNIT_TIME^2 * 1/(eV * UNIT_ENERGY/eV * UNIT_LENGTH)
259278
// = unitless, not weighted
260279
float_X const K
261280
= (pmacc::math::isApproxZero(
262281
superCellConstantInput.temperatureTimesk_Boltzman * superCellConstantInput.debyeLength))
263282
? 0._X
264-
: constFactor * (chargeState + 1)
283+
: constFactor * (chargeState_X + 1._X)
265284
/ (superCellConstantInput.temperatureTimesk_Boltzman * eV
266285
* superCellConstantInput.debyeLength);
267286

268287
// eV, not weighted
269-
return superCellConstantInput.temperatureTimesk_Boltzman
270-
* (math::pow((3 * (superCellConstantInput.zStar + 1) * K + 1), 2._X / 3._X) - 1._X)
271-
/ (2._X * (superCellConstantInput.zStar + 1._X));
288+
float_X const stewartPyattIPD
289+
= superCellConstantInput.temperatureTimesk_Boltzman
290+
* (math::pow((3._X * (superCellConstantInput.zStar + 1._X) * K + 1._X), 2._X / 3._X) - 1._X)
291+
/ (2._X * (superCellConstantInput.zStar + 1._X));
292+
293+
float_X result = 0._X;
294+
295+
// additional ipd contribution as implemented in scfly
296+
if constexpr(T_useSCFLYextendedStewartPyattIPD)
297+
{
298+
/** @details taken directly from the scfly source code,
299+
* https://github.com/ComputationalRadiationPhysics/scfly/blob/5d9b0a997bb6688b04e1a34a7fa8db0f2657106a/code/src/scdrv.f#L3519-L3547
300+
*
301+
* no source known, Brian Marre, 2025
302+
*/
303+
// (1/(1/unit_length^3))^(1/3) = unit_length
304+
float_X const ionSphereRadius = math::pow(
305+
0.75_X * (chargeState_X + 1._X) / (picongpu::PI * superCellConstantInput.freeElectronDensity),
306+
1. / 3.);
307+
308+
// (eV * cm) * m/cm * unit_length/m = eV * unit_length
309+
constexpr float_X constFactorSCFLY = 1.45e-7_X / sim.unit.length() * 1.e-2_X;
310+
311+
// 1 * eV * unit_length / sqrt(unit_length^2 + unit_length^2) = eV, not weighted
312+
float_X const scfly_secondBranch
313+
= (chargeState_X + 1._X) * constFactorSCFLY
314+
/ math::sqrt(
315+
(ionSphereRadius / 1.5_X) * (ionSphereRadius / 1.5_X)
316+
+ (superCellConstantInput.debyeLength) * (superCellConstantInput.debyeLength));
317+
318+
result = math::min(scfly_secondBranch, stewartPyattIPD);
319+
}
320+
else
321+
{
322+
result = stewartPyattIPD;
323+
}
324+
325+
return result;
272326
}
273327

274328
/** call a kernel, appending the IPD-input to the kernel's input
@@ -290,14 +344,18 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression
290344
= *dc.get<s_IPD::localHelperFields::TemperatureEnergyField<picongpu::MappingDesc>>(
291345
"TemperatureEnergyField");
292346
auto& zStarField = *dc.get<s_IPD::localHelperFields::ZStarField<picongpu::MappingDesc>>("ZStarField");
347+
auto& freeElectronDensityField
348+
= *dc.get<s_IPD::localHelperFields::FreeElectronDensityField<picongpu::MappingDesc>>(
349+
"FreeElectronDensityField");
293350

294351
PMACC_LOCKSTEP_KERNEL(T_Kernel())
295352
.template config<T_chunkSize>(mapper.getGridDim())(
296353
mapper,
297354
kernelInput...,
298355
debyeLengthField.getDeviceDataBox(),
299356
temperatureEnergyField.getDeviceDataBox(),
300-
zStarField.getDeviceDataBox());
357+
zStarField.getDeviceDataBox(),
358+
freeElectronDensityField.getDeviceDataBox());
301359
}
302360
};
303361
} // namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression

‎include/picongpu/particles/atomicPhysics/ionizationPotentialDepression/kernel/CalculateIPDInput.kernel‎

Lines changed: 19 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,8 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::ker
7070
* for all local superCells, unitless, not weighted
7171
* @param debyeLengthBox deviceDataBox giving access to the local debye length for all local superCells,
7272
* sim.unit.length(), not weighted
73+
* @param freeElectronDensityBox deviceDataBox giving access to the local free electron density for all local
74+
* superCells, unit 1/unit_length^3, not weighted
7375
*/
7476
template<
7577
typename T_Worker,
@@ -82,7 +84,8 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::ker
8284
typename T_SumChargeNumberSquaredIonsFieldDataBox,
8385
typename T_LocalTemperatureFieldDataBox,
8486
typename T_ZStarFieldDataBox,
85-
typename T_LocaDebyeLengthFieldDataBox>
87+
typename T_LocaDebyeLengthFieldDataBox,
88+
typename T_FreeElectronDensityFieldDataBox>
8689
HDINLINE void operator()(
8790
T_Worker const& worker,
8891
T_AreaMapping const areaMapping,
@@ -92,9 +95,10 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::ker
9295
T_SumWeightElectronFieldBox const sumWeightElectronNormalizedBox,
9396
T_SumChargeNumberIonsFieldDataBox const sumChargeNumberBox,
9497
T_SumChargeNumberSquaredIonsFieldDataBox const sumChargeNumberSquaredBox,
98+
T_LocaDebyeLengthFieldDataBox debyeLengthBox,
9599
T_LocalTemperatureFieldDataBox temperatureEnergyBox,
96100
T_ZStarFieldDataBox zStarBox,
97-
T_LocaDebyeLengthFieldDataBox debyeLengthBox) const
101+
T_FreeElectronDensityFieldDataBox freeElectronDensityBox) const
98102
{
99103
auto const superCellFieldIdx = KernelIndexation::getSuperCellFieldIndex(worker, areaMapping);
100104

@@ -130,8 +134,9 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::ker
130134
&sumChargeNumber,
131135
&sumChargeNumberSquared,
132136
&temperatureEnergyBox,
133-
&zStarBox,
134137
&debyeLengthBox,
138+
&zStarBox,
139+
&freeElectronDensityBox,
135140
&superCellFieldIdx]()
136141
{
137142
/* unit: (unitless * weight / typicalNumParticlesPerMacroParticle)
@@ -175,12 +180,17 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::ker
175180
* unit: sqrt(unit_length^2 / unit_energy / typicalNumParticlesPerMacroParticle
176181
* * unit_energy / (1/typicalNumParticlesPerMacroParticle)) = unit_length */
177182
debyeLengthBox(superCellFieldIdx)
178-
= (!pmacc::math::isApproxZero(sumChargeNumberSquared + sumWeightElectronNormalized))
179-
? math::sqrt(
180-
constFactorDebyeLength * temperatureEnergy_PIC
181-
/ (sumChargeNumberSquared + sumWeightElectronNormalized))
182-
: -1._X;
183-
;
183+
= ((!pmacc::math::isApproxZero(sumChargeNumberSquared + sumWeightElectronNormalized))
184+
? math::sqrt(
185+
constFactorDebyeLength * temperatureEnergy_PIC
186+
/ (sumChargeNumberSquared + sumWeightElectronNormalized))
187+
: -1._X);
188+
189+
constexpr float_X constFactorFreeElectronDensity = sim.unit.typicalNumParticlesPerMacroParticle();
190+
191+
// unit: 1 / unit_length^3
192+
freeElectronDensityBox(superCellFieldIdx)
193+
= constFactorFreeElectronDensity * sumWeightElectronNormalized / volumeSuperCell;
184194
});
185195
}
186196
};

‎include/picongpu/particles/atomicPhysics/ionizationPotentialDepression/localHelperFields/FreeElectronDensityField.hpp‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
/* Copyright 2024-2024 Brian Marre
2+
*
3+
* This file is part of PIConGPU.
4+
*
5+
* PIConGPU is free software: you can redistribute it and/or modify
6+
* it under the terms of the GNU General Public License as published by
7+
* the Free Software Foundation, either version 3 of the License, or
8+
* (at your option) any later version.
9+
*
10+
* PIConGPU is distributed in the hope that it will be useful,
11+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
12+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13+
* GNU General Public License for more details.
14+
*
15+
* You should have received a copy of the GNU General Public License
16+
* along with PIConGPU.
17+
* If not, see <http://www.gnu.org/licenses/>.
18+
*/
19+
20+
#pragma once
21+
22+
#include "picongpu/defines.hpp"
23+
#include "picongpu/particles/atomicPhysics/SuperCellField.hpp"
24+
25+
#include <cstdint>
26+
#include <string>
27+
28+
namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::localHelperFields
29+
{
30+
/**superCell field of the local free electron density
31+
*
32+
* @details unitless, not weighted
33+
*
34+
* @note required for calculating the local ionization potential depression(IPD) and filled by
35+
* calculateIPDInput kernel.
36+
*
37+
* @tparam T_MappingDescription description of local mapping from device to grid
38+
*/
39+
template<typename T_MappingDescription>
40+
struct FreeElectronDensityField : public SuperCellField<float_X, T_MappingDescription, /*no guards*/ false>
41+
{
42+
FreeElectronDensityField(T_MappingDescription const& mappingDesc)
43+
: SuperCellField<float_X, T_MappingDescription, /*no guards*/ false>(mappingDesc)
44+
{
45+
}
46+
47+
// required by ISimulationData
48+
std::string getUniqueId() override
49+
{
50+
return "FreeElectronDensityField";
51+
}
52+
};
53+
} // namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::localHelperFields

‎include/picongpu/particles/atomicPhysics/ionizationPotentialDepression/stage/ApplyIPDIonization.hpp‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,8 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::sta
7878
= *dc.get<s_IPD::localHelperFields::TemperatureEnergyField<picongpu::MappingDesc>>(
7979
"TemperatureEnergyField");
8080
auto& zStarField = *dc.get<s_IPD::localHelperFields::ZStarField<picongpu::MappingDesc>>("ZStarField");
81+
auto& freeElectronDensityField
82+
= *dc.get<s_IPD::localHelperFields::ZStarField<picongpu::MappingDesc>>("FreeElectronDensityField");
8183
auto idProvider = dc.get<IdProvider>("globalId");
8284

8385
auto& fieldE = *dc.get<FieldE>(FieldE::getName());
@@ -104,7 +106,8 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::sta
104106
fieldE.getDeviceDataBox(),
105107
debyeLengthField.getDeviceDataBox(),
106108
temperatureEnergyField.getDeviceDataBox(),
107-
zStarField.getDeviceDataBox());
109+
zStarField.getDeviceDataBox(),
110+
freeElectronDensityField.getDeviceDataBox());
108111
}
109112
else
110113
{
@@ -126,7 +129,8 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::sta
126129
fieldE.getDeviceDataBox(),
127130
debyeLengthField.getDeviceDataBox(),
128131
temperatureEnergyField.getDeviceDataBox(),
129-
zStarField.getDeviceDataBox());
132+
zStarField.getDeviceDataBox(),
133+
freeElectronDensityField.getDeviceDataBox());
130134
}
131135

132136
// no need to call fillAllGaps, since we do not leave any gaps

0 commit comments

Comments
 (0)