Added separated class for the charge density equation solver + file f…

include/ModelingMainDriver.hpp

@@ -43,7 +43,7 @@
 class ModelingMainDriver final
 {
 private:
-    std::string __expt_m_config_filename;
+    std::string m_config_filename; ///< Filename of the configuration to parse it.
 
     static constexpr short const kdefault_max_numparticles_to_anim{5'000}; ///< Maximal count of particles to do animation.
 
@@ -151,23 +151,6 @@ class ModelingMainDriver final
     template <typename Function, typename... Args>
     void _processWithThreads(unsigned int num_threads, Function &&function, std::launch launch_plicy, Args &&...args);
 
-    /**
-     * @brief Solves the equation for the system using node charge density and boundary conditions.
-     *
-     * This function solves the system of equations based on the provided node charge density map
-     * and boundary conditions. It updates the solution vector accordingly.
-     *
-     * @param nodeChargeDensityMap A reference to a map containing the charge density at each node.
-     * @param gsmAssembler A shared pointer to the GSM assembler that constructs the system of equations.
-     * @param solutionVector A shared pointer to the vector manager that holds the solution.
-     * @param boundaryConditions A reference to a map containing boundary conditions for the system.
-     * @param time The current simulation time.
-     */
-    void _solveEquation(std::map<GlobalOrdinal, double> &nodeChargeDensityMap,
-                        std::shared_ptr<GSMAssembler> &gsmAssembler,
-                        std::shared_ptr<VectorManager> &solutionVector,
-                        std::map<GlobalOrdinal, double> &boundaryConditions, double time);
-
     /**
      * @brief Processes particle-in-cell (PIC) and surface collision tracking within a particle range.
      *

include/ParticleInCellEngine/ChargeDensityEquationSolver.hpp

@@ -0,0 +1,44 @@
+#ifndef CHARGEDNSITYEQUATIONSOLVER_HPP
+#define CHARGEDNSITYEQUATIONSOLVER_HPP
+
+#include "FiniteElementMethod/BoundaryConditions/BoundaryConditionsManager.hpp"
+#include "FiniteElementMethod/GSMAssembler.hpp"
+#include "FiniteElementMethod/LinearAlgebraManagers/VectorManager.hpp"
+#include "PICTypes.hpp"
+#include "Utilities/ConfigParser.hpp"
+
+/**
+ * @class ChargeDensityEquationSolver
+ * @brief A static class providing methods to solve charge density equations.
+ *
+ * This class offers static methods to process charge density maps, apply boundary conditions, 
+ * and solve the resulting equations using a configurable matrix equation solver. The solution
+ * updates a vector manager, reflecting changes in node potentials and fields.
+ */
+class ChargeDensityEquationSolver {
+  public:
+   /**
+     * @brief Solves the charge density equation system at a given time moment.
+     *
+     * This method processes the provided node charge density map and boundary conditions, 
+     * constructs the system of equations using the specified assembler, and solves the system. 
+     * It also handles writing outputs such as electric potentials and fields to files.
+     *
+     * @param timeMoment The simulation time for which the system is solved.
+     * @param configFilename The configuration file specifying solver parameters and FEM settings.
+     * @param nodeChargeDensityMap A reference to a map containing charge density values per node.
+     * @param gsmAssembler A shared pointer to the GSM assembler that constructs the system matrix.
+     * @param solutionVector A shared pointer to the vector manager storing the computed solution.
+     * @param boundaryConditions A reference to a map defining boundary conditions for specific nodes.
+     *
+     * @throws std::runtime_error If a critical error occurs while parsing the configuration or solving equations.
+     * @throws std::exception For any generic error during processing.
+     */
+    static void solve(double timeMoment, std::string_view configFilename,
+                      NodeChargeDensitiesMap &nodeChargeDensityMap,
+                      std::shared_ptr<GSMAssembler> &gsmAssembler,
+                      std::shared_ptr<VectorManager> &solutionVector,
+                      BoundaryConditionsMap &boundaryConditions);
+};
+
+#endif // !CHARGEDNSITYEQUATIONSOLVER_HPP

include/ParticleInCellEngine/NodeChargeDensityProcessor.hpp

@@ -6,37 +6,9 @@
 
 #include "FiniteElementMethod/GSMAssembler.hpp"
 #include "Geometry/CubicGrid.hpp"
-#include "Particle/Particle.hpp"
+#include "ParticleInCellEngine/PICTypes.hpp"
 #include "Utilities/ConfigParser.hpp"
 
-/**
- * @brief Map of charge densities for nodes in a mesh.
- * @details This type associates each node ID with its corresponding charge density.
- *          The structure is as follows:
- *          - Key: Node ID (GlobalOrdinal)
- *          - Value: Charge density at the node (double)
- */
-using NodeChargeDensitiesMap = std::map<GlobalOrdinal, double>;
-
-/**
- * @brief Tracker for particles inside tetrahedrons over time.
- * @details This type organizes particles by simulation time and tetrahedron ID.
- *          The structure is as follows:
- *          - Key: Simulation time (double)
- *          - Value: Map of tetrahedron ID to the particles inside it:
- *            - Key: Tetrahedron ID (size_t)
- *            - Value: Vector of particles inside the tetrahedron (ParticleVector)
- */
-using ParticleTrackerMap = std::map<double, std::map<size_t, ParticleVector>>;
-
-/**
- * @brief Set of particle IDs that have settled on a 2D mesh.
- * @details This type stores the IDs of particles that are considered settled
- *          after colliding or interacting with the 2D mesh.
- *          - Element: Particle ID (size_t)
- */
-using ParticlesIDSet = std::set<size_t>;
-
 static std::mutex g_nodeChargeDensityMap_mutex;    ///< Mutex for synchronizing access to the charge densities in nodes.
 static std::mutex g_particleTrackerMap_mutex;      ///< Mutex for synchronizing access to the particles in tetrahedrons.
 static std::shared_mutex g_settledParticles_mutex; ///< Mutex for synchronizing access to settled particle IDs.

include/ParticleInCellEngine/PICTypes.hpp

@@ -0,0 +1,47 @@
+#ifndef PICTYPES_HPP
+#define PICTYPES_HPP
+
+#include <map>
+#include <set>
+
+#include "FiniteElementMethod/FEMTypes.hpp"
+#include "Particle/Particle.hpp"
+
+/**
+ * @brief Map of charge densities for nodes in a mesh.
+ * @details This type associates each node ID with its corresponding charge density.
+ *          The structure is as follows:
+ *          - Key: Node ID (GlobalOrdinal)
+ *          - Value: Charge density at the node (double)
+ */
+using NodeChargeDensitiesMap = std::map<GlobalOrdinal, double>;
+
+/**
+ * @brief Tracker for particles inside tetrahedrons over time.
+ * @details This type organizes particles by simulation time and tetrahedron ID.
+ *          The structure is as follows:
+ *          - Key: Simulation time (double)
+ *          - Value: Map of tetrahedron ID to the particles inside it:
+ *            - Key: Tetrahedron ID (size_t)
+ *            - Value: Vector of particles inside the tetrahedron (ParticleVector)
+ */
+using ParticleTrackerMap = std::map<double, std::map<size_t, ParticleVector>>;
+
+/**
+ * @brief Set of particle IDs that have settled on a 2D mesh.
+ * @details This type stores the IDs of particles that are considered settled
+ *          after colliding or interacting with the 2D mesh.
+ *          - Element: Particle ID (size_t)
+ */
+using ParticlesIDSet = std::set<size_t>;
+
+/**
+ * @brief Map of the boundary conditions.
+ * @details This type associates each node ID with its corresponding boundary condition value.
+ *          The structure is as follows:
+ *          - Key: Node ID (GlobalOrdinal)
+ *          - Value: Value of the boundary conditions (double)
+ */
+using BoundaryConditionsMap = std::map<GlobalOrdinal, double>;
+
+#endif // !PICTYPES_HPP

src/ModelingMainDriver.cpp

@@ -14,6 +14,7 @@ using json = nlohmann::json;
 #include "ModelingMainDriver.hpp"
 #include "Particle/CUDA/ParticleDeviceMemoryConverter.cuh"
 #include "Particle/PhysicsCore/CollisionModel/CollisionModelFactory.hpp"
+#include "ParticleInCellEngine/ChargeDensityEquationSolver.hpp"
 #include "ParticleInCellEngine/NodeChargeDensityProcessor.hpp"
 
 std::mutex ModelingMainDriver::m_PICTracker_mutex;
@@ -270,7 +271,7 @@ void ModelingMainDriver::_processWithThreads(unsigned int num_threads, Function
 }
 
 ModelingMainDriver::ModelingMainDriver(std::string_view config_filename) 
-    : __expt_m_config_filename(config_filename), 
+    : m_config_filename(config_filename), 
       m_config(config_filename)
 {
     // Checking mesh filename on validity and assign it to the class member.
@@ -289,44 +290,6 @@ ModelingMainDriver::ModelingMainDriver(std::string_view config_filename)
 
 ModelingMainDriver::~ModelingMainDriver() { _gfinalize(); }
 
-void ModelingMainDriver::_solveEquation(std::map<GlobalOrdinal, double> &nodeChargeDensityMap,
-                                        std::shared_ptr<GSMAssembler> &gsmAssembler,
-                                        std::shared_ptr<VectorManager> &solutionVector,
-                                        std::map<GlobalOrdinal, double> &boundaryConditions, double time)
-{
-    try
-    {
-        auto nonChangebleNodes{m_config.getNonChangeableNodes()};
-        for (auto const &[nodeId, nodeChargeDensity] : nodeChargeDensityMap)
-#if __cplusplus >= 202002L
-            if (std::ranges::find(nonChangebleNodes, nodeId) == nonChangebleNodes.cend())
-#else
-            if (std::find(nonChangebleNodes.cbegin(), nonChangebleNodes.cend(), nodeId) == nonChangebleNodes.cend())
-#endif
-                boundaryConditions[nodeId] = nodeChargeDensity;
-        BoundaryConditionsManager::set(solutionVector->get(), FEM_LIMITS_DEFAULT_POLYNOMIAL_ORDER, boundaryConditions);
-
-        MatrixEquationSolver solver(gsmAssembler, solutionVector);
-        auto solverParams{solver.createSolverParams(m_config.getSolverName(), m_config.getMaxIterations(), m_config.getConvergenceTolerance(),
-                                                    m_config.getVerbosity(), m_config.getOutputFrequency(), m_config.getNumBlocks(), m_config.getBlockSize(),
-                                                    m_config.getMaxRestarts(), m_config.getFlexibleGMRES(), m_config.getOrthogonalization(),
-                                                    m_config.getAdaptiveBlockSize(), m_config.getConvergenceTestFrequency())};
-        solver.solve(m_config.getSolverName(), solverParams);
-        solver.calculateElectricField(); // Getting electric field for the each cell.
-
-        solver.writeElectricPotentialsToPosFile(time);
-        solver.writeElectricFieldVectorsToPosFile(time);
-    }
-    catch (std::exception const &ex)
-    {
-        ERRMSG(util::stringify("Can't solve the equation: ", ex.what()));
-    }
-    catch (...)
-    {
-        ERRMSG("Some error occured while solving the matrix equation Ax=b");
-    }
-}
-
 void ModelingMainDriver::_processPIC_and_SurfaceCollisionTracker(size_t start_index, size_t end_index, double t,
                                                                  std::shared_ptr<CubicGrid> cubicGrid, std::shared_ptr<GSMAssembler> gsmAssembler)
 {
@@ -598,7 +561,7 @@ void ModelingMainDriver::startModeling()
 #endif
 
         NodeChargeDensityProcessor::gather(t,
-                                           __expt_m_config_filename,
+                                           m_config_filename,
                                            cubicGrid,
                                            gsmAssembler,
                                            m_particles,
@@ -607,7 +570,12 @@ void ModelingMainDriver::startModeling()
                                            nodeChargeDensityMap);
 
         // 2. Solve equation in the main thread.
-        _solveEquation(nodeChargeDensityMap, gsmAssembler, solutionVector, boundaryConditions, t);
+        ChargeDensityEquationSolver::solve(t,
+                                           m_config_filename,
+                                           nodeChargeDensityMap, 
+                                           gsmAssembler, 
+                                           solutionVector,
+                                           boundaryConditions);
 
         // 3. Process surface collision tracking in parallel.
 #ifdef USE_OMP

src/ParticleInCellEngine/CMakeLists.txt

@@ -1,6 +1,7 @@
 # src/ParticleInCellEngine/CMakeLists.txt
 
 set(PARTICLE_IN_CELL_ENGINE_SOURCES
+    ChargeDensityEquationSolver.cpp
     NodeChargeDensityProcessor.cpp
 )
 

src/ParticleInCellEngine/ChargeDensityEquationSolver.cpp

@@ -0,0 +1,49 @@
+#include "ParticleInCellEngine/ChargeDensityEquationSolver.hpp"
+#include "FiniteElementMethod/FEMLimits.hpp"
+#include "FiniteElementMethod/MatrixEquationSolver.hpp"
+
+void ChargeDensityEquationSolver::solve(
+    double timeMoment, std::string_view configFilename,
+    NodeChargeDensitiesMap &nodeChargeDensityMap,
+    std::shared_ptr<GSMAssembler> &gsmAssembler,
+    std::shared_ptr<VectorManager> &solutionVector,
+    BoundaryConditionsMap &boundaryConditions) {
+    try {
+        ConfigParser configParser(configFilename);
+
+        auto nonChangebleNodes{configParser.getNonChangeableNodes()};
+        for (auto const &[nodeId, nodeChargeDensity] : nodeChargeDensityMap)
+#if __cplusplus >= 202002L
+            if (std::ranges::find(nonChangebleNodes, nodeId) ==
+                nonChangebleNodes.cend())
+#else
+            if (std::find(nonChangebleNodes.cbegin(), nonChangebleNodes.cend(),
+                          nodeId) == nonChangebleNodes.cend())
+#endif
+                boundaryConditions[nodeId] = nodeChargeDensity;
+        BoundaryConditionsManager::set(solutionVector->get(),
+                                       FEM_LIMITS_DEFAULT_POLYNOMIAL_ORDER,
+                                       boundaryConditions);
+
+        MatrixEquationSolver solver(gsmAssembler, solutionVector);
+        auto solverParams{solver.createSolverParams(
+            configParser.getSolverName(), configParser.getMaxIterations(),
+            configParser.getConvergenceTolerance(), configParser.getVerbosity(),
+            configParser.getOutputFrequency(), configParser.getNumBlocks(),
+            configParser.getBlockSize(), configParser.getMaxRestarts(),
+            configParser.getFlexibleGMRES(),
+            configParser.getOrthogonalization(),
+            configParser.getAdaptiveBlockSize(),
+            configParser.getConvergenceTestFrequency())};
+        solver.solve(configParser.getSolverName(), solverParams);
+        solver.calculateElectricField(); // Getting electric field for the
+                                         // each cell.
+
+        solver.writeElectricPotentialsToPosFile(timeMoment);
+        solver.writeElectricFieldVectorsToPosFile(timeMoment);
+    } catch (std::exception const &ex) {
+        ERRMSG(util::stringify("Can't solve the equation: ", ex.what()));
+    } catch (...) {
+        ERRMSG("Some error occured while solving the matrix equation Ax=b");
+    }
+}