Skip to content

Commit

Permalink
update EField
Browse files Browse the repository at this point in the history
ci: picongpu
  • Loading branch information
BrianMarre committed Nov 5, 2024
1 parent 7bd0120 commit a062878
Show file tree
Hide file tree
Showing 4 changed files with 217 additions and 0 deletions.
Original file line number Diff line number Diff line change
@@ -0,0 +1,130 @@
/* Copyright 2024 Brian Marre
*
* This file is part of PIConGPU.
*
* PIConGPU is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* PIConGPU is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with PIConGPU.
* If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once

#include "picongpu/defines.hpp"
#include "picongpu/particles/atomicPhysics/KernelIndexation.hpp"

#include <pmacc/lockstep/ForEach.hpp>

#include <cstdint>

namespace picongpu::particles::atomicPhysics::kernel
{
/** update electric field due to field ionization energy use
*
* @tparam T_numberAtomicPhysicsIonSpecies specialization template parameter used to prevent compilation of all
* atomicPhysics kernels if no atomic physics species is present.
*/
template<uint32_t T_numberAtomicPhysicsIonSpecies>
struct UpdateElectricFieldKernel
{
/** call operator
*
* called by UpdateElectricField atomic physics sub-stage
*
* @param worker object containing the device and block information, passed by PMACC_KERNEL call
* @param areMapping mapping of blockIndex to block superCell index
* @param timeRemainingBox deviceDataBox giving access to the time remaining in the atomicPhysics step for each
* local superCell
* @param eFieldBox deviceDataBox giving access to eField values for all local superCells
* @param fieldEnergyUseCacheBox deviceDataBox giving access to the field energy use cache for each local
* superCell
*/
template<
typename T_Worker,
typename T_AreaMapping,
typename T_TimeRemainingBox,
typename T_EFieldBox,
typename T_FieldEnergyUseCacheBox>
HDINLINE void operator()(
T_Worker const& worker,
T_AreaMapping const areaMapping,
T_TimeRemainingBox const timeRemainingBox,
T_EFieldBox eFieldBox,
T_FieldEnergyUseCacheBox const fieldEnergyUseCacheBox) const
{
auto const superCellIdx = KernelIndexation::getSuperCellIndex(worker, areaMapping);
auto const superCellFieldIdx = KernelIndexation::getSuperCellFieldIndex(worker, areaMapping, superCellIdx);

auto const timeRemaining = timeRemainingBox(superCellFieldIdx);
if(timeRemaining <= 0._X)
return;

T_FieldEnergyUseCache const& eFieldEnergyUseCache = fieldEnergyUseCacheBox(superCellFieldIdx);
DataSpace<picongpu::simDim> const superCellCellOffset = superCellIdx * picongpu::SuperCellSize::toRT();

auto forEachCell = pmacc::lockstep::makeForEach<T_FieldEnergyUseCache::numberCells, T_Worker>(worker);
forEachCell(
[&worker,
&superCellCellOffset,
&eFieldBox,
&eFieldEnergyUseCache,
&sharedResourcesOverSubscribed,
&rejectionProbabilityCache](uint32_t const linearCellIndex)
{
DataSpace<picongpu::simDim> const localCellIndex
= pmacc::math::mapToND(picongpu::SuperCellSize::toRT(), static_cast<int>(linearCellIndex));
DataSpace<picongpu::simDim> const cellIndex = localCellIndex + superCellCellOffset;

// ((sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge()))^2
auto const oldEFieldVector = eFieldBox(cellIndex);
auto const eFieldNormSquared = pmacc::math::l2norm2(oldEFieldVector);

/// @note zero field also means at most zero energy use
if(eFieldNormSquared <= 0)
return;

// eV * 1 = eV * sim.unit.energy()/sim.unit.energy() = (ev / sim.unit.energy()) * sim.unit.energy()
// sim.unit.energy()
float_X const eFieldEnergyUse
= picongpu::sim.pic.get_eV() * eFieldEnergyUseCache.energyUsed(linearCellIndex);

// sim.unit.charge()^2 * sim.unit.time()^2 / (sim.unit.mass() * sim.unit.length()^3)
// * sim.unit.length()^3
// = sim.unit.charge()^2 * sim.unit.time()^2 / sim.unit.mass()
constexpr float_X eps0HalfTimesCellVolume
= (picongpu::sim.pic.getEps0() / 2._X) * picongpu::sim.pic.getCellSize().productOfComponents();

/* (((sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge()))^2
* - sim.unit.energy() / (sim.unit.charge()^2 * sim.unit.time()^2 / sim.unit.mass()))^1/2
* = (((sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge()))^2
* - (sim.unit.mass() * sim.unit.length()^2/sim.unit.time()^2
* / (sim.unit.charge()^2 * sim.unit.time()^2 / sim.unit.mass()))
* )^1/2
* = (((sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge()))^2
* - (sim.unit.mass()^2 * sim.unit.length()^2/(sim.unit.time()^4 * sim.unit.charge()^2))
* )^1/2
* = ((sim.unit.mass()^2 * sim.unit.length()^2)/(sim.unit.time()^4 * sim.unit.charge()^2)
* - (sim.unit.mass()^2 * sim.unit.length()^2/(sim.unit.time()^4 * sim.unit.charge()^2))
* )^1/2
* = ((sim.unit.mass()^2 * sim.unit.length()^2)/(sim.unit.time()^4 * sim.unit.charge()^2))^1/2
* = (sim.unit.mass() * sim.unit.length())/(sim.unit.time()^2 * sim.unit.charge())
* = sim.unit.eField()
*/
// sim.unit.eField()
float_X const newEFieldNorm
= pmacc::math::sqrt(eFieldNormSquared - eFieldEnergyUse / eps0HalfTimesCellVolume);

eFieldBox(cellIndex) = oldEFieldVector / pmacc::math::sqrt(eFieldNormSquared) * newEFieldNorm;
});
}
};
} // namespace picongpu::particles::atomicPhysics::kernel
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
#include "picongpu/fields/FieldE.hpp"
#include "picongpu/particles/atomicPhysics/electronDistribution/LocalHistogramField.hpp"
#include "picongpu/particles/atomicPhysics/kernel/CheckForOverSubscription.kernel"
#include "picongpu/particles/atomicPhysics/localHelperFields/FieldEnergyUseCacheField.hpp"
#include "picongpu/particles/atomicPhysics/localHelperFields/RejectionProbabilityCacheField.hpp"
#include "picongpu/particles/atomicPhysics/localHelperFields/SharedResourcesOverSubscribedField.hpp"
#include "picongpu/particles/atomicPhysics/localHelperFields/TimeRemainingField.hpp"
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
/* Copyright 2024 Brian Marre
*
* This file is part of PIConGPU.
*
* PIConGPU is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* PIConGPU is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with PIConGPU.
* If not, see <http://www.gnu.org/licenses/>.
*/

//! @file update electric field with field ionization energy use

#pragma once

// need picongpu::atomicPhysics::ElectronHistogram from atomicPhysics.param
#include "picongpu/defines.hpp"
#include "picongpu/fields/FieldE.hpp"
#include "picongpu/particles/atomicPhysics/kernel/UpdateElectricField.kernel"
#include "picongpu/particles/atomicPhysics/localHelperFields/FieldEnergyUseCacheField.hpp"
#include "picongpu/particles/atomicPhysics/localHelperFields/TimeRemainingField.hpp"

#include <pmacc/Environment.hpp>
#include <pmacc/mappings/kernel/AreaMapping.hpp>
#include <pmacc/type/Area.hpp>

#include <cstdint>

namespace picongpu::particles::atomicPhysics::stage
{
/** UpdateElectricField atomic physics sub-stage
*
* remove the used field energy from the electric field by reducing the norm of the local E-Field vector to match
* the energy use
*
* @tparam T_numberAtomicPhysicsIonSpecies specialization template parameter used to prevent compilation of all
* atomicPhysics kernels if no atomic physics species is present.
*/
template<uint32_t T_numberAtomicPhysicsIonSpecies>
struct UpdateElectricField
{
//! call of kernel for every superCell
HINLINE void operator()(picongpu::MappingDesc const mappingDesc) const
{
// full local domain, no guards
pmacc::AreaMapping<CORE + BORDER, MappingDesc> mapper(mappingDesc);
pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector();

auto& timeRemainingField = *dc.get<
picongpu::particles::atomicPhysics::localHelperFields::TimeRemainingField<picongpu::MappingDesc>>(
"TimeRemainingField");

auto& eField = *dc.get<FieldE>(FieldE::getName());

using FieldEnergyUseCacheField
= picongpu::particles::atomicPhysics::localHelperFields::FieldEnergyUseCacheField<
picongpu::MappingDesc>;
auto& fieldEnergyUseCacheField = *dc.get<FieldEnergyUseCacheField>("FieldEnergyUseCacheField");

// macro for call of kernel for every superCell, see pull request #4321
PMACC_LOCKSTEP_KERNEL(
picongpu::particles::atomicPhysics::kernel::UpdateElectricField<T_numberAtomicPhysicsIonSpecies>())
.template config<FieldEnergyUseCacheField::ElementType::numberCells>(mapper.getGridDim())(
mapper,
timeRemainingField.getDeviceDataBox(),
eField.getDeviceDataBox(),
fieldEnergyUseCacheField.getDeviceDataBox());
}
};
} // namespace picongpu::particles::atomicPhysics::stage
8 changes: 8 additions & 0 deletions include/picongpu/simulation/stage/AtomicPhysics.x.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@
#include "picongpu/particles/atomicPhysics/stage/ResetTimeStepField.hpp"
#include "picongpu/particles/atomicPhysics/stage/RollForOverSubscription.hpp"
#include "picongpu/particles/atomicPhysics/stage/SpawnIonizationElectrons.hpp"
#include "picongpu/particles/atomicPhysics/stage/UpdateElectricField.hpp"
#include "picongpu/particles/atomicPhysics/stage/UpdateTimeRemaining.hpp"
#include "picongpu/particles/filter/filter.hpp"

Expand Down Expand Up @@ -440,6 +441,12 @@ namespace picongpu::simulation::stage
ForEachIonSpeciesSpawnIonizationElectrons{}(mappingDesc, currentStep);
}

HINLINE static void updateElectricField(picongpu::MappingDesc const& mappingDesc)
{
picongpu::particles::atomicPhysics::stage::UpdateElectricField<T_numberAtomicPhysicsIonSpecies>()(
mappingDesc);
}

template<typename T_DeviceReduce>
HINLINE static void doIPDIonization(
picongpu::MappingDesc const& mappingDesc,
Expand Down Expand Up @@ -576,6 +583,7 @@ namespace picongpu::simulation::stage

recordChanges(mappingDesc);
updateElectrons(mappingDesc, currentStep);
updateElectricField(mappingDesc);
doIPDIonization(mappingDesc, currentStep, deviceLocalReduce);
updateTimeRemaining(mappingDesc);
isSubSteppingComplete = isSubSteppingFinished(mappingDesc, deviceLocalReduce);
Expand Down

0 comments on commit a062878

Please sign in to comment.