add poisson solver to generate electrical field from non neutral charged gas densities

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

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"

include/picongpu/fields/poissonSolver/BICGStab.hpp

@@ -0,0 +1,35 @@
+/* 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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragme once
+
+#include "picongpu/defines.hpp"
+#include "picongpu/fields/FieldTmpOperations.hpp"
+
+struct BICGStab
+{
+    // return residual
+    // return number of iterations
+    void operator()(FieldTmp& fieldV, FieldTmp& fiedlRho, MappingDesk* cellDescription)
+    {
+        // set boundary conditions on fieldV (Dirichlet or Neuman)
+
+        // normalize the problem based on norm(fieldRho)
+    }
+};

include/picongpu/fields/poissonSolver/BoundaryConditions.hpp

@@ -0,0 +1,118 @@
+/* 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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "picongpu/defines.hpp"
+#include "picongpu/fields/FieldTmpOperations.hpp"
+#include "picongpu/fields/poissonSolver/FieldV.hpp"
+
+#include <pmacc/lockstep/lockstep.hpp>
+#include <pmacc/mappings/kernel/ExchangeMapping.hpp>
+#include <pmacc/memory/dataTypes/Mask.hpp>
+
+namespace picongpu::fields::poissonSolver
+{
+    struct SolutionFunction
+    {
+        HDINLINE auto operator()(math::Vector<double, simDim> const& totalCellCoordinate) const
+        {
+            if constexpr(simDim == 3u)
+            {
+                return math::sin(totalCellCoordinate.x()) + math::cos(totalCellCoordinate.y())
+                       + 3.0 * math::sin(totalCellCoordinate.z())
+                       + totalCellCoordinate.x() * totalCellCoordinate.productOfComponents() + 10.0;
+            }
+            else if constexpr(simDim == 2u)
+            {
+                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<simDim> cellOffsetToTotalOrigin,
+            auto const mapper) const -> void
+        {
+            // including guards
+            DataSpace<simDim> const superCellIdx(mapper.getSuperCellIndex(worker.blockDomIdxND()));
+
+            DataSpace<simDim> numGuardCells = mapper.getGuardingSuperCells() * SuperCellSize::toRT();
+
+            // no guards included
+            DataSpace<simDim> superCellTotalCellOffset
+                = cellOffsetToTotalOrigin + superCellIdx * SuperCellSize::toRT() - numGuardCells;
+
+            constexpr uint32_t cellsPerSuperCell = pmacc::math::CT::volume<SuperCellSize>::type::value;
+
+            auto forEachCellInSupercell = lockstep::makeForEach<cellsPerSuperCell>(worker);
+
+            forEachCellInSupercell(
+                [&](int32_t const linearCellIdx)
+                {
+                    /* cell index within the superCell */
+                    DataSpace<simDim> const cellIdx = pmacc::math::mapToND(SuperCellSize::toRT(), linearCellIdx);
+                    // without guards
+                    DataSpace<simDim> const totalCellIdx = superCellTotalCellOffset + cellIdx;
+
+                    auto totalDistance = precisionCast<float_64>(totalCellIdx)
+                                         * precisionCast<float_64>(sim.pic.getCellSize().shrink<simDim>());
+
+                    fieldVBox(superCellIdx * SuperCellSize::toRT() + cellIdx) = boundaryFunction(totalDistance);
+                });
+        }
+    };
+
+    struct BoundaryConditionsDirichlet
+    {
+        // return residual
+        // return number of iterations
+        void operator()(FieldV& fieldV, MappingDesc cellDescription) const
+        {
+            SubGrid<simDim> const& subGrid = Environment<simDim>::get().SubGrid();
+            auto globalDomain = subGrid.getGlobalDomain();
+            auto localDomain = subGrid.getLocalDomain();
+
+            auto cellOffsetToTotalOrigin = globalDomain.offset + localDomain.offset;
+
+
+            for(uint32_t i = 1; i < NumberOfExchanges<simDim>::value; ++i)
+            {
+                /* only call for planes: left right top bottom back front*/
+                if(FRONT % i == 0 && !(Environment<simDim>::get().GridController().getCommunicationMask().isSet(i)))
+                {
+                    ExchangeMapping<GUARD, MappingDesc> mapper(cellDescription, i);
+
+                    PMACC_LOCKSTEP_KERNEL(ApplyDirichletBCsFromFunctionKernel{})
+                        .config(mapper.getGridDim(), SuperCellSize{})(
+                            fieldV.fieldVBuffer->getDeviceBuffer().getDataBox(),
+                            SolutionFunction{},
+                            cellOffsetToTotalOrigin,
+                            mapper);
+                }
+            }
+        }
+    };
+} // namespace picongpu::fields::poissonSolver

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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once

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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "picongpu/defines.hpp"
+#include "picongpu/fields/YeeCell.hpp"
+#include "picongpu/traits/FieldPosition.hpp"
+#include "picongpu/traits/SIBaseUnits.hpp"
+
+#include <pmacc/memory/buffers/GridBuffer.hpp>
+
+namespace picongpu::fields::poissonSolver
+{
+    struct FieldV : pmacc::ISimulationData
+    {
+        using ValueType = float_64;
+
+        std::shared_ptr<GridBuffer<ValueType, simDim>> fieldVBuffer;
+
+        //! Unit type of field components
+        using UnitValueType = float1_64;
+
+        using DataBoxType = pmacc::DataBox<PitchedBox<ValueType, simDim>>;
+
+        //! Number of components of ValueType, for serialization
+        static constexpr int numComponents = 1u;
+
+        FieldV(picongpu::MappingDesc const mappingDesc)
+            : fieldVBuffer{std::make_shared<GridBuffer<float_64, simDim>>(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<simDim> 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<float_64> getUnitDimension()
+        {
+            /* V is in volts: V  = kg * m^2 / (A * s^3)
+             *   -> L^2 * M * T^-3 * I^-1
+             */
+            std::vector<float_64> 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<fields::YeeCell, fields::poissonSolver::FieldV, simDim>
+    {
+        using VectorVectorDD = ::pmacc::math::Vector<floatD_X, simDim> const;
+
+        HDINLINE FieldPosition() = default;
+
+        HDINLINE VectorVectorDD operator()() const
+        {
+            return VectorVectorDD::create(floatD_X::create(0.0));
+        }
+    };
+} // namespace picongpu::traits