From c9a41b4918fe500161cac7814a8d72df0e9cf680 Mon Sep 17 00:00:00 2001 From: Edgar Marquardt Date: Fri, 24 Jul 2026 09:42:46 +0200 Subject: [PATCH 1/7] add poisson solver to generate electrical field from non neutral charged gas densities MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit With this commit PIConGPU provides a poisson solver to generate the electric field during the starting conditions. It can be activated by providing the `--poisson.activate` flag when starting a picongpu simulation. See TBG docu and `picongpu --help` for more options. A test for the Poisson solver is also provided. Co-authored-by: René Widera Co-authered-by: Luca --- docs/TBG_macros.cfg | 11 + include/picongpu/fields/Fields.def | 4 + include/picongpu/fields/Fields.hpp | 1 + .../fields/poissonSolver/BICGStab.hpp | 38 + .../poissonSolver/BoundaryConditions.hpp | 116 +++ .../fields/poissonSolver/ChargeDeposition.hpp | 20 + .../picongpu/fields/poissonSolver/FieldV.hpp | 121 ++++ .../RightHandSideNormalization.hpp | 95 +++ .../picongpu/fields/poissonSolver/Stencil.hpp | 159 ++++ .../simulation/control/Simulation.hpp | 10 + include/picongpu/simulation/stage/Poisson.hpp | 104 +++ .../picongpu/simulation/stage/Poisson.x.cpp | 680 ++++++++++++++++++ include/pmacc/algorithms/ForEachCell.hpp | 108 +++ share/picongpu/tests/PoissonSolver/README.rst | 10 + share/picongpu/tests/PoissonSolver/bin/ci.sh | 157 ++++ .../include/picongpu/param/density.param | 32 + .../include/picongpu/param/dimension.param | 31 + .../include/picongpu/param/particle.param | 70 ++ .../include/picongpu/param/simulation.param | 84 +++ .../picongpu/param/speciesDefinition.param | 91 +++ .../param/speciesInitialization.param | 47 ++ .../lib/python/test/validate_results.py | 117 +++ 22 files changed, 2106 insertions(+) create mode 100644 include/picongpu/fields/poissonSolver/BICGStab.hpp create mode 100644 include/picongpu/fields/poissonSolver/BoundaryConditions.hpp create mode 100644 include/picongpu/fields/poissonSolver/ChargeDeposition.hpp create mode 100644 include/picongpu/fields/poissonSolver/FieldV.hpp create mode 100644 include/picongpu/fields/poissonSolver/RightHandSideNormalization.hpp create mode 100644 include/picongpu/fields/poissonSolver/Stencil.hpp create mode 100644 include/picongpu/simulation/stage/Poisson.hpp create mode 100644 include/picongpu/simulation/stage/Poisson.x.cpp create mode 100644 include/pmacc/algorithms/ForEachCell.hpp create mode 100644 share/picongpu/tests/PoissonSolver/README.rst create mode 100755 share/picongpu/tests/PoissonSolver/bin/ci.sh create mode 100644 share/picongpu/tests/PoissonSolver/include/picongpu/param/density.param create mode 100644 share/picongpu/tests/PoissonSolver/include/picongpu/param/dimension.param create mode 100644 share/picongpu/tests/PoissonSolver/include/picongpu/param/particle.param create mode 100644 share/picongpu/tests/PoissonSolver/include/picongpu/param/simulation.param create mode 100644 share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesDefinition.param create mode 100644 share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesInitialization.param create mode 100755 share/picongpu/tests/PoissonSolver/lib/python/test/validate_results.py diff --git a/docs/TBG_macros.cfg b/docs/TBG_macros.cfg index facbf957a11..db2f7c00bf8 100644 --- a/docs/TBG_macros.cfg +++ b/docs/TBG_macros.cfg @@ -113,6 +113,17 @@ TBG_particleBoundaries="--e_boundary periodic absorbing thermal --e_boundaryOffs # In case the species has multiple such FromOpenPMDImpl invocations, same runtime value is used for all of them. TBG_runtimeDensityFile="--e_runtimeDensityFile /bigdata/hplsim/production/someDirectory/density.bp" +# Activate the Poisson solver and set its parameters. The Poisson solver is used to calculate the electrostatic +# potential from the charge density. +# --poisson.maxSteps sets the maximum number of iterations for the solver, and --poisson.epsilon sets the convergence +# criterion: the maximum acceptable relative error. +TBG_poissonSolver="--poisson.activate --poisson.maxSteps 2000 --poisson.epsilon 1e-8" + +# The Poisson solver is by default preconditioned. The preconditioner can be disabled with --poisson.preconditioner.disable. +# The maximum number of iterations for the preconditioner can be set with --poisson.preconditioner.maxSteps. +# The preconditioner is used to improve the convergence of the Poisson solver. +TBG_poissonSolverPreconditioner="--poisson.preconditioner.disable --poisson.preconditioner.maxSteps 20" + # Set absorber type of absorbing boundaries. # Supported values: exponential, pml (default). diff --git a/include/picongpu/fields/Fields.def b/include/picongpu/fields/Fields.def index 2744b6aac71..777c6e94d83 100644 --- a/include/picongpu/fields/Fields.def +++ b/include/picongpu/fields/Fields.def @@ -68,4 +68,8 @@ namespace picongpu /** Current Density j, @see FieldJ.hpp */ class FieldJ; + namespace fields::poissonSolver + { + struct FieldV; + } // namespace fields::poissonSolver } // namespace picongpu diff --git a/include/picongpu/fields/Fields.hpp b/include/picongpu/fields/Fields.hpp index 8ab7dd93050..dee16e1fe7f 100644 --- a/include/picongpu/fields/Fields.hpp +++ b/include/picongpu/fields/Fields.hpp @@ -25,3 +25,4 @@ #include "picongpu/fields/FieldJ.hpp" #include "picongpu/fields/FieldTmp.hpp" #include "picongpu/fields/Fields.def" +#include "picongpu/fields/poissonSolver/FieldV.hpp" diff --git a/include/picongpu/fields/poissonSolver/BICGStab.hpp b/include/picongpu/fields/poissonSolver/BICGStab.hpp new file mode 100644 index 00000000000..85ef7c092f5 --- /dev/null +++ b/include/picongpu/fields/poissonSolver/BICGStab.hpp @@ -0,0 +1,38 @@ +/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera + * + * 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 . + */ + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/fields/FieldTmpOperations.hpp" + +namespace picongpu::fields::poissonSolver +{ + struct BICGStab + { + // return residual + // return number of iterations + void operator()(FieldTmp& fieldV, FieldTmp& fiedlRho, MappingDesc* cellDescription) + { + // set boundary conditions on fieldV (Dirichlet or Neuman) + + // normalize the problem based on norm(fieldRho) + } + }; +} // namespace picongpu::fields::poissonSolver diff --git a/include/picongpu/fields/poissonSolver/BoundaryConditions.hpp b/include/picongpu/fields/poissonSolver/BoundaryConditions.hpp new file mode 100644 index 00000000000..085c4036856 --- /dev/null +++ b/include/picongpu/fields/poissonSolver/BoundaryConditions.hpp @@ -0,0 +1,116 @@ +/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera + * + * 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 . + */ + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/fields/FieldTmpOperations.hpp" +#include "picongpu/fields/poissonSolver/FieldV.hpp" + +#include +#include +#include + +namespace picongpu::fields::poissonSolver +{ + struct SolutionFunction + { + HDINLINE auto operator()(math::Vector const& totalCellCoordinate) const + { + return math::sin(totalCellCoordinate.x()) + math::cos(totalCellCoordinate.y()) + + 3.0 * math::sin(totalCellCoordinate.z()) + + totalCellCoordinate.x() * totalCellCoordinate.productOfComponents() + 10.0; + } + + HDINLINE auto operator()(math::Vector const& totalCellCoordinate) const + { + return math::sin(totalCellCoordinate.x()) + math::cos(totalCellCoordinate.y()) + + totalCellCoordinate.x() * totalCellCoordinate.productOfComponents() + 10.0; + } + }; + + struct ApplyDirichletBCsFromFunctionKernel + { + DINLINE auto operator()( + auto const& worker, + auto fieldVBox, + auto const boundaryFunction, + DataSpace cellOffsetToTotalOrigin, + auto const mapper) const -> void + { + // including guards + DataSpace const superCellIdx(mapper.getSuperCellIndex(worker.blockDomIdxND())); + + DataSpace numGuardCells = mapper.getGuardingSuperCells() * SuperCellSize::toRT(); + + // no guards included + DataSpace superCellTotalCellOffset + = cellOffsetToTotalOrigin + superCellIdx * SuperCellSize::toRT() - numGuardCells; + + constexpr uint32_t cellsPerSuperCell = pmacc::math::CT::volume::type::value; + + auto forEachCellInSupercell = lockstep::makeForEach(worker); + + forEachCellInSupercell( + [&](int32_t const linearCellIdx) + { + /* cell index within the superCell */ + DataSpace const cellIdx = pmacc::math::mapToND(SuperCellSize::toRT(), linearCellIdx); + // without guards + DataSpace const totalCellIdx = superCellTotalCellOffset + cellIdx; + + auto totalDistance = precisionCast(totalCellIdx) + * precisionCast(sim.pic.getCellSize().shrink()); + + fieldVBox(superCellIdx * SuperCellSize::toRT() + cellIdx) = boundaryFunction(totalDistance); + }); + } + }; + + struct BoundaryConditionsDirichlet + { + // return residual + // return number of iterations + void operator()(FieldV& fieldV, MappingDesc cellDescription) const + { + SubGrid const& subGrid = Environment::get().SubGrid(); + auto globalDomain = subGrid.getGlobalDomain(); + auto localDomain = subGrid.getLocalDomain(); + + auto cellOffsetToTotalOrigin = globalDomain.offset + localDomain.offset; + + + for(uint32_t i = 1; i < NumberOfExchanges::value; ++i) + { + /* only call for planes: left right top bottom back front*/ + if(FRONT % i == 0 && !(Environment::get().GridController().getCommunicationMask().isSet(i))) + { + ExchangeMapping mapper(cellDescription, i); + + PMACC_LOCKSTEP_KERNEL(ApplyDirichletBCsFromFunctionKernel{}) + .config(mapper.getGridDim(), SuperCellSize{})( + fieldV.fieldVBuffer->getDeviceBuffer().getDataBox(), + SolutionFunction{}, + cellOffsetToTotalOrigin, + mapper); + } + } + } + }; +} // namespace picongpu::fields::poissonSolver diff --git a/include/picongpu/fields/poissonSolver/ChargeDeposition.hpp b/include/picongpu/fields/poissonSolver/ChargeDeposition.hpp new file mode 100644 index 00000000000..e6b0bda3e46 --- /dev/null +++ b/include/picongpu/fields/poissonSolver/ChargeDeposition.hpp @@ -0,0 +1,20 @@ +/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera + * + * 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 . + */ + +#pragma once diff --git a/include/picongpu/fields/poissonSolver/FieldV.hpp b/include/picongpu/fields/poissonSolver/FieldV.hpp new file mode 100644 index 00000000000..24f8374c906 --- /dev/null +++ b/include/picongpu/fields/poissonSolver/FieldV.hpp @@ -0,0 +1,121 @@ +/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera + * + * 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 . + */ + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/fields/YeeCell.hpp" +#include "picongpu/traits/FieldPosition.hpp" +#include "picongpu/traits/SIBaseUnits.hpp" + +#include + +namespace picongpu::fields::poissonSolver +{ + struct FieldV : pmacc::ISimulationData + { + using ValueType = float_64; + + std::shared_ptr> fieldVBuffer; + + //! Unit type of field components + using UnitValueType = float1_64; + + using DataBoxType = pmacc::DataBox>; + + //! Number of components of ValueType, for serialization + static constexpr int numComponents = 1u; + + FieldV(picongpu::MappingDesc const mappingDesc) + : fieldVBuffer{std::make_shared>(mappingDesc.getGridLayout())} + { + } + + void synchronize() override + { + fieldVBuffer->getHostBuffer().copyFrom(fieldVBuffer->getDeviceBuffer()); + } + + //! Get the host data box for the field values + DataBoxType getHostDataBox() + { + return fieldVBuffer->getHostBuffer().getDataBox(); + } + + //! Get the device data box for the field values + DataBoxType getDeviceDataBox() + { + return fieldVBuffer->getDeviceBuffer().getDataBox(); + } + + GridLayout getGridLayout() + { + return fieldVBuffer->getGridLayout(); + } + + /** + * Return the globally unique identifier for this simulation data. + * + * @return globally unique identifier + */ + SimulationDataId getUniqueId() override + { + return "FieldV"; + } + + static std::string getName() + { + return "FieldV"; + } + + static UnitValueType getUnit() + { + return UnitValueType{sim.unit.eField() * sim.unit.length()}; + } + + static std::vector getUnitDimension() + { + /* V is in volts: V = kg * m^2 / (A * s^3) + * -> L^2 * M * T^-3 * I^-1 + */ + std::vector unitDimension(7, 0.0); + unitDimension.at(SIBaseUnits::length) = 2.0; + unitDimension.at(SIBaseUnits::mass) = 1.0; + unitDimension.at(SIBaseUnits::time) = -3.0; + unitDimension.at(SIBaseUnits::electricCurrent) = -1.0; + return unitDimension; + } + }; +} // namespace picongpu::fields::poissonSolver + +namespace picongpu::traits +{ + template<> + struct FieldPosition + { + using VectorVectorDD = ::pmacc::math::Vector const; + + HDINLINE FieldPosition() = default; + + HDINLINE VectorVectorDD operator()() const + { + return VectorVectorDD::create(floatD_X::create(0.0)); + } + }; +} // namespace picongpu::traits diff --git a/include/picongpu/fields/poissonSolver/RightHandSideNormalization.hpp b/include/picongpu/fields/poissonSolver/RightHandSideNormalization.hpp new file mode 100644 index 00000000000..a5fe7d6840a --- /dev/null +++ b/include/picongpu/fields/poissonSolver/RightHandSideNormalization.hpp @@ -0,0 +1,95 @@ +/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera + * + * 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 . + */ + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/fields/FieldTmpOperations.hpp" +#include "picongpu/fields/poissonSolver/FieldV.hpp" + +#include +#include +#include +#include + +namespace picongpu::fields::poissonSolver +{ + struct RightHandSideNormalizationKernel + { + DINLINE auto operator()(auto const& worker, auto const fieldVBox, auto fieldRhoBox, auto const mapper) const + -> void + { + DataSpace const relExchangeDir = Mask::getRelativeDirections(mapper.getExchangeType()); + DataSpace const superCellIdx(mapper.getSuperCellIndex(worker.blockDomIdxND())); + DataSpace superCellCellOffset = superCellIdx * SuperCellSize::toRT(); + + DataSpace numGuardCells = mapper.getGuardingSuperCells() * SuperCellSize::toRT(); + DataSpace adjustedSuperCellCellOffset = superCellCellOffset - (relExchangeDir * numGuardCells); + + DataSpace numCellsLocalDomain = mapper.getGridSuperCells() * SuperCellSize::toRT(); + + constexpr uint32_t cellsPerSuperCell = pmacc::math::CT::volume::type::value; + + auto forEachCellInSupercell = lockstep::makeForEach(worker); + + forEachCellInSupercell( + [&](int32_t const linearCellIdx) + { + /* cell index within the superCell */ + DataSpace const cellIdx = pmacc::math::mapToND(SuperCellSize::toRT(), linearCellIdx); + + DataSpace const localCellIdx = superCellCellOffset + cellIdx; + DataSpace const dataCellIdx = adjustedSuperCellCellOffset + cellIdx; + + for(uint32_t d = 0; d < simDim; ++d) + { + if(relExchangeDir[d] != 0 + && (localCellIdx[d] == 0u || localCellIdx[d] == numCellsLocalDomain[d] - 1)) + { + fieldRhoBox(dataCellIdx) += fieldVBox(dataCellIdx + relExchangeDir) + / (sim.pic.getCellSize()[d] * sim.pic.getCellSize()[d]); + } + } + }); + } + }; + + struct RightHandSideNormalization + { + // return residual + // return number of iterations + void operator()(FieldV& fieldV, FieldTmp& fieldRho, MappingDesc cellDescription) + { + /* only call for planes: left right top bottom back front*/ + for(uint32_t i = 1; i < NumberOfExchanges::value; ++i) + { + if(FRONT % i == 0 && !(Environment::get().GridController().getCommunicationMask().isSet(i))) + { + ExchangeMapping mapper(cellDescription, i); + + PMACC_LOCKSTEP_KERNEL(RightHandSideNormalizationKernel{}) + .config(mapper.getGridDim(), SuperCellSize{})( + fieldV.fieldVBuffer->getDeviceBuffer().getDataBox(), + fieldRho.getDeviceDataBox(), + mapper); + } + } + } + }; +} // namespace picongpu::fields::poissonSolver diff --git a/include/picongpu/fields/poissonSolver/Stencil.hpp b/include/picongpu/fields/poissonSolver/Stencil.hpp new file mode 100644 index 00000000000..e84690eb5bb --- /dev/null +++ b/include/picongpu/fields/poissonSolver/Stencil.hpp @@ -0,0 +1,159 @@ +/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera + * + * 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 . + */ + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/fields/FieldTmpOperations.hpp" +#include "picongpu/fields/poissonSolver/FieldV.hpp" + +#include +#include +#include + +namespace picongpu::fields::poissonSolver +{ + struct StencilFunc + { + HDINLINE auto operator()(auto fieldIn) const + { + constexpr auto cellSize = sim.pic.getCellSize().shrink(); + constexpr auto cellSizeSquared = cellSize * cellSize; + float_64 const fc0 = 2.0 * (1.0 / (cellSizeSquared)).sumOfComponents(); + return stencilFunction(fieldIn, fc0, cellSizeSquared); + } + + HDINLINE auto stencilFunction( + auto fieldIn, + float_64 const fc0, + pmacc::math::Vector const cellSizeSquared) const + { + constexpr DataSpace<3u> xDir(1, 0, 0); + constexpr DataSpace<3u> yDir(0, 1, 0); + constexpr DataSpace<3u> zDir(0, 0, 1); + + + return *fieldIn * fc0 - (fieldIn(-xDir) + fieldIn(xDir)) / cellSizeSquared.x() + - (fieldIn(-yDir) + fieldIn(yDir)) / cellSizeSquared.y() + - (fieldIn(-zDir) + fieldIn(zDir)) / cellSizeSquared.z(); + } + + HDINLINE auto stencilFunction( + auto fieldIn, + float_64 const fc0, + pmacc::math::Vector const cellSizeSquared) const + { + constexpr DataSpace<2u> xDir(1, 0); + constexpr DataSpace<2u> yDir(0, 1); + + return *fieldIn * fc0 - (fieldIn(-xDir) + fieldIn(xDir)) / cellSizeSquared.x() + - (fieldIn(-yDir) + fieldIn(yDir)) / cellSizeSquared.y(); + } + }; + + struct GetFieldEStencil + { + HDINLINE auto operator()(auto fieldIn) const + { + constexpr auto cellSize = sim.pic.getCellSize().shrink(); + return getFieldE(fieldIn, cellSize); + } + + HDINLINE auto getFieldE(auto fieldIn, pmacc::math::Vector const cellSize) const + { + float3_X eValue; + constexpr DataSpace<3u> xDir(1, 0, 0); + constexpr DataSpace<3u> yDir(0, 1, 0); + constexpr DataSpace<3u> zDir(0, 0, 1); + + eValue.x() = (*fieldIn - fieldIn(xDir)) / cellSize.x(); + eValue.y() = (*fieldIn - fieldIn(yDir)) / cellSize.y(); + eValue.z() = (*fieldIn - fieldIn(zDir)) / cellSize.z(); + + return eValue; + } + + HDINLINE auto getFieldE(auto fieldIn, pmacc::math::Vector const cellSize) const + { + float3_X eValue; + constexpr DataSpace<2u> xDir(1, 0); + constexpr DataSpace<2u> yDir(0, 1); + + eValue.x() = (*fieldIn - fieldIn(xDir)) / cellSize.x(); + eValue.y() = (*fieldIn - fieldIn(yDir)) / cellSize.y(); + eValue.z() = 0.0_X; + + return eValue; + } + }; + + struct Stencil + { + DINLINE auto operator()( + auto const& worker, + auto const mapper, + auto const stencilFunctor, + auto fieldOut, + auto fieldIn) const -> void + { + DataSpace const superCellIdx(mapper.getSuperCellIndex(worker.blockDomIdxND())); + DataSpace superCellCellOffset = superCellIdx * SuperCellSize::toRT(); + + using Type = typename decltype(fieldIn)::ValueType; + using BlockArea = pmacc::SuperCellDescription< + SuperCellSize, + typename pmacc::math::CT::make_Int::type, + typename pmacc::math::CT::make_Int::type>; + + constexpr uint32_t cellsPerSuperCell = pmacc::math::CT::volume::type::value; + + // use the cached buffer, beacuse I am doing multiple reads, moves the blockArea to shared memory + auto cache = pmacc::CachedBox::create<0, Type>(worker, BlockArea()); + auto buffShifted = fieldIn.shift(superCellCellOffset); + + // the thread collective is a convenience wrapper for lockstep make for each + // it deals with the guard offset, subtracts the origin offset + auto collective = pmacc::makeThreadCollective(); + + pmacc::math::operation::Assign assign; + collective(worker, assign, cache, buffShifted); + + worker.sync(); + + auto forEachCellInSupercell = lockstep::makeForEach(worker); + + forEachCellInSupercell( + [&](int32_t const linearCellIdx) + { + /* cell index within the superCell */ + DataSpace const cellIdx = pmacc::math::mapToND(SuperCellSize::toRT(), linearCellIdx); + fieldOut[superCellCellOffset + cellIdx] = stencilFunctor(cache.shift(cellIdx)); + }); + } + }; + + inline void stencil(auto mapper, auto const& functor, auto& outBuffer, auto& inBuffer) + { + // poisson stencil + PMACC_LOCKSTEP_KERNEL(fields::poissonSolver::Stencil{}) + .config( + mapper.getGridDim(), + SuperCellSize{})(mapper, functor, outBuffer.getDataBox(), inBuffer.getDataBox()); + } +} // namespace picongpu::fields::poissonSolver diff --git a/include/picongpu/simulation/control/Simulation.hpp b/include/picongpu/simulation/control/Simulation.hpp index 265dbb7daa5..0b3e0f902bd 100644 --- a/include/picongpu/simulation/control/Simulation.hpp +++ b/include/picongpu/simulation/control/Simulation.hpp @@ -53,6 +53,7 @@ #include "picongpu/simulation/stage/ParticleInit.hpp" #include "picongpu/simulation/stage/ParticleIonization.hpp" #include "picongpu/simulation/stage/ParticlePush.hpp" +#include "picongpu/simulation/stage/Poisson.hpp" #include "picongpu/simulation/stage/RuntimeDensityFile.hpp" #include "picongpu/simulation/stage/SynchrotronRadiation.hpp" #include "picongpu/versionFormat.hpp" @@ -157,6 +158,9 @@ namespace picongpu fieldBackground->registerHelp(desc); particleBoundaries.registerHelp(desc); runtimeDensityFile.registerHelp(desc); + + poissonSolver = std::make_shared(); + poissonSolver->registerHelp(desc); } void startSimulation() override @@ -343,6 +347,9 @@ namespace picongpu // initialize runtime density file paths runtimeDensityFile.init(); + // create memory for poisson solver + poissonSolver->init(*cellDescription); + // create factory for the random number generator uint32_t const userSeed = random::seed::ISeed{}(); uint32_t const seed = std::hash{}(std::to_string(userSeed)); @@ -460,6 +467,7 @@ namespace picongpu { simulation::stage::ParticleInit{}(step); (*atomicPhysics).fixAtomicStateInit(*cellDescription); + (*poissonSolver)(step); // Check Debye resolution particles::debyeLength::check(*cellDescription); } @@ -592,6 +600,8 @@ namespace picongpu InitialiserController* initialiserController{nullptr}; + std::shared_ptr poissonSolver; + std::unique_ptr cellDescription; // layout parameter diff --git a/include/picongpu/simulation/stage/Poisson.hpp b/include/picongpu/simulation/stage/Poisson.hpp new file mode 100644 index 00000000000..4b5d9d98ddd --- /dev/null +++ b/include/picongpu/simulation/stage/Poisson.hpp @@ -0,0 +1,104 @@ +/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera + * + * 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 . + */ + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/fields/FieldJ.hpp" +#include "picongpu/fields/FieldJ.kernel" +#include "picongpu/fields/FieldTmpOperations.hpp" +#include "picongpu/fields/currentDeposition/Deposit.hpp" +#include "picongpu/fields/poissonSolver/FieldV.hpp" +#include "picongpu/particles/filter/filter.hpp" +#include "picongpu/particles/param.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +namespace picongpu::simulation::stage +{ + //! Functor for the stage of the PIC loop performing charge deposition + struct Poisson + { + Poisson(); + + void init(picongpu::MappingDesc const mappingDesc); + + void registerHelp(po::options_description& desc); + + /** Compute the current created by particles and add it to the current + * density + * + * @param currentStep index of time iteration + */ + void operator()(uint32_t const currentStep); + + void preconditioner( + std::unique_ptr>& xBuffer, + std::unique_ptr>& bBuffer); + + private: + void participate(bool status) + { + mpiReduce.participate(status); + } + + template + auto reduceGlobal(DataSpace fieldSize, auto dataBoxIn, T_ReduceFunc reduceFunctor = T_ReduceFunc{}); + + std::unique_ptr> pkBuffer; + std::unique_ptr> rkBuffer; + std::unique_ptr> r0Buffer; + std::unique_ptr> mpkBuffer; + std::unique_ptr> ampkBuffer; + std::unique_ptr> zkBuffer; + std::unique_ptr> azkBuffer; + + std::unique_ptr> yBuffer; + std::unique_ptr> wBuffer; + std::unique_ptr> zBuffer; + + std::shared_ptr fieldV; + + std::optional m_mappingDesc; + + mpi::MPIReduce mpiReduce; + std::unique_ptr localReduce; + + // defaults will be overwritten by command line arguments + bool m_useSolver = false; + uint32_t m_maxSolverSteps = 20; + float_64 m_solverEpsilon = 1.0e-8; + + bool m_disablePreconditioner = false; + uint32_t m_maxPreconditionerSteps = 20; + }; +} // namespace picongpu::simulation::stage diff --git a/include/picongpu/simulation/stage/Poisson.x.cpp b/include/picongpu/simulation/stage/Poisson.x.cpp new file mode 100644 index 00000000000..43ee543080a --- /dev/null +++ b/include/picongpu/simulation/stage/Poisson.x.cpp @@ -0,0 +1,680 @@ +/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera + * + * 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 . + */ + +#include "picongpu/simulation/stage/Poisson.hpp" + +#include "picongpu/defines.hpp" +#include "picongpu/fields/FieldJ.hpp" +#include "picongpu/fields/FieldJ.kernel" +#include "picongpu/fields/FieldTmpOperations.hpp" +#include "picongpu/fields/currentDeposition/Deposit.hpp" +#include "picongpu/fields/poissonSolver/BoundaryConditions.hpp" +#include "picongpu/fields/poissonSolver/RightHandSideNormalization.hpp" +#include "picongpu/fields/poissonSolver/Stencil.hpp" +#include "picongpu/particles/filter/filter.hpp" +#include "picongpu/particles/param.hpp" +#include "picongpu/particles/particleToGrid/CombinedDerive.hpp" +#include "picongpu/particles/particleToGrid/ComputeGridValuePerFrame.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +namespace picongpu +{ + namespace simulation + { + namespace stage + { + namespace detail + { + + template + struct ComputeChargeDensity + { + using SpeciesType = pmacc::particles::meta::FindByNameOrType_t; + static uint32_t const area = T_Area::value; + + HINLINE void operator()(FieldTmp& fieldTmp, uint32_t const currentStep) const + { + DataConnector& dc = Environment<>::get().DataConnector(); + + /* load species without copying the particle data to the host */ + auto speciesTmp = dc.get(SpeciesType::FrameType::getName()); + + /* run algorithm */ + using ChargeDensitySolver = typename particles::particleToGrid::CreateFieldTmpOperation_t< + SpeciesType, + particles::particleToGrid::derivedAttributes::ChargeDensity>::Solver; + + computeFieldTmpValue(fieldTmp, *speciesTmp, currentStep); + } + }; + + } // namespace detail + + namespace deriveField = particles::particleToGrid; + template + using SpeciesEligibleForChargeConservation = typename particles::traits:: + SpeciesEligibleForSolver::type; + + Poisson::Poisson() : localReduce{std::make_unique(1024)} + { + } + + void Poisson::registerHelp(po::options_description& desc) + { + namespace po = boost::program_options; + po::options_description solverDesc("Poisson solver"); + + solverDesc.add_options()( + "poisson.activate", + po::value(&m_useSolver)->zero_tokens(), + "enable poisson solver"); + solverDesc.add_options()( + "poisson.maxSteps", + po::value(&m_maxSolverSteps)->default_value(2000), + "maximum number of steps for the solver"); + solverDesc.add_options()( + "poisson.epsilon", + po::value(&m_solverEpsilon)->default_value(1.0e-8), + "maximal allowed error of the poisson solver"); + // preconitioner + solverDesc.add_options()( + "poisson.preconditioner.disable", + po::value(&m_disablePreconditioner)->zero_tokens(), + "disable poisson solver preconditioner"); + solverDesc.add_options()( + "poisson.preconditioner.maxSteps", + po::value(&m_maxPreconditionerSteps)->default_value(20), + "maximum number of steps for the preconditioner"); + desc.add(solverDesc); + } + + void Poisson::init(MappingDesc const mappingDesc) + { + m_mappingDesc = std::make_optional(mappingDesc); + + if(m_useSolver) + { + pkBuffer = std::make_unique>(m_mappingDesc->getGridLayout()); + rkBuffer = std::make_unique>(m_mappingDesc->getGridLayout()); + r0Buffer = std::make_unique>(m_mappingDesc->getGridLayout()); + mpkBuffer = std::make_unique>(m_mappingDesc->getGridLayout()); + ampkBuffer = std::make_unique>(m_mappingDesc->getGridLayout()); + zkBuffer = std::make_unique>(m_mappingDesc->getGridLayout()); + azkBuffer = std::make_unique>(m_mappingDesc->getGridLayout()); + fieldV = std::make_shared(m_mappingDesc.value()); + + yBuffer = std::make_unique>(m_mappingDesc->getGridLayout()); + wBuffer = std::make_unique>(m_mappingDesc->getGridLayout()); + zBuffer = std::make_unique>(m_mappingDesc->getGridLayout()); + + auto const commTag0 = pmacc::traits::getUniqueId(); + auto const commTag1 = pmacc::traits::getUniqueId(); + auto const commTagPk = pmacc::traits::getUniqueId(); + auto const commTagRk = pmacc::traits::getUniqueId(); + auto const commTagFieldV = pmacc::traits::getUniqueId(); + + auto const commTagY = pmacc::traits::getUniqueId(); + /*go over all directions*/ + for(uint32_t i = 1; i < NumberOfExchanges::value; ++i) + { + // set communication only for planes + if(FRONT % i == 0) + { + DataSpace relativeMask = Mask::getRelativeDirections(i); + /* guarding cells depend on direction + * for negative direction use originGuard else endGuard (relative direction ZERO is + * ignored) don't switch end and origin because this is a read buffer and no send buffer + */ + auto guardingCells = DataSpace::create(0); + for(uint32_t d = 0; d < simDim; ++d) + guardingCells[d] = (relativeMask[d] == 0 ? 0 : 1); + mpkBuffer->addExchange(GUARD, i, guardingCells, commTag0); + zkBuffer->addExchange(GUARD, i, guardingCells, commTag1); + pkBuffer->addExchange(GUARD, i, guardingCells, commTagPk); + rkBuffer->addExchange(GUARD, i, guardingCells, commTagRk); + fieldV->fieldVBuffer->addExchange(GUARD, i, guardingCells, commTagFieldV); + + yBuffer->addExchange(GUARD, i, guardingCells, commTagY); + } + } + DataConnector& dc = Environment<>::get().DataConnector(); + dc.share(fieldV); + + participate(true); + } + } + + template + class TransformDataBox : private T_TranformFunctor + { + public: + using ValueType = decltype(std::declval()(DataSpace::create(0))); + + static constexpr std::uint32_t Dim = simDim; + + HDINLINE TransformDataBox() = default; + + HDINLINE TransformDataBox(T_TranformFunctor transformFunc) : T_TranformFunctor(transformFunc) + { + } + + HDINLINE TransformDataBox(TransformDataBox const&) = default; + + HDINLINE ValueType operator()(DataSpace const& idx) const + { + return T_TranformFunctor::operator()(idx + m_offset); + } + + HDINLINE ValueType operator[](DataSpace const idx) const + { + return T_TranformFunctor::operator()(idx + m_offset); + } + + HDINLINE TransformDataBox shift(DataSpace const& offset) const + { + TransformDataBox result(*this); + result.m_offset += offset; + return result; + } + + DataSpace m_offset = DataSpace::create(0); + }; + + template + inline auto Poisson::reduceGlobal( + DataSpace fieldSize, + auto dataBoxIn, + T_TranformFunctor reduceFunctor) + { + DataBoxDim1Access d1Access(dataBoxIn, fieldSize); + + float_64 resultLocal = (*localReduce)(reduceFunctor, d1Access, fieldSize.productOfComponents()); + + // avoid deadlock between not finished pmacc tasks and mpi blocking collectives + eventSystem::getTransactionEvent().waitForFinished(); + float_64 resultGlobal; + mpiReduce(reduceFunctor, &resultGlobal, &resultLocal, 1, mpi::reduceMethods::AllReduce()); + + return resultGlobal; + } + + void Poisson::preconditioner( + std::unique_ptr>& xBuffer, + std::unique_ptr>& bBuffer) + { + yBuffer->getDeviceBuffer().setValue(0.0); + wBuffer->getDeviceBuffer().setValue(0.0); + zBuffer->getDeviceBuffer().setValue(0.0); + + SubGrid const& subGrid = Environment::get().SubGrid(); + auto globalDomain = subGrid.getGlobalDomain().size; + auto cellSizeSquared = sim.pic.getCellSize() * sim.pic.getCellSize(); + + float_64 eigenMin = 0.0; + float_64 eigenMax = 0.0; + + for(uint32_t d = 0; d < simDim; ++d) + { + eigenMin += 4.0 * math::sin(1.0 * pmacc::math::Pi::halfValue / (globalDomain[d] + 1)) + * math::sin(1.0 * pmacc::math::Pi::halfValue / (globalDomain[d] + 1)) + / (cellSizeSquared[d]); + + eigenMax + += 4.0 + * math::sin(globalDomain[d] * pmacc::math::Pi::halfValue / (globalDomain[d] + 1)) + * math::sin(globalDomain[d] * pmacc::math::Pi::halfValue / (globalDomain[d] + 1)) + / (cellSizeSquared[d]); + } + + float_64 const theta = 0.5 * (eigenMax + eigenMin); + float_64 const delta = 0.5 * (eigenMax - eigenMin); + float_64 const sigma = theta / delta; + + float_64 rhoOld = 1. / sigma; + float_64 rhoCurrent = 1. / (2. * sigma - rhoOld); + + bBuffer->communication(); + + auto coreBorderMapper = makeAreaMapper(m_mappingDesc.value()); + + namespace poi = fields::poissonSolver; + + pmacc::algorithms::forEachCell( + coreBorderMapper, + [theta] DEVICEONLY(auto, auto const bValue) -> float_64 { return bValue / theta; }, + zBuffer->getDeviceBuffer(), + bBuffer->getDeviceBuffer()); + + poi::stencil( + coreBorderMapper, + poi::StencilFunc{}, + yBuffer->getDeviceBuffer(), + bBuffer->getDeviceBuffer()); + + pmacc::algorithms::forEachCell( + coreBorderMapper, + [rhoCurrent, theta, delta] DEVICEONLY(auto const yValue) -> float_64 + { return -2.0 * rhoCurrent * yValue / theta / delta; }, + yBuffer->getDeviceBuffer()); + + pmacc::algorithms::forEachCell( + coreBorderMapper, + [rhoCurrent, delta] DEVICEONLY(auto const yValue, auto const bValue) -> float_64 + { return yValue + 4.0 * bValue * rhoCurrent / delta; }, + yBuffer->getDeviceBuffer(), + bBuffer->getDeviceBuffer()); + + + uint32_t iterMax = m_maxPreconditionerSteps; + for(uint32_t i = 2; i < iterMax; ++i) + { + rhoOld = rhoCurrent; + rhoCurrent = 1. / (2. * sigma - rhoOld); + yBuffer->communication(); + + poi::stencil( + coreBorderMapper, + poi::StencilFunc{}, + wBuffer->getDeviceBuffer(), + yBuffer->getDeviceBuffer()); + + pmacc::algorithms::forEachCell( + coreBorderMapper, + [rhoCurrent, delta] DEVICEONLY(auto const wValue) -> float_64 + { return -2.0 * rhoCurrent * wValue / delta; }, + wBuffer->getDeviceBuffer()); + + pmacc::algorithms::forEachCell( + coreBorderMapper, + [rhoCurrent, delta] DEVICEONLY(auto const wValue, auto const bValue) -> float_64 + { return wValue + 2.0 * rhoCurrent * bValue / delta; }, + wBuffer->getDeviceBuffer(), + bBuffer->getDeviceBuffer()); + + pmacc::algorithms::forEachCell( + coreBorderMapper, + [rhoCurrent, sigma] DEVICEONLY(auto const wValue, auto const yValue) -> float_64 + { return wValue + 2.0 * rhoCurrent * sigma * yValue; }, + wBuffer->getDeviceBuffer(), + yBuffer->getDeviceBuffer()); + + pmacc::algorithms::forEachCell( + coreBorderMapper, + [rhoCurrent, rhoOld] DEVICEONLY(auto const wValue, auto const zValue) -> float_64 + { return wValue - rhoCurrent * rhoOld * zValue; }, + wBuffer->getDeviceBuffer(), + zBuffer->getDeviceBuffer()); + + zBuffer->getDeviceBuffer().copyFrom(yBuffer->getDeviceBuffer()); + yBuffer->getDeviceBuffer().copyFrom(wBuffer->getDeviceBuffer()); + } // loop + xBuffer->getDeviceBuffer().copyFrom(wBuffer->getDeviceBuffer()); + } + + void Poisson::operator()(uint32_t const currentStep) + { + namespace poi = fields::poissonSolver; + + // only one rank is writing status initial information of the poisson solver + bool mainRank = Environment::get().GridController().getGlobalRank() == 0; + if(mainRank) + log("Poisson solver:"); + if(!m_useSolver) + { + if(mainRank) + log(" - disabled"); + return; + } + eventSystem::getTransactionEvent().waitForFinished(); + auto beginT = std::chrono::high_resolution_clock::now(); + + using namespace pmacc; + constexpr uint fieldRhoSlot = 0; + DataConnector& dc = Environment<>::get().DataConnector(); + auto& fieldRho = *dc.get(FieldTmp::getUniqueId(fieldRhoSlot)); + + DataSpace numGuardCells = fieldRho.getGridLayout().guardSizeND(); + DataSpace coreBorderSize = fieldRho.getGridLayout().sizeWithoutGuardND(); + + using EligibleSpecies = pmacc::mp_filter; + /* calculate and add the charge density values from all species in FieldTmp */ + meta::ForEach< + EligibleSpecies, + detail::ComputeChargeDensity>, + boost::mpl::_1> + computeChargeDensity; + + fieldRho.getGridBuffer().getDeviceBuffer().setValue(FieldTmp::ValueType(0.0)); + computeChargeDensity(fieldRho, currentStep); + + /* add results of all species that are still in GUARD to next GPUs BORDER */ + EventTask fieldTmpEvent = fieldRho.asyncCommunication(eventSystem::getTransactionEvent()); + eventSystem::setTransactionEvent(fieldTmpEvent); + + auto boundaryConditionsDirichlet = poi::BoundaryConditionsDirichlet{}; + boundaryConditionsDirichlet(*fieldV.get(), m_mappingDesc.value()); + + auto rightHandSideNormalization = poi::RightHandSideNormalization{}; + rightHandSideNormalization(*fieldV.get(), fieldRho, m_mappingDesc.value()); + + + float_64 normRho; + { + auto rhoDeviceBox = fieldRho.getDeviceDataBox().shift(numGuardCells); + TransformDataBox fieldTransform( + [rhoDeviceBox] DEVICEONLY(DataSpace const& idx) -> float_64 + { return precisionCast(rhoDeviceBox[idx].x() * rhoDeviceBox[idx].x()); }); + + normRho = std::sqrt(reduceGlobal(coreBorderSize, fieldTransform)); + } + // recalculate rho + fieldRho.getGridBuffer().getDeviceBuffer().setValue(FieldTmp::ValueType(0.0)); + computeChargeDensity(fieldRho, currentStep); + /* add results of all species that are still in GUARD to next GPUs BORDER */ + eventSystem::setTransactionEvent(fieldRho.asyncCommunication(eventSystem::getTransactionEvent())); + + auto coreBorderMapper = makeAreaMapper(m_mappingDesc.value()); + + // normalize rho + pmacc::algorithms::forEachCell( + coreBorderMapper, + [normRho] DEVICEONLY(auto rhoValue) -> float_64 { return rhoValue.x() / normRho; }, + fieldRho.getGridBuffer().getDeviceBuffer()); + + auto vMapper = makeAreaMapper(m_mappingDesc.value()); + + // normalize v + pmacc::algorithms::forEachCell( + vMapper, + [normRho] DEVICEONLY(auto const vValue) -> float_64 { return vValue / normRho; }, + fieldV->fieldVBuffer->getDeviceBuffer()); + + poi::stencil( + coreBorderMapper, + poi::StencilFunc{}, + r0Buffer->getDeviceBuffer(), + fieldV->fieldVBuffer->getDeviceBuffer()); + + pmacc::algorithms::forEachCell( + coreBorderMapper, + [] DEVICEONLY(auto const r0Value, auto const rhoValue) -> float_64 + { return rhoValue.x() - r0Value; }, + r0Buffer->getDeviceBuffer(), + fieldRho.getGridBuffer().getDeviceBuffer()); + + pkBuffer->getDeviceBuffer().copyFrom(r0Buffer->getDeviceBuffer()); + rkBuffer->getDeviceBuffer().copyFrom(r0Buffer->getDeviceBuffer()); + + float_64 rho0; + /* rho reduction */ + { + auto r0Box = r0Buffer->getDeviceBuffer().getDataBox(); + auto r0BoxBorderGuard = r0Box.shift(numGuardCells); + + TransformDataBox fieldTransform( + [r0BoxBorderGuard] DEVICEONLY(DataSpace const& idx) -> float_64 + { return r0BoxBorderGuard[idx] * r0BoxBorderGuard[idx]; }); + + rho0 = reduceGlobal(coreBorderSize, fieldTransform); + } + + float_64 rho1 = rho0; + float_64 totalSum2; + + int maxIterations = m_maxSolverSteps; + + bool foundSolution = false; + int iteration = 0; + for(; iteration < maxIterations; ++iteration) + { + // preconditioner + if(m_disablePreconditioner) + mpkBuffer->getDeviceBuffer().copyFrom(pkBuffer->getDeviceBuffer()); + else + preconditioner(mpkBuffer, pkBuffer); + + mpkBuffer->communication(); + + // w = Ap + poi::stencil( + coreBorderMapper, + poi::StencilFunc{}, + ampkBuffer->getDeviceBuffer(), + mpkBuffer->getDeviceBuffer()); + + float_64 totalSum1; + /* local p = rw */ + { + auto r0Box = r0Buffer->getDeviceBuffer().getDataBox(); + auto r0BoxBorderGuard = r0Box.shift(numGuardCells); + + auto ampkBox = ampkBuffer->getDeviceBuffer().getDataBox(); + auto ampkBoxBorderGuard = ampkBox.shift(numGuardCells); + + TransformDataBox fieldTransform( + [r0BoxBorderGuard, ampkBoxBorderGuard] DEVICEONLY(DataSpace const& idx) -> float_64 + { return r0BoxBorderGuard[idx] * ampkBoxBorderGuard[idx]; }); + + totalSum1 = reduceGlobal(coreBorderSize, fieldTransform); + } + float_64 alpha = rho0 / totalSum1; + // r = r - alpha * w + pmacc::algorithms::forEachCell( + coreBorderMapper, + [alpha] DEVICEONLY(auto const rkValue, auto const ampkValue) -> float_64 + { return rkValue - alpha * ampkValue; }, + rkBuffer->getDeviceBuffer(), + ampkBuffer->getDeviceBuffer()); + + // preconditioner + if(m_disablePreconditioner) + zkBuffer->getDeviceBuffer().copyFrom(rkBuffer->getDeviceBuffer()); + else + preconditioner(zkBuffer, rkBuffer); + + zkBuffer->communication(); + + // t = A * r + poi::stencil( + coreBorderMapper, + poi::StencilFunc{}, + azkBuffer->getDeviceBuffer(), + zkBuffer->getDeviceBuffer()); + + /* totalSum1 = azk * rk */ + { + auto azkBox = azkBuffer->getDeviceBuffer().getDataBox(); + auto azkBoxBorderGuard = azkBox.shift(numGuardCells); + + auto rkBox = rkBuffer->getDeviceBuffer().getDataBox(); + auto rkBoxBorderGuard = rkBox.shift(numGuardCells); + + TransformDataBox fieldTransform( + [azkBoxBorderGuard, rkBoxBorderGuard] DEVICEONLY(DataSpace const& idx) -> float_64 + { return azkBoxBorderGuard[idx] * rkBoxBorderGuard[idx]; }); + + totalSum1 = reduceGlobal(coreBorderSize, fieldTransform); + } + + /* totalSum1 = azk * azk */ + { + auto azkBox = azkBuffer->getDeviceBuffer().getDataBox(); + auto azkBoxBorderGuard = azkBox.shift(numGuardCells); + + TransformDataBox fieldTransform( + [azkBoxBorderGuard] DEVICEONLY(DataSpace const& idx) -> float_64 + { return azkBoxBorderGuard[idx] * azkBoxBorderGuard[idx]; }); + + totalSum2 = reduceGlobal(coreBorderSize, fieldTransform); + } + + float_64 omega = totalSum1 / totalSum2; + + // v = v + alpha * mpk + omega * zk + pmacc::algorithms::forEachCell( + coreBorderMapper, + [alpha, + omega] DEVICEONLY(auto const vValue, auto const mpkValue, auto const zkValue) -> float_64 + { return vValue + alpha * mpkValue + omega * zkValue; }, + fieldV->fieldVBuffer->getDeviceBuffer(), + mpkBuffer->getDeviceBuffer(), + zkBuffer->getDeviceBuffer()); + + // rk = rk - omega * azk + pmacc::algorithms::forEachCell( + coreBorderMapper, + [omega] DEVICEONLY(auto const rkValue, auto const azkValue) -> float_64 + { return rkValue - omega * azkValue; }, + rkBuffer->getDeviceBuffer(), + azkBuffer->getDeviceBuffer()); + + /* totalSum1 = r0 * rk */ + { + auto r0Box = r0Buffer->getDeviceBuffer().getDataBox(); + auto r0BoxBorderGuard = r0Box.shift(numGuardCells); + + auto rkBox = rkBuffer->getDeviceBuffer().getDataBox(); + auto rkBoxBorderGuard = rkBox.shift(numGuardCells); + + TransformDataBox fieldTransform( + [r0BoxBorderGuard, rkBoxBorderGuard] DEVICEONLY(DataSpace const& idx) -> float_64 + { return r0BoxBorderGuard[idx] * rkBoxBorderGuard[idx]; }); + + totalSum1 = reduceGlobal(coreBorderSize, fieldTransform); + } + /* totalSum2 = rk * rk */ + { + auto rkBox = rkBuffer->getDeviceBuffer().getDataBox(); + auto rkBoxBorderGuard = rkBox.shift(numGuardCells); + + TransformDataBox fieldTransform( + [rkBoxBorderGuard] DEVICEONLY(DataSpace const& idx) -> float_64 + { return rkBoxBorderGuard[idx] * rkBoxBorderGuard[idx]; }); + + totalSum2 = reduceGlobal(coreBorderSize, fieldTransform); + } + + rho1 = totalSum1; + float_64 beta = rho1 / rho0 * alpha / omega; + rho0 = rho1; + if(std::sqrt(totalSum2) < m_solverEpsilon) + { + foundSolution = true; + break; + } + // pk = rk + beta * (pk - omega * ampk) + pmacc::algorithms::forEachCell( + coreBorderMapper, + [beta, + omega] DEVICEONLY(auto const pkValue, auto const rkValue, auto const ampkValue) -> float_64 + { return rkValue + beta * (pkValue - omega * ampkValue); }, + pkBuffer->getDeviceBuffer(), + rkBuffer->getDeviceBuffer(), + ampkBuffer->getDeviceBuffer()); + + } // for loop + + // avid deadlock due blocking colective MPI operation + eventSystem::getTransactionEvent().waitForFinished(); + + bool ioRank = mpiReduce.hasResult(mpi::reduceMethods::Reduce()); + + int maxGlobalIterations = 0; + mpiReduce( + pmacc::math::operation::Max(), + &maxGlobalIterations, + &iteration, + 1, + mpi::reduceMethods::Reduce()); + + float_64 maxGlobalNormRho = 0.0; + mpiReduce(pmacc::math::operation::Max(), &maxGlobalNormRho, &normRho, 1, mpi::reduceMethods::Reduce()); + + float_64 maxGlobalTotalSum2 = 0.0; + mpiReduce( + pmacc::math::operation::Max(), + &maxGlobalTotalSum2, + &totalSum2, + 1, + mpi::reduceMethods::Reduce()); + + if(foundSolution) + { + if(ioRank) + log(" - converged after %1%/%2% iterations with norm=%3%, total epsilon=%4%") + % maxGlobalIterations % maxIterations % maxGlobalNormRho % std::sqrt(maxGlobalTotalSum2); + + // normalize v back + pmacc::algorithms::forEachCell( + coreBorderMapper, + [normRho] DEVICEONLY(auto const vValue) -> float_64 + { return vValue * normRho / sim.pic.getEps0(); }, + fieldV->fieldVBuffer->getDeviceBuffer()); + fieldV->fieldVBuffer->communication(); + + // compute fieldE + auto& eField = *dc.get(FieldE::getName()); + + poi::stencil( + coreBorderMapper, + poi::GetFieldEStencil{}, + eField.getGridBuffer().getDeviceBuffer(), + fieldV->fieldVBuffer->getDeviceBuffer()); + setTransactionEvent(eField.asyncCommunication(eventSystem::getTransactionEvent())); + } + + eventSystem::getTransactionEvent().waitForFinished(); + auto endT = std::chrono::high_resolution_clock::now(); + double duration = std::chrono::duration(endT - beginT).count(); + + // avid deadlock due blocking collective MPI operation + eventSystem::getTransactionEvent().waitForFinished(); + double globalMaxDuration = 0.0; + mpiReduce( + pmacc::math::operation::Max(), + &globalMaxDuration, + &duration, + 1, + mpi::reduceMethods::Reduce()); + if(ioRank) + log(" - duration %1% sec") % globalMaxDuration; + + if(ioRank && !foundSolution) + { + log(" - not converge after %1% iterations with norm=%2%, total epsilon=%3%") + % maxIterations % normRho % std::sqrt(rho1); + throw std::runtime_error("Poisson solver did not converge after max iterations"); + } + } + } // namespace stage + } // namespace simulation +} // namespace picongpu diff --git a/include/pmacc/algorithms/ForEachCell.hpp b/include/pmacc/algorithms/ForEachCell.hpp new file mode 100644 index 00000000000..f88030d96b5 --- /dev/null +++ b/include/pmacc/algorithms/ForEachCell.hpp @@ -0,0 +1,108 @@ +/* Copyright 2025-2025 , Rene Widera, Edgar Marquardt + * + * This file is part of PMacc. + * + * PMacc is free software: you can redistribute it and/or modify + * it under the terms of either the GNU General Public License or + * the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * PMacc 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 and the GNU Lesser General Public License + * for more details. + * + * You should have received a copy of the GNU General Public License + * and the GNU Lesser General Public License along with PMacc. + * If not, see . + */ + +#pragma once + +#include "pmacc/Environment.hpp" +#include "pmacc/attribute/FunctionSpecifier.hpp" +#include "pmacc/dimensions/DataSpace.hpp" +#include "pmacc/lockstep/ForEach.hpp" +#include "pmacc/lockstep/Kernel.hpp" +#include "pmacc/types.hpp" + +#include +#include + +#include "pmacc/math/vector/compile-time/Vector.hpp" + +namespace pmacc::algorithms +{ + template + struct DeviceLambda + { + T_Func const func; + + template + DEVICEONLY auto operator()(T&&... args) const + { + return func(std::forward(args)...); + } + + template + DEVICEONLY auto operator()(T&&... args) + { + return func(std::forward(args)...); + } + }; + + template + DeviceLambda(T_Func const) -> DeviceLambda; +} // namespace pmacc::algorithms + +namespace alpaka +{ + template + struct IsKernelArgumentTriviallyCopyable, void> : std::true_type + { + }; +} // namespace alpaka + +namespace pmacc +{ + namespace algorithms + { + template + struct ForEachCellKernel + { + TMapper mapper; + + DINLINE auto operator()(auto const& worker, auto const& func, auto outBox, auto... boxes) const -> void + { + constexpr auto simDim = TMapper::Dim; + using SuperCellSize = typename TMapper::SuperCellSize; + DataSpace const superCellIdx(mapper.getSuperCellIndex(worker.blockDomIdxND())); + DataSpace superCellCellOffset = superCellIdx * SuperCellSize::toRT(); + + constexpr uint32_t cellsPerSuperCell = pmacc::math::CT::volume::type::value; + + auto forEachCellInSupercell = lockstep::makeForEach(worker); + + forEachCellInSupercell( + [&](int32_t const linearCellIdx) + { + DataSpace const cellIdx = pmacc::math::mapToND(SuperCellSize::toRT(), linearCellIdx); + DataSpace const dataCellOffset = superCellCellOffset + cellIdx; + outBox[dataCellOffset] = func(outBox[dataCellOffset], boxes[dataCellOffset]...); + }); + } + }; + + template + inline void forEachCell(TMapper mapper, auto const& functor, auto& outBuffer, auto&... buffers) + { + PMACC_LOCKSTEP_KERNEL(ForEachCellKernel{mapper}) + .config(mapper.getGridDim(), typename TMapper::SuperCellSize{})( + DeviceLambda{functor}, + outBuffer.getDataBox(), + buffers.getDataBox()...); + } + } // namespace algorithms +} // namespace pmacc diff --git a/share/picongpu/tests/PoissonSolver/README.rst b/share/picongpu/tests/PoissonSolver/README.rst new file mode 100644 index 00000000000..66135e56d21 --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/README.rst @@ -0,0 +1,10 @@ +=========================== +Test for the Poisson solver +=========================== + +This test is for validating the PoissonSolver in the PIConGPU initialization. + +To run this test, one has to execute ci.sh with the location of the input and output directory. + +..code-block:: bash +./picongpu/share/picongpu/tests/PoissonSolver/bin/ci.sh picongpu/share/picongpu/tests/PoissonSolver/ ./run02 diff --git a/share/picongpu/tests/PoissonSolver/bin/ci.sh b/share/picongpu/tests/PoissonSolver/bin/ci.sh new file mode 100755 index 00000000000..1bb7a63cbdb --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/bin/ci.sh @@ -0,0 +1,157 @@ +#!/bin/bash +# +# This file is part of PIConGPU. +# +# Copyright 2026-2026 PIConGPU contributors +# Author: Edgar Marquardt +# +# 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 . +# + +function absolute_path() +{ + builtin cd -- $1 && pwd +} + +help() +{ + echo "Tests the Poisson solver with a simple setup." + echo "" + echo "Usage: ci.sh [-d dataPath] [inputSetPath] [destinationPath]" + echo "" + echo "Options" + echo "-h | --help - show help" + echo "" + echo "inputSetPath - path to the simulation input set" + echo " Default: current directory" + echo "destinationPath - path to the destination where the input set is cloned to via" + echo " 'pic-create'" + echo " Default: current directory" +} + +##################### +## option handling ## +##################### +# options may be followed by +# - one colon to indicate they have a required argument +OPTS=`getopt -o h -l help -- "$@"` +if [ $? != 0 ] ; then + # something went wrong, getopt will put out an error message for us + exit 1 +fi + +eval set -- "$OPTS" + +# parser +while true ; do + case "$1" in + -h|--help) + echo -e "$(help)" + shift + exit 0 + ;; + --) shift; break;; + esac + shift +done + + +############################ +## build and run picongpu ## +############################ +if [ $# -eq 2 ] ; then + inputSetPath=$1 + inputDestinationPath=$2 +else + echo "Two arguments are required, $# given!" >&2 + echo -e "$(help)" >&2 + exit 1 +fi + +if [ -d "$inputSetPath/include" ] ; then + if [ -d "$inputDestinationPath" ] ; then + echo "Output directory $inputDestinationPath exists" >&2 + echo "removing now" >&2 + rm -r $inputDestinationPath + fi + echo "start setting up" + pic-create $inputSetPath $inputDestinationPath + ret_create=$? + if [ $ret_create -ne 0 ] ; then + echo "pic-create failed" + exit $ret_create + fi + + inputDestinationPath=$(absolute_path $inputDestinationPath) + cd $inputDestinationPath + + echo "building" + pic-build + ret_build=$? + +else + echo "Input path $inputSetPath does not contain an include directory" >&2 + exit 2 +fi + +if [ $ret_build -eq 0 ] ; then + ## create simulation data directory + simPath="./simOutput" + + if [ -d "$simPath" ] ; then + echo "Simulation path already in use, cannot create new folder" >&2 + exit 3 + fi + + mkdir -p $simPath + + # use absolut path's + simPath=$(absolute_path $simPath) + + cd $simPath + + # run the simulation + echo "Simulation path: " $simPath"/" + mpiexec -n 1 ../bin/picongpu -d 1 1 1 -g 64 64 64 --periodic 1 1 1 -s 4 \ + --openPMD.period 8 --openPMD.ext bp --openPMD.file simData --poisson.activate + + ret_sim=$? + if [ $ret_sim -ne 0 ] ; then + echo "running simulation failed" >&2 + exit $ret_sim + fi + + cd .. + # validate the results + if [ -d "./lib/python/test/" ] ; then + python3 ./lib/python/test/validate_results.py -r $simPath + + ret=$? + if [ $ret -eq 0 ] ; then + echo "test successfully validated" + else + echo "test validation failed" + fi + exit $ret + + else + echo "Input path $inputSetPath does not contain an lib/python/test directory" >&2 + exit 2 + fi + +else + echo "build failed" >&2 + exit $ret_build +fi diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/density.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/density.param new file mode 100644 index 00000000000..f354524edd6 --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/density.param @@ -0,0 +1,32 @@ +/* Copyright 2013-2024 Axel Huebl, Heiko Burau, Rene Widera, Felix Schmitt, + * Richard Pausch + * + * 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 . + */ + +#pragma once + +#include "picongpu/particles/densityProfiles/profiles.def" + +namespace picongpu +{ + namespace densityProfiles + { + /* definition of homogenous profile */ + using Homogenous = HomogenousImpl; + } // namespace densityProfiles +} // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/dimension.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/dimension.param new file mode 100644 index 00000000000..7ba41c27a77 --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/dimension.param @@ -0,0 +1,31 @@ +/* Copyright 2014-2024 Axel Huebl, Rene Widera + * + * 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 . + */ + +#pragma once + +#ifndef PARAM_DIMENSION +# define PARAM_DIMENSION DIM3 +#endif + +#define SIMDIM PARAM_DIMENSION + +namespace picongpu +{ + constexpr uint32_t simDim = SIMDIM; +} // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/particle.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/particle.param new file mode 100644 index 00000000000..e11790099ad --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/particle.param @@ -0,0 +1,70 @@ +/* Copyright 2013-2024 Axel Huebl, Rene Widera, Benjamin Worpitz, + * Richard Pausch, Klaus Steiniger + * + * 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 . + */ + +/** @file + * + * Configurations for particle manipulators. Set up and declare functors that + * can be used in speciesInitialization.param for particle species + * initialization and manipulation, such as temperature distributions, drifts, + * pre-ionization and in-cell position. + */ + +#pragma once + +#include "picongpu/particles/manipulators/manipulators.def" +#include "picongpu/particles/startPosition/functors.def" + +#include + +namespace picongpu +{ + namespace particles + { + /** a particle with a weighting below MIN_WEIGHTING will not + * be created / will be deleted + * + * unit: none */ + constexpr float_X MIN_WEIGHTING = 10.0; + + namespace startPosition + { + /** Define target number for marco-particles per cell + * to be used in Random start position functor. + */ + struct RandomParameter + { + /** Maximum number of macro-particles per cell during density profile evaluation. + * + * Determines the weighting of a macro particle as well as the number of + * macro-particles which sample the evolution of the particle distribution + * function in phase space. + * + * unit: none + */ + static constexpr uint32_t numParticlesPerCell = 2u; + }; + + /** Definition of start position functor that randomly distributes macro-particles within a cell. */ + + using Random = RandomImpl; + + } // namespace startPosition + } // namespace particles +} // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/simulation.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/simulation.param new file mode 100644 index 00000000000..e54091dbb9a --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/simulation.param @@ -0,0 +1,84 @@ +/* Copyright 2013-2024 Axel Huebl, Rene Widera, Benjamin Worpitz + * + * 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 . + */ + +/** @file + * + * Definition of cell sizes and time step. Our cells are defining a regular, + * cartesian grid. Our explicit FDTD field solvers require an upper bound for + * the time step value in relation to the cell size for convergence. Make + * sure to resolve important wavelengths of your simulation, e.g. shortest + * plasma wavelength, Debye length and central laser wavelength both spatially + * and temporarily. + * + * **Units in reduced dimensions** + * + * In 2D3V simulations, the CELL_DEPTH_SI (Z) cell length + * is still used for normalization of densities, etc.. + * + * A 2D3V simulation in a cartesian PIC simulation such as + * ours only changes the degrees of freedom in motion for + * (macro) particles and all (field) information in z + * travels instantaneously, making the 2D3V simulation + * behave like the interaction of infinite "wire particles" + * in fields with perfect symmetry in Z. + * + */ + +#pragma once + +namespace picongpu +{ + namespace SI + { + /** Duration of one timestep + * unit: seconds */ + /** Duration of one timestep + * unit: seconds */ + constexpr float_64 DELTA_T_SI = 2.0564e-16; + + /** equals X + * unit: meter */ + constexpr float_64 CELL_WIDTH_SI = 2.0e-6; + /** equals Y - the laser & moving window propagation direction + * unit: meter */ + constexpr float_64 CELL_HEIGHT_SI = 6.202e-08; + /** equals Z + * unit: meter */ + constexpr float_64 CELL_DEPTH_SI = CELL_WIDTH_SI; + + /** Base density in particles per m^3 in the density profiles. + * + * This is often taken as reference maximum density in normalized profiles. + * Individual particle species can define a `densityRatio` flag relative + * to this value. + * + * unit: ELEMENTS/m^3 + */ +#ifndef PARAM_BASE_DENSITY_SI +# define PARAM_BASE_DENSITY_SI 1.e25 +#endif + constexpr float_64 BASE_DENSITY_SI = PARAM_BASE_DENSITY_SI; + } // namespace SI + + /** Approximate number of maximum macro-particles per cell. + * + * Used internally for unit normalization. + */ + constexpr uint32_t TYPICAL_PARTICLES_PER_CELL = 2u; +} // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesDefinition.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesDefinition.param new file mode 100644 index 00000000000..72b3bdf7b44 --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesDefinition.param @@ -0,0 +1,91 @@ +/* Copyright 2013-2024 Rene Widera, Marco Garten, Richard Pausch, + * Benjamin Worpitz, Axel Huebl + * + * 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 . + */ + +#pragma once + + +#include "picongpu/particles/Particles.hpp" + +#include +#include +#include +#include +#include + +namespace picongpu +{ + /*########################### define particle attributes #####################*/ + + /** describe attributes of a particle*/ + using DefaultParticleAttributes = MakeSeq_t, momentum, weighting>; + + /*########################### end particle attributes ########################*/ + + /*########################### define species #################################*/ + + /*--------------------------- electrons --------------------------------------*/ + + /* ratio relative to base charge and base mass */ + value_identifier(float_X, MassRatioElectrons, 1.0); + value_identifier(float_X, ChargeRatioElectrons, 1.0); + + using ParticleFlagsElectrons = MakeSeq_t< + particlePusher, + shape, + interpolation, + current, + massRatio, + chargeRatio>; + + /* define species electrons */ + using PIC_Electrons = Particles; + + /*--------------------------- ions -------------------------------------------*/ + + /* ratio relative to base charge and base mass */ + value_identifier(float_X, MassRatioIons, 1836.152672); + value_identifier(float_X, ChargeRatioIons, -1.0); + + /* ratio relative to BASE_DENSITY */ + value_identifier(float_X, DensityRatioIons, 1.0); + + using ParticleFlagsIons = MakeSeq_t< + particlePusher, + shape, + interpolation, + current, + massRatio, + chargeRatio, + densityRatio, + atomicNumbers>; + + /* define species ions */ + using PIC_Ions = Particles; + + /*########################### end species ####################################*/ + + /** All known particle species of the simulation + * + * List all defined particle species from above in this list + * to make them available to the PIC algorithm. + */ + using VectorAllSpecies = MakeSeq_t; + +} // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesInitialization.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesInitialization.param new file mode 100644 index 00000000000..610c2b56c5b --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesInitialization.param @@ -0,0 +1,47 @@ +/* Copyright 2015-2024 Rene Widera, Axel Huebl + * + * 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 . + */ + +/** @file + * + * Initialize particles inside particle species. This is the final step in + * setting up particles (defined in `speciesDefinition.param`) via density + * profiles (defined in `density.param`). One can then further derive particles + * from one species to another and manipulate attributes with "manipulators" + * and "filters" (defined in `particle.param` and `particleFilters.param`). + */ + +#pragma once + +#include "picongpu/particles/InitFunctors.hpp" + +namespace picongpu +{ + namespace particles + { + /** InitPipeline define in which order species are initialized + * + * the functors are called in order (from first to last functor) + */ + using InitPipeline = pmacc::mp_list< + + CreateDensity, + CreateDensity>; + + } // namespace particles +} // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/lib/python/test/validate_results.py b/share/picongpu/tests/PoissonSolver/lib/python/test/validate_results.py new file mode 100755 index 00000000000..a181cac6952 --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/lib/python/test/validate_results.py @@ -0,0 +1,117 @@ +#!/usr/bin/env python3 +# +# Copyright 2026-2026 Edgar Marquardt +# +# 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 . +# + +""" +Test for the Poisson solver. It checks if the divergence of the electric field +matches the charge density. +""" + +import argparse +import sys + +import openpmd_api as io +from scipy.constants import epsilon_0 +import numpy as np + +parser = argparse.ArgumentParser(description="1") + + +parser.add_argument( + "-r", + help="Path to the simulation results", + dest="path", + type=str, +) + +args = parser.parse_args() + +# Get the simulation data +filename = args.path + "/openPMD/simData_000000.bp" + +series = io.Series(filename, io.Access.read_only) +m = series.iterations[0].meshes + +# collect the total charge density +rho_e = m["e_all_chargeDensity"][:, :, :] +rho_i = m["i_all_chargeDensity"][:, :, :] +series.flush() +rho_e *= m["e_all_chargeDensity"].get_attribute("unitSI") +rho_i *= m["i_all_chargeDensity"].get_attribute("unitSI") + +rho = rho_e + rho_i +rho = np.transpose(rho) + +# get the grid parameters +spacing = np.array(m["e_all_chargeDensity"].get_attribute("gridSpacing"))[::-1] * m[ + "e_all_chargeDensity" +].get_attribute("gridUnitSI") + +x = (np.arange(rho.shape[0]) - (rho.shape[0] - 1) / 2) * spacing[0] +y = (np.arange(rho.shape[1]) - (rho.shape[1] - 1) / 2) * spacing[1] +z = (np.arange(rho.shape[2]) - (rho.shape[2] - 1) / 2) * spacing[2] + +# collect electric field +E_x = m["E"]["x"][:, :, :] +E_y = m["E"]["y"][:, :, :] +E_z = m["E"]["z"][:, :, :] +series.flush() + +E_x *= m["E"]["x"].get_attribute("unitSI") +E_x = np.transpose(E_x) +E_y *= m["E"]["y"].get_attribute("unitSI") +E_y = np.transpose(E_y) +E_z *= m["E"]["z"].get_attribute("unitSI") +E_z = np.transpose(E_z) + +# get the divergence of the electric field, which should equal the charge density +diff = np.zeros(E_x.shape) +for ix in range(E_x.shape[0] - 2): + d0 = (E_x[ix + 1, :, :] - E_x[ix, :, :]) / spacing[0] + diff[ix + 1, :, :] += d0 + +for iy in range(E_y.shape[1] - 2): + d0 = (E_y[:, iy + 1, :] - E_y[:, iy, :]) / spacing[1] + diff[:, iy + 1, :] += d0 + +for iz in range(E_z.shape[2] - 2): + d0 = (E_z[:, :, iz + 1] - E_z[:, :, iz]) / spacing[2] + diff[:, :, iz + 1] += d0 + +diff *= epsilon_0 + +# comparing the two +rerr = np.abs((diff[1:-1, 1:-1, 1:-1] - rho[1:-1, 1:-1, 1:-1]) / diff[1:-1, 1:-1, 1:-1]) + +print("maximum relative error:", np.max(rerr)) +print( + "maximum relative error for absolute values above 1e-3 of the maximum charge density:", + np.max(rerr[np.abs(diff[1:-1, 1:-1, 1:-1]) > 1e-3 * np.max(np.abs(diff))]), +) + +if np.max(rerr[np.abs(diff[1:-1, 1:-1, 1:-1]) > 1e-3 * np.max(np.abs(diff))]) > 0.005: + print("relative error too high!") + sys.exit(1) + +if np.max(rerr) > 0.9: + print("relative error too high!") + sys.exit(1) + +sys.exit(0) From c31e9ecc5ab11b1469d9913e6ff4a44b653405e7 Mon Sep 17 00:00:00 2001 From: Edgar Marquardt Date: Thu, 30 Jul 2026 16:14:15 +0200 Subject: [PATCH 2/7] Add the poisson solver to pypicongpu and picmi There is now a way to use it in both frameworks. It probably needs some more documentation though. --- lib/python/picongpu/picmi/simulation.py | 6 ++++ .../picongpu/pypicongpu/poissonsolver.py | 30 +++++++++++++++++++ lib/python/picongpu/pypicongpu/simulation.py | 4 +++ .../templates/etc/picongpu/N.cfg.mustache | 19 ++++++++++++ 4 files changed, 59 insertions(+) create mode 100644 lib/python/picongpu/pypicongpu/poissonsolver.py diff --git a/lib/python/picongpu/picmi/simulation.py b/lib/python/picongpu/picmi/simulation.py index e84a3a07e9b..6b3b4c45cf7 100644 --- a/lib/python/picongpu/picmi/simulation.py +++ b/lib/python/picongpu/picmi/simulation.py @@ -43,6 +43,7 @@ from picongpu.pypicongpu.species.constant.synchrotron import SynchrotronParams from picongpu.pypicongpu.util import UnpackChain, unique from picongpu.pypicongpu.walltime import Walltime +from picongpu.pypicongpu.poissonsolver import PoissonSolver as PIConGPUPoissonSolver class _DensityImpl(BaseModel): @@ -199,6 +200,9 @@ def _validate_typical_ppc(value: int | None) -> int | None: picongpu_base_density: float | None = Field(default=None) """value to normalise densities with""" + picongpu_poisson_solver: PIConGPUPoissonSolver | None = Field(default=None) + """Poisson solver to use for electrostatic calculations for the starting conditions, set to None to disable""" + picongpu_walltime: datetime.timedelta | None = Field(default=None) """time after which the cluster scheduler will stop the simulation""" @@ -208,6 +212,7 @@ def _validate_typical_ppc(value: int | None) -> int | None: model_config = ConfigDict(arbitrary_types_allowed=True) + @model_validator(mode="after") def _post_init(self): # additional PICMI stuff checks, @todo move to picmistandard, Brian Marre, 2024 @@ -429,6 +434,7 @@ def get_as_pypicongpu(self) -> pypicongpu.simulation.Simulation: grid=self.solver.grid.get_as_pypicongpu(), binomial_current_interpolation=self.solver.source_smoother is not None, moving_window=moving_window, + poisson_solver=self.picongpu_poisson_solver, walltime=walltime or Walltime(walltime=datetime.timedelta(hours=1)), time_steps=time_steps, laser=[ll.get_as_pypicongpu() for ll in self.lasers] or None, diff --git a/lib/python/picongpu/pypicongpu/poissonsolver.py b/lib/python/picongpu/pypicongpu/poissonsolver.py new file mode 100644 index 00000000000..c83a9b70c05 --- /dev/null +++ b/lib/python/picongpu/pypicongpu/poissonsolver.py @@ -0,0 +1,30 @@ +""" +This file is part of the PIConGPU. +Copyright 2025-2026 PIConGPU contributors +Authors: Edgar Marquardt +License: GPLv3+ +""" + +from typing import Annotated + +from pydantic import BaseModel, Field + +from .rendering import RenderedObject + + +class PoissonSolver(RenderedObject, BaseModel): + """ + Poisson solver for the electric field in the starting condition. + """ + + max_steps: Annotated[int, Field(..., gt=0.0)] | None + """maximum number of iterations for the Poisson solver""" + + epsilon: Annotated[float, Field(..., gt=0.0)] | None + """tolerance for the Poisson solver""" + + preconditioner_disabled: Annotated[bool, Field(...)] | None + """disable preconditioner for the Poisson solver""" + + preconditioner_max_steps: Annotated[int, Field(..., gt=0.0)] | None + """maximum number of iterations for the preconditioner""" diff --git a/lib/python/picongpu/pypicongpu/simulation.py b/lib/python/picongpu/pypicongpu/simulation.py index d84817e5b93..00ca4d86c05 100644 --- a/lib/python/picongpu/pypicongpu/simulation.py +++ b/lib/python/picongpu/pypicongpu/simulation.py @@ -23,6 +23,7 @@ from .grid import Grid3D from .laser import AnyLaser from .movingwindow import MovingWindow +from .poissonsolver import PoissonSolver from .output import AnyPlugin, OpenPMDPlugin from .rendering import RenderedObject from .walltime import Walltime @@ -73,6 +74,9 @@ class Simulation(RenderedObject, BaseModel): moving_window: MovingWindow | None """used moving Window, set to None to disable""" + poisson_solver: PoissonSolver | None + """used poisson solver for electrostatic calculations for the starting conditions, set to None to disable""" + walltime: Walltime """time limit of the simulation run""" diff --git a/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache b/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache index 4dfe71e5cc8..bbe96c1ce0e 100644 --- a/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache +++ b/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache @@ -81,6 +81,22 @@ TBG_steps="{{{time_steps}}}" {{/stop_iteration}} {{/moving_window}} +{{#poisson_solver}} + TBG_poissonsolver="--poisson.activate" + {{#max_steps}} + TBG_poissonsolver="$TBG_poissonsolver --poisson.maxSteps {{{max_steps}}}" + {{/max_steps}} + {{#epsilon}} + TBG_poissonsolver="$TBG_poissonsolver --poisson.epsilon {{{epsilon}}}" + {{/epsilon}} + {{#preconditioner_disabled}} + TBG_poissonsolver="$TBG_poissonsolver --poisson.preconditioner.disable" + {{/preconditioner_disabled}} + {{#preconditioner_max_steps}} + TBG_poissonsolver="$TBG_poissonsolver --poisson.preconditioner.maxSteps {{{preconditioner_max_steps}}}" + {{/preconditioner_max_steps}} +{{/poisson_solver}} + {{#binomial_current_interpolation}} TBG_currentInterpolation="--currentInterpolation binomial" {{/binomial_current_interpolation}} @@ -226,6 +242,9 @@ TBG_programParams="-d !TBG_deviceDist \ !TBG_stopWindow \ {{/stop_iteration}} {{/moving_window}} + {{#poisson_solver}} + !TBG_poissonsolver \ + {{/poisson_solver}} --versionOnce" # TOTAL number of devices From 37df26abf9675077127d68880f2a1820cf3ad37f Mon Sep 17 00:00:00 2001 From: Edgar Marquardt Date: Mon, 7 Sep 2026 11:03:43 +0200 Subject: [PATCH 3/7] improved the interface --- docs/TBG_macros.cfg | 4 ++-- include/picongpu/simulation/stage/Poisson.hpp | 2 +- include/picongpu/simulation/stage/Poisson.x.cpp | 6 +++--- lib/python/picongpu/pypicongpu/poissonsolver.py | 10 +++++----- .../picongpu/templates/etc/picongpu/N.cfg.mustache | 6 +++--- 5 files changed, 14 insertions(+), 14 deletions(-) diff --git a/docs/TBG_macros.cfg b/docs/TBG_macros.cfg index db2f7c00bf8..a6a7ce4b6c3 100644 --- a/docs/TBG_macros.cfg +++ b/docs/TBG_macros.cfg @@ -115,9 +115,9 @@ TBG_runtimeDensityFile="--e_runtimeDensityFile /bigdata/hplsim/production/someDi # Activate the Poisson solver and set its parameters. The Poisson solver is used to calculate the electrostatic # potential from the charge density. -# --poisson.maxSteps sets the maximum number of iterations for the solver, and --poisson.epsilon sets the convergence +# --poisson.maxSteps sets the maximum number of iterations for the solver, and --poisson.tolerance sets the convergence # criterion: the maximum acceptable relative error. -TBG_poissonSolver="--poisson.activate --poisson.maxSteps 2000 --poisson.epsilon 1e-8" +TBG_poissonSolver="--poisson.activate --poisson.maxSteps 2000 --poisson.tolerance 1e-8" # The Poisson solver is by default preconditioned. The preconditioner can be disabled with --poisson.preconditioner.disable. # The maximum number of iterations for the preconditioner can be set with --poisson.preconditioner.maxSteps. diff --git a/include/picongpu/simulation/stage/Poisson.hpp b/include/picongpu/simulation/stage/Poisson.hpp index 4b5d9d98ddd..b296feb4c02 100644 --- a/include/picongpu/simulation/stage/Poisson.hpp +++ b/include/picongpu/simulation/stage/Poisson.hpp @@ -96,7 +96,7 @@ namespace picongpu::simulation::stage // defaults will be overwritten by command line arguments bool m_useSolver = false; uint32_t m_maxSolverSteps = 20; - float_64 m_solverEpsilon = 1.0e-8; + float_64 m_solverTolerance = 1.0e-8; bool m_disablePreconditioner = false; uint32_t m_maxPreconditionerSteps = 20; diff --git a/include/picongpu/simulation/stage/Poisson.x.cpp b/include/picongpu/simulation/stage/Poisson.x.cpp index 43ee543080a..4efc714475b 100644 --- a/include/picongpu/simulation/stage/Poisson.x.cpp +++ b/include/picongpu/simulation/stage/Poisson.x.cpp @@ -102,8 +102,8 @@ namespace picongpu po::value(&m_maxSolverSteps)->default_value(2000), "maximum number of steps for the solver"); solverDesc.add_options()( - "poisson.epsilon", - po::value(&m_solverEpsilon)->default_value(1.0e-8), + "poisson.tolerance", + po::value(&m_solverTolerance)->default_value(1.0e-8), "maximal allowed error of the poisson solver"); // preconitioner solverDesc.add_options()( @@ -586,7 +586,7 @@ namespace picongpu rho1 = totalSum1; float_64 beta = rho1 / rho0 * alpha / omega; rho0 = rho1; - if(std::sqrt(totalSum2) < m_solverEpsilon) + if(std::sqrt(totalSum2) < m_solverTolerance) { foundSolution = true; break; diff --git a/lib/python/picongpu/pypicongpu/poissonsolver.py b/lib/python/picongpu/pypicongpu/poissonsolver.py index c83a9b70c05..8f1f5042c5a 100644 --- a/lib/python/picongpu/pypicongpu/poissonsolver.py +++ b/lib/python/picongpu/pypicongpu/poissonsolver.py @@ -17,14 +17,14 @@ class PoissonSolver(RenderedObject, BaseModel): Poisson solver for the electric field in the starting condition. """ - max_steps: Annotated[int, Field(..., gt=0.0)] | None + max_steps: Annotated[int, Field(..., gt=0)] | None = None """maximum number of iterations for the Poisson solver""" - epsilon: Annotated[float, Field(..., gt=0.0)] | None - """tolerance for the Poisson solver""" + tolerance: Annotated[float, Field(..., gt=0.0)] | None = None + """maximum tolerance for the Poisson solver""" - preconditioner_disabled: Annotated[bool, Field(...)] | None + preconditioner_disabled: Annotated[bool, Field(...)] | None = None """disable preconditioner for the Poisson solver""" - preconditioner_max_steps: Annotated[int, Field(..., gt=0.0)] | None + preconditioner_max_steps: Annotated[int, Field(..., gt=0)] | None = None """maximum number of iterations for the preconditioner""" diff --git a/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache b/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache index bbe96c1ce0e..80f176d4c29 100644 --- a/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache +++ b/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache @@ -86,9 +86,9 @@ TBG_steps="{{{time_steps}}}" {{#max_steps}} TBG_poissonsolver="$TBG_poissonsolver --poisson.maxSteps {{{max_steps}}}" {{/max_steps}} - {{#epsilon}} - TBG_poissonsolver="$TBG_poissonsolver --poisson.epsilon {{{epsilon}}}" - {{/epsilon}} + {{#tolerance}} + TBG_poissonsolver="$TBG_poissonsolver --poisson.tolerance {{{tolerance}}}" + {{/tolerance}} {{#preconditioner_disabled}} TBG_poissonsolver="$TBG_poissonsolver --poisson.preconditioner.disable" {{/preconditioner_disabled}} From 4b96fa36bad47e42f2ccfcdb6a116a7061d123b9 Mon Sep 17 00:00:00 2001 From: Edgar Marquardt Date: Wed, 9 Sep 2026 09:41:57 +0200 Subject: [PATCH 4/7] removed unneccessary files and parts --- .../fields/poissonSolver/BICGStab.hpp | 38 ------------------- .../fields/poissonSolver/ChargeDeposition.hpp | 20 ---------- include/pmacc/algorithms/ForEachCell.hpp | 6 --- 3 files changed, 64 deletions(-) delete mode 100644 include/picongpu/fields/poissonSolver/BICGStab.hpp delete mode 100644 include/picongpu/fields/poissonSolver/ChargeDeposition.hpp diff --git a/include/picongpu/fields/poissonSolver/BICGStab.hpp b/include/picongpu/fields/poissonSolver/BICGStab.hpp deleted file mode 100644 index 85ef7c092f5..00000000000 --- a/include/picongpu/fields/poissonSolver/BICGStab.hpp +++ /dev/null @@ -1,38 +0,0 @@ -/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera - * - * 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 . - */ - -#pragma once - -#include "picongpu/defines.hpp" -#include "picongpu/fields/FieldTmpOperations.hpp" - -namespace picongpu::fields::poissonSolver -{ - struct BICGStab - { - // return residual - // return number of iterations - void operator()(FieldTmp& fieldV, FieldTmp& fiedlRho, MappingDesc* cellDescription) - { - // set boundary conditions on fieldV (Dirichlet or Neuman) - - // normalize the problem based on norm(fieldRho) - } - }; -} // namespace picongpu::fields::poissonSolver diff --git a/include/picongpu/fields/poissonSolver/ChargeDeposition.hpp b/include/picongpu/fields/poissonSolver/ChargeDeposition.hpp deleted file mode 100644 index e6b0bda3e46..00000000000 --- a/include/picongpu/fields/poissonSolver/ChargeDeposition.hpp +++ /dev/null @@ -1,20 +0,0 @@ -/* Copyright 2025 Tapish Narwal, Luca Pennati, Rene Widera - * - * 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 . - */ - -#pragma once diff --git a/include/pmacc/algorithms/ForEachCell.hpp b/include/pmacc/algorithms/ForEachCell.hpp index f88030d96b5..f278932912f 100644 --- a/include/pmacc/algorithms/ForEachCell.hpp +++ b/include/pmacc/algorithms/ForEachCell.hpp @@ -45,12 +45,6 @@ namespace pmacc::algorithms { return func(std::forward(args)...); } - - template - DEVICEONLY auto operator()(T&&... args) - { - return func(std::forward(args)...); - } }; template From 98fa7d34573fc5572a8d9bd43f5d272b0a69a8fd Mon Sep 17 00:00:00 2001 From: Edgar Marquardt Date: Wed, 9 Sep 2026 09:47:11 +0200 Subject: [PATCH 5/7] updated the test setup to an analytically solvable one --- share/picongpu/tests/PoissonSolver/README.rst | 4 +- share/picongpu/tests/PoissonSolver/bin/ci.sh | 2 +- .../include/picongpu/param/density.param | 24 +++++++- .../include/picongpu/param/fieldSolver.param | 59 +++++++++++++++++++ .../include/picongpu/param/particle.param | 4 +- .../include/picongpu/param/simulation.param | 2 +- .../picongpu/param/speciesDefinition.param | 23 +------- .../param/speciesInitialization.param | 6 +- 8 files changed, 91 insertions(+), 33 deletions(-) create mode 100644 share/picongpu/tests/PoissonSolver/include/picongpu/param/fieldSolver.param diff --git a/share/picongpu/tests/PoissonSolver/README.rst b/share/picongpu/tests/PoissonSolver/README.rst index 66135e56d21..222059c1a5a 100644 --- a/share/picongpu/tests/PoissonSolver/README.rst +++ b/share/picongpu/tests/PoissonSolver/README.rst @@ -2,7 +2,9 @@ Test for the Poisson solver =========================== -This test is for validating the PoissonSolver in the PIConGPU initialization. +This test is for validating the PoissonSolver in the PIConGPU initialization. The test setup is a thick, charged tube, +for which the electric field can be calculated analytically. The test compares the analytical solution with the numerical solution +from the Poisson solver. To run this test, one has to execute ci.sh with the location of the input and output directory. diff --git a/share/picongpu/tests/PoissonSolver/bin/ci.sh b/share/picongpu/tests/PoissonSolver/bin/ci.sh index 1bb7a63cbdb..e1bf72e1a90 100755 --- a/share/picongpu/tests/PoissonSolver/bin/ci.sh +++ b/share/picongpu/tests/PoissonSolver/bin/ci.sh @@ -124,7 +124,7 @@ if [ $ret_build -eq 0 ] ; then # run the simulation echo "Simulation path: " $simPath"/" - mpiexec -n 1 ../bin/picongpu -d 1 1 1 -g 64 64 64 --periodic 1 1 1 -s 4 \ + mpiexec -n 1 ../bin/picongpu -d 1 1 1 -g 128 128 128 --periodic 1 1 1 -s 4 \ --openPMD.period 8 --openPMD.ext bp --openPMD.file simData --poisson.activate ret_sim=$? diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/density.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/density.param index f354524edd6..d0ced8b0be6 100644 --- a/share/picongpu/tests/PoissonSolver/include/picongpu/param/density.param +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/density.param @@ -26,7 +26,27 @@ namespace picongpu { namespace densityProfiles { - /* definition of homogenous profile */ - using Homogenous = HomogenousImpl; + /* definition of the functor that defines the density profile */ + struct CoaxialWireFunctor + { + static constexpr float3_64 center_SI{1.28e-4, 0.0, 1.28e-4}; + + static constexpr float_64 innerRadius_SI = 2.0e-5; + static constexpr float_64 outerRadius_SI = 4.0e-5; + + HDINLINE float_X operator()(floatD_64 const& position_SI, float3_64 const& cellSize_SI) const + { + float_64 const dx = position_SI.x() - center_SI.x(); + float_64 const dz = position_SI.z() - center_SI.z(); + float_64 const radiusSquared = dx * dx + dz * dz; + + return float_X( + radiusSquared >= innerRadius_SI * innerRadius_SI && radiusSquared < outerRadius_SI * outerRadius_SI + && position_SI.y() > 16 * cellSize_SI.y() && position_SI.y() < (128 - 16) * cellSize_SI.y()); + } + }; + + /* definition of coaxial wire profile */ + using CoaxialWire = FreeFormulaImpl; } // namespace densityProfiles } // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/fieldSolver.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/fieldSolver.param new file mode 100644 index 00000000000..99e79eb34bb --- /dev/null +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/fieldSolver.param @@ -0,0 +1,59 @@ +/* Copyright 2013-2024 Axel Huebl, Heiko Burau, Rene Widera, Sergei Bastrakov, + * Klaus Steiniger, Lennert Sprenger + * + * 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 . + */ + +/** @file + * + * Configure the field solver. + * + * Select the numerical Maxwell solver (e.g. Yee's method). + * + * \attention + * Currently, the laser initialization in PIConGPU is implemented to work with the standard Yee solver. + * Using a solver of higher order will result in a slightly increased laser amplitude and energy than expected. + * + */ + +#pragma once + +#include "picongpu/fields/MaxwellSolver/Solvers.def" + +namespace picongpu +{ + namespace fields + { + /** FieldSolver + * + * Field Solver Selection (note <> for some solvers), all in namespace maxwellSolver: + * - Yee: Standard Yee solver approximating derivatives with respect to time and + * space by second order finite differences. + * - CKC: Cole-Karkkainen-Cowan Solver, Dispersion free solver in the direction of the smallest + * grid size. + * - Lehe<>: Num. Cherenkov free field solver in a chosen direction + * - ArbitraryOrderFDTD<4>: Solver using 4 neighbors to each direction to approximate + * *spatial* derivatives by finite differences. The number of neighbors can be changed from 4 to any positive, + * integer number. The order of the solver will be twice the number of neighbors in each direction. Yee's + * method is a special case of this using one neighbor to each direction. + * - Substepping: use the given Solver (Yee, etc.) and substep each time step by factor 4 + * - None: disable the vacuum update of E and B, including no J contribution to E + */ + using Solver = maxwellSolver::None; + + } // namespace fields +} // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/particle.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/particle.param index e11790099ad..5db8c18cf48 100644 --- a/share/picongpu/tests/PoissonSolver/include/picongpu/param/particle.param +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/particle.param @@ -41,7 +41,7 @@ namespace picongpu * be created / will be deleted * * unit: none */ - constexpr float_X MIN_WEIGHTING = 10.0; + constexpr float_X MIN_WEIGHTING = 1.0e-2; namespace startPosition { @@ -58,7 +58,7 @@ namespace picongpu * * unit: none */ - static constexpr uint32_t numParticlesPerCell = 2u; + static constexpr uint32_t numParticlesPerCell = ::picongpu::TYPICAL_PARTICLES_PER_CELL; }; /** Definition of start position functor that randomly distributes macro-particles within a cell. */ diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/simulation.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/simulation.param index e54091dbb9a..c88f0fd0fc9 100644 --- a/share/picongpu/tests/PoissonSolver/include/picongpu/param/simulation.param +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/simulation.param @@ -80,5 +80,5 @@ namespace picongpu * * Used internally for unit normalization. */ - constexpr uint32_t TYPICAL_PARTICLES_PER_CELL = 2u; + constexpr uint32_t TYPICAL_PARTICLES_PER_CELL = 50u; } // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesDefinition.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesDefinition.param index 72b3bdf7b44..d4b2785ec97 100644 --- a/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesDefinition.param +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesDefinition.param @@ -57,27 +57,6 @@ namespace picongpu /* define species electrons */ using PIC_Electrons = Particles; - /*--------------------------- ions -------------------------------------------*/ - - /* ratio relative to base charge and base mass */ - value_identifier(float_X, MassRatioIons, 1836.152672); - value_identifier(float_X, ChargeRatioIons, -1.0); - - /* ratio relative to BASE_DENSITY */ - value_identifier(float_X, DensityRatioIons, 1.0); - - using ParticleFlagsIons = MakeSeq_t< - particlePusher, - shape, - interpolation, - current, - massRatio, - chargeRatio, - densityRatio, - atomicNumbers>; - - /* define species ions */ - using PIC_Ions = Particles; /*########################### end species ####################################*/ @@ -86,6 +65,6 @@ namespace picongpu * List all defined particle species from above in this list * to make them available to the PIC algorithm. */ - using VectorAllSpecies = MakeSeq_t; + using VectorAllSpecies = MakeSeq_t; } // namespace picongpu diff --git a/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesInitialization.param b/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesInitialization.param index 610c2b56c5b..a706a185e93 100644 --- a/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesInitialization.param +++ b/share/picongpu/tests/PoissonSolver/include/picongpu/param/speciesInitialization.param @@ -38,10 +38,8 @@ namespace picongpu * * the functors are called in order (from first to last functor) */ - using InitPipeline = pmacc::mp_list< - - CreateDensity, - CreateDensity>; + using InitPipeline + = pmacc::mp_list>; } // namespace particles } // namespace picongpu From c203ce554c23b0abca9ec5fda517211c8220df4d Mon Sep 17 00:00:00 2001 From: Edgar Marquardt Date: Wed, 9 Sep 2026 10:33:03 +0200 Subject: [PATCH 6/7] update the picmi implementation to follow the standard better --- lib/python/picongpu/picmi/simulation.py | 11 +++---- lib/python/picongpu/picmi/solver.py | 29 ++++++++++++++++++- .../picongpu/pypicongpu/poissonsolver.py | 18 +++++++----- .../templates/etc/picongpu/N.cfg.mustache | 12 ++------ 4 files changed, 48 insertions(+), 22 deletions(-) diff --git a/lib/python/picongpu/picmi/simulation.py b/lib/python/picongpu/picmi/simulation.py index 6b3b4c45cf7..9e05cbe2e42 100644 --- a/lib/python/picongpu/picmi/simulation.py +++ b/lib/python/picongpu/picmi/simulation.py @@ -24,6 +24,7 @@ from picongpu.picmi.diagnostics.field_dump import NativeFieldDump, _FieldDump from picongpu.picmi.diagnostics.particle_dump import ParticleDump from picongpu.picmi.grid import Cartesian3DGrid +from picongpu.picmi.solver import ElectrostaticSolver from picongpu.picmi.interaction import Interaction, Synchrotron from picongpu.picmi.interaction.collision import Collision, CollisionalPhysicsSetup from picongpu.picmi.layout import AnyLayout @@ -43,7 +44,6 @@ from picongpu.pypicongpu.species.constant.synchrotron import SynchrotronParams from picongpu.pypicongpu.util import UnpackChain, unique from picongpu.pypicongpu.walltime import Walltime -from picongpu.pypicongpu.poissonsolver import PoissonSolver as PIConGPUPoissonSolver class _DensityImpl(BaseModel): @@ -200,8 +200,8 @@ def _validate_typical_ppc(value: int | None) -> int | None: picongpu_base_density: float | None = Field(default=None) """value to normalise densities with""" - picongpu_poisson_solver: PIConGPUPoissonSolver | None = Field(default=None) - """Poisson solver to use for electrostatic calculations for the starting conditions, set to None to disable""" + picongpu_electrostatic_solver: ElectrostaticSolver | None = Field(default=None) + """Electrostatic solver to use for electrostatic calculations for the starting conditions""" picongpu_walltime: datetime.timedelta | None = Field(default=None) """time after which the cluster scheduler will stop the simulation""" @@ -212,7 +212,6 @@ def _validate_typical_ppc(value: int | None) -> int | None: model_config = ConfigDict(arbitrary_types_allowed=True) - @model_validator(mode="after") def _post_init(self): # additional PICMI stuff checks, @todo move to picmistandard, Brian Marre, 2024 @@ -434,7 +433,9 @@ def get_as_pypicongpu(self) -> pypicongpu.simulation.Simulation: grid=self.solver.grid.get_as_pypicongpu(), binomial_current_interpolation=self.solver.source_smoother is not None, moving_window=moving_window, - poisson_solver=self.picongpu_poisson_solver, + poisson_solver=self.picongpu_electrostatic_solver.get_as_pypicongpu() + if self.picongpu_electrostatic_solver is not None + else None, walltime=walltime or Walltime(walltime=datetime.timedelta(hours=1)), time_steps=time_steps, laser=[ll.get_as_pypicongpu() for ll in self.lasers] or None, diff --git a/lib/python/picongpu/picmi/solver.py b/lib/python/picongpu/picmi/solver.py index 71000421f27..728313f0f7b 100644 --- a/lib/python/picongpu/picmi/solver.py +++ b/lib/python/picongpu/picmi/solver.py @@ -7,11 +7,13 @@ from collections.abc import Sequence from typing import Annotated, Literal +from pydantic import Field -from picmistandard import PICMI_BinomialSmoother, PICMI_ElectromagneticSolver +from picmistandard import PICMI_BinomialSmoother, PICMI_ElectromagneticSolver, PICMI_ElectrostaticSolver from picongpu.pypicongpu import util from picongpu.pypicongpu.field_solver import AnySolver, LeheSolver, YeeSolver +from picongpu.pypicongpu.poissonsolver import PoissonSolver class BinomialSmoother(PICMI_BinomialSmoother): @@ -50,3 +52,28 @@ class ElectromagneticSolver(PICMI_ElectromagneticSolver): def get_as_pypicongpu(self) -> AnySolver: return YeeSolver() if self.method == "Yee" else LeheSolver() + + +class ElectrostaticSolver(PICMI_ElectrostaticSolver): + """ + PICMI Electrostatic Solver + + See PICMI spec for full documentation. + + Only the Poisson solver is supported; solver options that PIConGPU + does not implement are rejected at construction time. + """ + + method: Literal["Poisson"] = "Poisson" + required_precision: Annotated[float, Field(..., gt=0.0)] = 1e-8 + maximum_iterations: Annotated[int, Field(..., gt=0)] = 2000 + preconditioner: Literal["default", "none"] = "default" + preconditioner_maximum_iterations: Annotated[int, Field(..., gt=0)] = 20 + + def get_as_pypicongpu(self) -> PoissonSolver: + return PoissonSolver( + tolerance=self.required_precision, + max_steps=self.maximum_iterations, + preconditioner=self.preconditioner, + preconditioner_max_steps=self.preconditioner_maximum_iterations, + ) diff --git a/lib/python/picongpu/pypicongpu/poissonsolver.py b/lib/python/picongpu/pypicongpu/poissonsolver.py index 8f1f5042c5a..44affdf5d56 100644 --- a/lib/python/picongpu/pypicongpu/poissonsolver.py +++ b/lib/python/picongpu/pypicongpu/poissonsolver.py @@ -5,9 +5,9 @@ License: GPLv3+ """ -from typing import Annotated +from typing import Annotated, Literal -from pydantic import BaseModel, Field +from pydantic import BaseModel, Field, computed_field from .rendering import RenderedObject @@ -17,14 +17,18 @@ class PoissonSolver(RenderedObject, BaseModel): Poisson solver for the electric field in the starting condition. """ - max_steps: Annotated[int, Field(..., gt=0)] | None = None + max_steps: Annotated[int, Field(..., gt=0)] = 2000 """maximum number of iterations for the Poisson solver""" - tolerance: Annotated[float, Field(..., gt=0.0)] | None = None + tolerance: Annotated[float, Field(..., gt=0.0)] = 1e-8 """maximum tolerance for the Poisson solver""" - preconditioner_disabled: Annotated[bool, Field(...)] | None = None - """disable preconditioner for the Poisson solver""" + preconditioner: Literal["default", "none"] = "default" + """preconditioner for the Poisson solver""" - preconditioner_max_steps: Annotated[int, Field(..., gt=0)] | None = None + preconditioner_max_steps: Annotated[int, Field(..., gt=0)] = 20 """maximum number of iterations for the preconditioner""" + + @computed_field + def preconditioner_disabled(self) -> bool: + return self.preconditioner == "none" diff --git a/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache b/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache index 80f176d4c29..524f708a713 100644 --- a/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache +++ b/lib/python/picongpu/templates/etc/picongpu/N.cfg.mustache @@ -83,18 +83,12 @@ TBG_steps="{{{time_steps}}}" {{#poisson_solver}} TBG_poissonsolver="--poisson.activate" - {{#max_steps}} - TBG_poissonsolver="$TBG_poissonsolver --poisson.maxSteps {{{max_steps}}}" - {{/max_steps}} - {{#tolerance}} - TBG_poissonsolver="$TBG_poissonsolver --poisson.tolerance {{{tolerance}}}" - {{/tolerance}} + TBG_poissonsolver="$TBG_poissonsolver --poisson.maxSteps {{{max_steps}}}" + TBG_poissonsolver="$TBG_poissonsolver --poisson.tolerance {{{tolerance}}}" {{#preconditioner_disabled}} TBG_poissonsolver="$TBG_poissonsolver --poisson.preconditioner.disable" {{/preconditioner_disabled}} - {{#preconditioner_max_steps}} - TBG_poissonsolver="$TBG_poissonsolver --poisson.preconditioner.maxSteps {{{preconditioner_max_steps}}}" - {{/preconditioner_max_steps}} + TBG_poissonsolver="$TBG_poissonsolver --poisson.preconditioner.maxSteps {{{preconditioner_max_steps}}}" {{/poisson_solver}} {{#binomial_current_interpolation}} From a620ed349b3aabf3ed8bcfdb83e2bb2a140d6813 Mon Sep 17 00:00:00 2001 From: Edgar Marquardt Date: Wed, 9 Sep 2026 14:07:13 +0200 Subject: [PATCH 7/7] made the picmi implementation run properly --- lib/python/picongpu/picmi/__init__.py | 3 +- lib/python/picongpu/picmi/solver.py | 38 ++++++++++++++++--- .../picongpu/pypicongpu/poissonsolver.py | 8 ++-- 3 files changed, 38 insertions(+), 11 deletions(-) diff --git a/lib/python/picongpu/picmi/__init__.py b/lib/python/picongpu/picmi/__init__.py index 0e24cbbb032..d4e7c0ac9b1 100644 --- a/lib/python/picongpu/picmi/__init__.py +++ b/lib/python/picongpu/picmi/__init__.py @@ -40,7 +40,7 @@ from .layout import GriddedLayout, OnePositionLayout, PseudoRandomLayout from .particle_functor import FilteredSpecies, ParticleFilter, ParticleFunctor from .simulation import Simulation -from .solver import BinomialSmoother, ElectromagneticSolver +from .solver import BinomialSmoother, ElectromagneticSolver, ElectrostaticSolver from .species import Species assert sys.version_info.major > 3 or sys.version_info.minor >= 11, "Python 3.11 is required for PIConGPU PICMI" @@ -50,6 +50,7 @@ "ParticleFunctor", "Cartesian3DGrid", "ElectromagneticSolver", + "ElectrostaticSolver", "BinomialSmoother", "DispersivePulseLaser", "FromOpenPMDPulseLaser", diff --git a/lib/python/picongpu/picmi/solver.py b/lib/python/picongpu/picmi/solver.py index 728313f0f7b..1edff1ebd14 100644 --- a/lib/python/picongpu/picmi/solver.py +++ b/lib/python/picongpu/picmi/solver.py @@ -1,15 +1,17 @@ """ This file is part of PIConGPU. -Copyright 2021-2024 PIConGPU contributors -Authors: Hannes Troepgen, Brian Edward Marre, Richard Pausch +Copyright 2021-2026 PIConGPU contributors +Authors: Hannes Troepgen, Brian Edward Marre, Richard Pausch, Edgar Marquardt License: GPLv3+ """ from collections.abc import Sequence -from typing import Annotated, Literal -from pydantic import Field +from typing import Annotated, Literal, get_args +from pydantic import Field, computed_field -from picmistandard import PICMI_BinomialSmoother, PICMI_ElectromagneticSolver, PICMI_ElectrostaticSolver +from picmistandard import PICMI_BinomialSmoother, PICMI_ElectromagneticSolver +from picmistandard.base import _PICMIModel +from picmistandard.fields import PICMI_AnyGrid from picongpu.pypicongpu import util from picongpu.pypicongpu.field_solver import AnySolver, LeheSolver, YeeSolver @@ -54,6 +56,31 @@ def get_as_pypicongpu(self) -> AnySolver: return YeeSolver() if self.method == "Yee" else LeheSolver() +class PICMI_ElectrostaticSolver(_PICMIModel): + """ + Electrostatic field solver + """ + + @computed_field + def methods_list(self) -> list[str]: + # Retained for backwards compatibility reasons. + # The type annotation of `method` is the ground-truth. + return list(get_args(type(self).__annotations__["method"])) + + grid: PICMI_AnyGrid = Field(description="Grid object for the diagnostic") + + method: Literal["FFT", "Multigrid"] | None = Field( + default=None, + description="The advance method use to solve the poisson equation. The default method is code dependent.", + ) + + required_precision: float | None = Field(default=None, description="The required precision for iterative solvers.") + + maximum_iterations: int | None = Field( + default=None, description="The maximum number of iterations for iterative solvers." + ) + + class ElectrostaticSolver(PICMI_ElectrostaticSolver): """ PICMI Electrostatic Solver @@ -64,7 +91,6 @@ class ElectrostaticSolver(PICMI_ElectrostaticSolver): does not implement are rejected at construction time. """ - method: Literal["Poisson"] = "Poisson" required_precision: Annotated[float, Field(..., gt=0.0)] = 1e-8 maximum_iterations: Annotated[int, Field(..., gt=0)] = 2000 preconditioner: Literal["default", "none"] = "default" diff --git a/lib/python/picongpu/pypicongpu/poissonsolver.py b/lib/python/picongpu/pypicongpu/poissonsolver.py index 44affdf5d56..efec24f71a9 100644 --- a/lib/python/picongpu/pypicongpu/poissonsolver.py +++ b/lib/python/picongpu/pypicongpu/poissonsolver.py @@ -17,16 +17,16 @@ class PoissonSolver(RenderedObject, BaseModel): Poisson solver for the electric field in the starting condition. """ - max_steps: Annotated[int, Field(..., gt=0)] = 2000 + max_steps: Annotated[int, Field(..., gt=0)] """maximum number of iterations for the Poisson solver""" - tolerance: Annotated[float, Field(..., gt=0.0)] = 1e-8 + tolerance: Annotated[float, Field(..., gt=0.0)] """maximum tolerance for the Poisson solver""" - preconditioner: Literal["default", "none"] = "default" + preconditioner: Literal["default", "none"] """preconditioner for the Poisson solver""" - preconditioner_max_steps: Annotated[int, Field(..., gt=0)] = 20 + preconditioner_max_steps: Annotated[int, Field(..., gt=0)] """maximum number of iterations for the preconditioner""" @computed_field