Merge pull request #105 from LemurPwned/fix/electrical-connection

Reworking optional synchronisation option

Author: LemurPwned

CHANGELOG.md

@@ -1,5 +1,17 @@
 # Changelog
 
+# 1.8.0
+
+- Fixing bugs in electrical coupling
+- You can now have useKCL = True Stacks which will implicitly preserve Kirchoff's current law, and will simulate perfect voltage sources
+
+# 1.7.0
+
+- Adaptive simulation is now available again with Dormand-Price method
+- Adding curated simulation examples
+- Memory cells are now supported (you can specify two polarisation vectors per Layer now)
+- Import fixes and other minor bugfixes
+
 # 1.6.2-1.6.6
 
 - Fixed a bug in the `ScalarDriver` class which prevented the use of custom drivers with no arguments.

cmtj/stack/__init__.pyi

@@ -9,6 +9,7 @@ class ParallelStack:
         topId: str = "free",
         bottomId: str = "bottom",
         phaseOffset: float = 0,
+        useKCL: bool = True,
     ) -> None:
         """
         Initialises a parallel connection of junctions.
@@ -104,6 +105,7 @@ class SeriesStack:
         topId: str = "free",
         bottomId: str = "bottom",
         phaseOffset: float = 0,
+        useKCL: bool = True,
     ) -> None:
         """
         Initialises a series connection of junctions.

core/junction.hpp

@@ -108,6 +108,19 @@ template <typename T> class EnergyDriver {
   }
 };
 
+// Forward declarations if needed
+template <typename T> class Layer;
+template <typename T> class Junction;
+
+// Typedef declarations outside the class
+template <typename T>
+using SolverFn = bool (Layer<T>::*)(T t, T &timeStep, const CVector<T> &bottom,
+                                    const CVector<T> &top);
+
+template <typename T>
+using RunnerFn = void (Junction<T>::*)(SolverFn<T> &functor, T &t, T &timeStep,
+                                       bool &step_accepted);
+
 template <typename T> struct AdaptiveIntegrationParams {
   T abs_tol = 1e-6;             // Absolute error tolerance
   T rel_tol = 1e-3;             // Relative error tolerance
@@ -1641,10 +1654,6 @@ template <typename T> class Junction {
     logFile.close();
   }
 
-  typedef bool (Layer<T>::*solverFn)(T t, T &timeStep, const CVector<T> &bottom,
-                                     const CVector<T> &top);
-  typedef void (Junction<T>::*runnerFn)(solverFn &functor, T &t, T &timeStep,
-                                        bool &step_accepted);
   /**
    * @brief Run Euler-Heun or RK4 method for a single layer.
    *
@@ -1655,7 +1664,7 @@ template <typename T> class Junction {
    * @param t: current time
    * @param timeStep: integration step
    */
-  void runSingleLayerSolver(solverFn &functor, T &t, T &timeStep,
+  void runSingleLayerSolver(SolverFn<T> &functor, T &t, T &timeStep,
                             bool &step_accepted) {
     const CVector<T> dummy(0, 0, 0);
     step_accepted = (layers[0].*functor)(t, timeStep, dummy, dummy);
@@ -1669,7 +1678,7 @@ template <typename T> class Junction {
    * @param t: current time
    * @param timeStep: integration step
    * */
-  void runMultiLayerSolver(solverFn &functor, T &t, T &timeStep,
+  void runMultiLayerSolver(SolverFn<T> &functor, T &t, T &timeStep,
                            bool &step_accepted) {
     // Run solver for each layer and check if all steps were accepted
     step_accepted = true;
@@ -1696,7 +1705,7 @@ template <typename T> class Junction {
     }
   }
 
-  void eulerHeunSolverStep(solverFn &functor, T &t, T &timeStep,
+  void eulerHeunSolverStep(SolverFn<T> &functor, T &t, T &timeStep,
                            bool &step_accepted) {
     /*
         Euler Heun method (stochastic heun)
@@ -1743,7 +1752,7 @@ template <typename T> class Junction {
     }
   }
 
-  void heunSolverStep(solverFn &functor, T &t, T &timeStep,
+  void heunSolverStep(SolverFn<T> &functor, T &t, T &timeStep,
                       bool &step_accepted) {
     /*
         Heun method
@@ -1878,7 +1887,7 @@ template <typename T> class Junction {
     }
   }
 
-  std::tuple<runnerFn, solverFn, SolverMode>
+  std::tuple<RunnerFn<T>, SolverFn<T>, SolverMode>
   getSolver(SolverMode mode, unsigned int totalIterations) {
     SolverMode localMode = mode;
     for (auto &l : this->layers) {

core/reservoir.hpp

@@ -32,12 +32,6 @@ void comb(int N, int K) {
 
 typedef std::array<CVector<double>, 3> tensor;
 typedef std::vector<tensor> tensorList;
-template <typename T>
-using solverFn = bool (Layer<T>::*)(T t, T &timeStep, const CVector<T> &bottom,
-                                    const CVector<T> &top);
-template <typename T>
-using runnerFn = void (Junction<T>::*)(solverFn<T> &functor, T &t, T &timeStep,
-                                       bool &step_accepted);
 
 const tensor getDipoleTensorFromRelPositions(const CVector<double> &r1,
                                              const CVector<double> &r2) {
@@ -134,7 +128,7 @@ class GroupInteraction {
 
   void stepFunctionalSolver(double time, double &timeStep,
                             interactionFunction interaction,
-                            runnerFn<double> runner, solverFn<double> solver,
+                            RunnerFn<double> runner, SolverFn<double> solver,
                             bool &step_accepted) {
     // collect all frozen states
     // for each element, compute the extra field from all other elements