let's see what tests fail now

CHANGELOG.rst

@@ -30,6 +30,8 @@ Change Log
 - [FIXED] a bug in init from pypowsybl when some object were disconnected. It raises
   an error (because they are not connected to a bus): now this function properly handles
   these cases.
+- [FIXED] a bug leading to not propagate correctly the "compute_results" flag when the 
+  environment was copied (for example)
 - [ADDED] sets of methods to extract the main component of a grid and perform powerflow only on this
   one.
 - [ADDED] possibility to set / retrieve the names of each elements of the grid.

benchmarks/benchmark_solvers.py

@@ -43,10 +43,10 @@
                 lightsim2grid.SolverType.SparseLU: "NR (SLU)",
                 lightsim2grid.SolverType.KLU: "NR (KLU)",
                 lightsim2grid.SolverType.NICSLU: "NR (NICSLU *)",
+                lightsim2grid.SolverType.CKTSO: "NR (CKTSO *)",
                 lightsim2grid.SolverType.SparseLUSingleSlack: "NR single (SLU)",
                 lightsim2grid.SolverType.KLUSingleSlack: "NR single (KLU)",
                 lightsim2grid.SolverType.NICSLUSingleSlack: "NR single (NICSLU *)",
-                lightsim2grid.SolverType.CKTSO: "NR (CKTSO *)",
                 lightsim2grid.SolverType.CKTSOSingleSlack: "NR single (CKTSO *)",
                 lightsim2grid.SolverType.FDPF_XB_SparseLU: "FDPF XB (SLU)",
                 lightsim2grid.SolverType.FDPF_BX_SparseLU: "FDPF BX (SLU)",

lightsim2grid/tests/test_GridModel.py

@@ -103,15 +103,16 @@ def make_v0(self, net):
         return V0
 
     def run_me_pf(self, V0):
-        return self.model.compute_newton(V0, self.max_it, self.tol)
+        res = self.model.ac_pf(V0, self.max_it, self.tol)
+        return res
 
     def run_ref_pf(self, net):
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore")
             pp.runpp(net,
                      init="flat",
                      lightsim2grid=False,
-                     numba=True, 
+                     numba=False, 
                      distributed_slack=False)
 
     def do_i_skip(self, func_name):

lightsim2grid/tests/test_SameResPP.py

@@ -131,7 +131,7 @@ def _aux_test(self, pn_net):
 
         V = backend._grid.ac_pf(v_tmp, 10, 1e-5)
         assert V.shape[0], "? lightsim diverge when initialized with pp final voltage ?"
-        backend._grid.tell_topo_changed()
+        backend._grid.tell_solver_need_reset()
 
         Y_pp = backend.init_pp_backend._grid._ppc["internal"]["Ybus"]
         Sbus = backend.init_pp_backend._grid._ppc["internal"]["Sbus"]
@@ -218,7 +218,7 @@ def _aux_test(self, pn_net):
         backend._grid.deactivate_result_computation()
         Vdc = backend._grid.dc_pf(Vinit, max_iter, tol_this)
         backend._grid.reactivate_result_computation()
-        backend._grid.tell_topo_changed()
+        backend._grid.tell_solver_need_reset()
         Ydc_me = copy.deepcopy(backend._grid.get_dcYbus())
         Sdc_me = copy.deepcopy(backend._grid.get_dcSbus())
         assert np.max(np.abs(V_init_ref[pp_vect_converter] - Vdc[:nb_sub])) <= 100.*self.tol,\

lightsim2grid/tests/test_case118.py

@@ -131,7 +131,7 @@ def test_neurips_track2(self):
             with warnings.catch_warnings():
                 warnings.filterwarnings("ignore")
                 ls_grid = init(self.pp_net)
-            ls_grid.tell_topo_changed()
+            ls_grid.tell_solver_need_reset()
             V = np.ones(2 * self.nb_bus_total, dtype=np.complex_)
             V = ls_grid.ac_pf(V, self.max_it, self.tol)
             self.check_results(V[:self.nb_bus_total], ls_grid, self.pp_net)

lightsim2grid/tests/test_issue_56.py

@@ -49,7 +49,7 @@ def test_dc(self):
                 grid_model.deactivate_powerline(l_id)
             else:
                 grid_model.deactivate_trafo(l_id - nb_powerline)
-            grid_model.tell_topo_changed()
+            grid_model.tell_solver_need_reset()
             V = 1.0 * self.env.backend.V  # np.ones(2 * self.env.n_sub, dtype=np.complex_)
             res = grid_model.dc_pf(V, 10, 1e-8)
             if len(res):

lightsim2grid/tests/test_ptdf.py

@@ -36,7 +36,8 @@ def setUp(self) -> None:
         if solver_type not in self.gridmodel.available_solvers():
             self.skipTest("Solver type not supported on this platform")
         self.gridmodel.change_solver(solver_type)
-        self.gridmodel.dc_pf(self.V_init, 1, 1e-8)
+        V = self.gridmodel.dc_pf(self.V_init, 1, 1e-8)
+        assert len(V), f"dc pf has diverged with error {self.gridmodel.get_dc_solver().get_error()}"
         self.dcYbus = 1.0 * self.gridmodel.get_dcYbus()
         self.dcSbus = 1.0 * self.gridmodel.get_dcSbus().real
         self.Bbus = 1.0 * self.dcYbus.real

src/BaseNRSolver.h

@@ -136,8 +136,16 @@ class BaseNRSolver : public BaseSolver
 
         void fill_value_map(Eigen::Index slack_bus_id,
                             const Eigen::VectorXi & pq,
-                            const Eigen::VectorXi & pvpq);
-
+                            const Eigen::VectorXi & pvpq,
+                            bool reset_J);
+
+        void reset_if_needed(){
+            if(_solver_control.need_reset_solver() || 
+               _solver_control.has_dimension_changed() ||
+               _solver_control.has_ybus_some_coeffs_zero()){
+               reset();
+            }
+        }
     protected:
         // used linear solver
         LinearSolver _linear_solver;
@@ -151,6 +159,8 @@ class BaseNRSolver : public BaseSolver
         // to store the mapping from the element of J_ in dS_dVm_ and dS_dVa_
         // it does not own any memory at all !
         std::vector<cplx_type*> value_map_;
+        // std::vector<int> col_map_;
+        // std::vector<int> row_map_;
 
         // timers
         double timer_initialize_;

src/BaseNRSolver.tpp

@@ -46,20 +46,14 @@ bool BaseNRSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
         exc_ << "V  (" << V.size() << ") and Ybus (" << Ybus.rows()<< ", " << Ybus.cols() << ").";
         throw std::runtime_error(exc_.str());
     }
-    reset_timer();
-    // std::cout << "dist slack" << std::endl;
-
-    if(_solver_control.need_reset_solver() || 
-       _solver_control.has_dimension_changed()){
-        reset();
-    }
-    auto timer = CustTimer();
     if(!is_linear_solver_valid()) {
         // err_ = ErrorType::NotInitError;
         return false;
     }
-
+    reset_timer();
+    reset_if_needed();
     err_ = ErrorType::NoError;  // reset the error if previous error happened
+    auto timer = CustTimer();
 
     Eigen::VectorXi my_pv = retrieve_pv_with_slack(slack_ids, pv);  // retrieve_pv_with_slack (not all), add_slack_to_pv (all)
     real_type slack_absorbed = std::real(Sbus.sum());  // initial guess for slack_absorbed
@@ -94,9 +88,11 @@ bool BaseNRSolver<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type>
     bool has_just_been_initialized = false;  // to avoid a call to klu_refactor follow a call to klu_factor in the same loop
     // std::cout << "iter " << nr_iter_ << " dx(0): " << -F(0) << " dx(1): " << -F(1) << std::endl;
     // std::cout << "slack_absorbed " << slack_absorbed << std::endl;
-    BaseNRSolver<LinearSolver>::value_map_.clear();  // TODO smarter solver: only needed if ybus has changed
-    BaseNRSolver<LinearSolver>::dS_dVm_.resize(0,0);  // TODO smarter solver: only needed if ybus has changed
-    BaseNRSolver<LinearSolver>::dS_dVa_.resize(0,0);  // TODO smarter solver: only needed if ybus has changed
+    value_map_.clear();  // TODO smarter solver: only needed if ybus has changed
+    // BaseNRSolver<LinearSolver>::col_map_.clear();  // TODO smarter solver: only needed if ybus has changed
+    // BaseNRSolver<LinearSolver>::row_map_.clear();  // TODO smarter solver: only needed if ybus has changed
+    dS_dVm_.resize(0,0);  // TODO smarter solver: only needed if ybus has changed
+    dS_dVa_.resize(0,0);  // TODO smarter solver: only needed if ybus has changed
     // BaseNRSolver<LinearSolver>::dS_dVm_.setZero();  // TODO smarter solver: only needed if ybus has changed
     // BaseNRSolver<LinearSolver>::dS_dVa_.setZero();  // TODO smarter solver: only needed if ybus has changed
     while ((!converged) & (nr_iter_ < max_iter)){
@@ -181,6 +177,7 @@ void BaseNRSolver<LinearSolver>::reset(){
 template<class LinearSolver>
 void BaseNRSolver<LinearSolver>::_dSbus_dV(const Eigen::Ref<const Eigen::SparseMatrix<cplx_type> > & Ybus,
                                            const Eigen::Ref<const CplxVect > & V){
+    // std::cout << "Ybus.nonZeros(): " << Ybus.nonZeros() << std::endl;
     auto timer = CustTimer();
     const auto size_dS = V.size();
     const CplxVect Vnorm = V.array() / V.array().abs();
@@ -198,6 +195,12 @@ void BaseNRSolver<LinearSolver>::_dSbus_dV(const Eigen::Ref<const Eigen::SparseM
 
     cplx_type * ds_dvm_x_ptr = dS_dVm_.valuePtr();
     cplx_type * ds_dva_x_ptr = dS_dVa_.valuePtr();
+
+    // there can be more value in dS_dVm_ or dS_dVa_ than in Ybus, if a line is disconnected for example
+    // TODO smarter solver
+    // const Eigen::Index * col_ptr = dS_dVm_.outerIndexPtr();
+    // const Eigen::Index * row_ptr = dS_dVm_.innerIndexPtr();
+
     unsigned int pos_el = 0;
     for (int col_id=0; col_id < size_dS; ++col_id){
         for (Eigen::Ref<const Eigen::SparseMatrix<cplx_type> >::InnerIterator it(Ybus, col_id); it; ++it)
@@ -298,14 +301,14 @@ void BaseNRSolver<LinearSolver>::_get_values_J(int & nb_obj_this_col,
 
 template<class LinearSolver>
 void BaseNRSolver<LinearSolver>::fill_jacobian_matrix(const Eigen::SparseMatrix<cplx_type> & Ybus,
-                                        const CplxVect & V,
-                                        Eigen::Index slack_bus_id,
-                                        const RealVect & slack_weights,
-                                        const Eigen::VectorXi & pq,
-                                        const Eigen::VectorXi & pvpq,
-                                        const std::vector<int> & pq_inv,
-                                        const std::vector<int> & pvpq_inv
-                                        )
+                                                      const CplxVect & V,
+                                                      Eigen::Index slack_bus_id,
+                                                      const RealVect & slack_weights,
+                                                      const Eigen::VectorXi & pq,
+                                                      const Eigen::VectorXi & pvpq,
+                                                      const std::vector<int> & pq_inv,
+                                                      const std::vector<int> & pvpq_inv
+                                                      )
 {
     /**
     Remember, J has the shape:
@@ -346,7 +349,7 @@ void BaseNRSolver<LinearSolver>::fill_jacobian_matrix(const Eigen::SparseMatrix<
         #endif  // __COUT_TIMES
         // first time i initialized the matrix, so i need to compute its sparsity pattern
         fill_jacobian_matrix_unkown_sparsity_pattern(Ybus, V, slack_bus_id, slack_weights, pq, pvpq, pq_inv, pvpq_inv);
-        fill_value_map(slack_bus_id, pq, pvpq);
+        fill_value_map(slack_bus_id, pq, pvpq, false);
         #ifdef __COUT_TIMES
             std::cout << "\t\t fill_jacobian_matrix_unkown_sparsity_pattern : " << timer2.duration() << std::endl;
         #endif  // __COUT_TIMES
@@ -356,7 +359,7 @@ void BaseNRSolver<LinearSolver>::fill_jacobian_matrix(const Eigen::SparseMatrix<
         #ifdef __COUT_TIMES
             auto timer3 = CustTimer();
         #endif  // 
-        if (BaseNRSolver<LinearSolver>::value_map_.size() == 0) fill_value_map(slack_bus_id,pq, pvpq);
+        if (BaseNRSolver<LinearSolver>::value_map_.size() == 0) fill_value_map(slack_bus_id, pq, pvpq, true);
         fill_jacobian_matrix_kown_sparsity_pattern(slack_bus_id,
                                                    pq, pvpq
                                                    );
@@ -541,11 +544,15 @@ template<class LinearSolver>
 void BaseNRSolver<LinearSolver>::fill_value_map(
         Eigen::Index slack_bus_id,
         const Eigen::VectorXi & pq,
-        const Eigen::VectorXi & pvpq
+        const Eigen::VectorXi & pvpq,
+        bool reset_J
         )
 {
     const int n_pvpq = static_cast<int>(pvpq.size());
-    value_map_ = std::vector<cplx_type*> (J_.nonZeros());
+    value_map_ = std::vector<cplx_type*> ();
+    value_map_.reserve(BaseNRSolver<LinearSolver>::J_.nonZeros());
+    // col_map_ = std::vector<int> (J_.nonZeros());
+    // row_map_ = std::vector<int> (J_.nonZeros());
 
     const auto n_row = J_.cols();
     unsigned int pos_el = 0;
@@ -554,16 +561,18 @@ void BaseNRSolver<LinearSolver>::fill_value_map(
         {
             auto row_id = it.row();
             const auto col_id = it.col() - 1;  // it's equal to "col_"
+            if(reset_J) it.valueRef() = 0.; // "forget" previous J value in this setting
+
             if(row_id==0){
                 // this is the row of the slack bus
                 const Eigen::Index row_id_dS_dVx_r = slack_bus_id;  // same for both matrices
                 if(col_id < n_pvpq){
                     const int col_id_dS_dVa_r = pvpq[col_id];
-                    value_map_[pos_el] = &dS_dVa_.coeffRef(row_id_dS_dVx_r, col_id_dS_dVa_r);
+                    value_map_.push_back(&dS_dVa_.coeffRef(row_id_dS_dVx_r, col_id_dS_dVa_r));
                 }
                 else{
                     const int col_id_dS_dVm_r = pq[col_id - n_pvpq];
-                    value_map_[pos_el] = &dS_dVm_.coeffRef(row_id_dS_dVx_r, col_id_dS_dVm_r);
+                    value_map_.push_back(&dS_dVm_.coeffRef(row_id_dS_dVx_r, col_id_dS_dVm_r));
                 }
             }else{
                 row_id -= 1;  // "do not consider" the row for slack bus (handled above)
@@ -573,7 +582,7 @@ void BaseNRSolver<LinearSolver>::fill_value_map(
                     const int row_id_dS_dVa_r = pvpq[row_id];
                     const int col_id_dS_dVa_r = pvpq[col_id];
                     // this_el = dS_dVa_r.coeff(row_id_dS_dVa_r, col_id_dS_dVa_r);
-                    value_map_[pos_el] = &dS_dVa_.coeffRef(row_id_dS_dVa_r, col_id_dS_dVa_r);
+                    value_map_.push_back(&dS_dVa_.coeffRef(row_id_dS_dVa_r, col_id_dS_dVa_r));
 
                     // I don't need to perform these checks: if they failed, the element would not be in J_ in the first place
                     // const int is_row_non_null = pq_inv[row_id_dS_dVa_r];
@@ -587,25 +596,27 @@ void BaseNRSolver<LinearSolver>::fill_value_map(
                     const int row_id_dS_dVa_i = pq[row_id - n_pvpq];
                     const int col_id_dS_dVa_i = pvpq[col_id];
                     // this_el = dS_dVa_i.coeff(row_id_dS_dVa_i, col_id_dS_dVa_i);
-                    value_map_[pos_el] = &dS_dVa_.coeffRef(row_id_dS_dVa_i, col_id_dS_dVa_i);
+                    value_map_.push_back(&dS_dVa_.coeffRef(row_id_dS_dVa_i, col_id_dS_dVa_i));
                 }else if((col_id >= n_pvpq) && (row_id < n_pvpq)){
                     // this is the J12 part (dS_dVm_r)
                     const int row_id_dS_dVm_r = pvpq[row_id];
                     const int col_id_dS_dVm_r = pq[col_id - n_pvpq];
                     // this_el = dS_dVm_r.coeff(row_id_dS_dVm_r, col_id_dS_dVm_r);
-                    value_map_[pos_el] = &dS_dVm_.coeffRef(row_id_dS_dVm_r, col_id_dS_dVm_r);
+                    value_map_.push_back(&dS_dVm_.coeffRef(row_id_dS_dVm_r, col_id_dS_dVm_r));
                 }else if((col_id >= n_pvpq) && (row_id >= n_pvpq)){
                     // this is the J22 part (dS_dVm_i)
                     const int row_id_dS_dVm_i = pq[row_id - n_pvpq];
                     const int col_id_dS_dVm_i = pq[col_id - n_pvpq];
                     // this_el = dS_dVm_i.coeff(row_id_dS_dVm_i, col_id_dS_dVm_i);
-                    value_map_[pos_el] = &dS_dVm_.coeffRef(row_id_dS_dVm_i, col_id_dS_dVm_i);
+                    value_map_.push_back(&dS_dVm_.coeffRef(row_id_dS_dVm_i, col_id_dS_dVm_i));
                 }
             }
             // go to the next element
             ++pos_el;
         }
     }
+    dS_dVa_.makeCompressed();
+    dS_dVm_.makeCompressed();
 }
 
 template<class LinearSolver>

src/BaseNRSolverSingleSlack.h

@@ -58,7 +58,8 @@ class BaseNRSolverSingleSlack : public BaseNRSolver<LinearSolver>
                  );
 
         void fill_value_map(const Eigen::VectorXi & pq,
-                            const Eigen::VectorXi & pvpq);
+                            const Eigen::VectorXi & pvpq,
+                            bool reset_J);
 
 };