Added Jacobian Assembly. Fixed and Added New Functionality to COO Mat…

…rices

+Matrix Assembly
+Refactored COO Matrix
+COO can return CSR Data format

ComponentLib/PowerElectronicsComponents/Capacitor/Capacitor.cpp

@@ -91,7 +91,7 @@ int Capacitor<ScalarT, IdxT>::evaluateJacobian()
     COO_Matrix<ScalarT,IdxT> Jacder = COO_Matrix<ScalarT, IdxT>(rcordder, ccordder, valsder,3,3);
 
     //Perform dF/dy + \alpha dF/dy'
-    this->J_.AXPY(this->alpha_, &Jacder);
+    this->J_.AXPY(this->alpha_, Jacder);
 
     return 0;
 }

ComponentLib/PowerElectronicsComponents/DiscreteGenerator/DiscreteGenerator.cpp

@@ -237,7 +237,7 @@ int DiscreteGenerator<ScalarT, IdxT>::evaluateJacobian()
 
     //Perform dF/dy + \alpha dF/dy'
 
-    this->J_.AXPY(-this->alpha_, &Jacder);
+    this->J_.AXPY(-this->alpha_, Jacder);
 
     return 0;
 }

ComponentLib/PowerElectronicsComponents/Inductor/Inductor.cpp

@@ -93,7 +93,7 @@ int Inductor<ScalarT, IdxT>::evaluateJacobian()
     COO_Matrix<ScalarT,IdxT> Jacder = COO_Matrix<ScalarT, IdxT>(rcordder, ccordder, valsder,3,3);
 
     //Perform dF/dy + \alpha dF/dy'
-    this->J_.AXPY(this->alpha_, &Jacder);
+    this->J_.AXPY(this->alpha_, Jacder);
 
     return 0;
 }

Examples/RLCircuit/RLCircuit.cpp

@@ -80,9 +80,12 @@ int main(int argc, char const *argv[])
 	}
 	std::cout << "}\n";
 
-	induct->updateTime(0.0, 1.0);
-	induct->evaluateJacobian();
-	induct->getJacobian().printMatrix(true);
+
+	sysmodel->updateTime(0.0, 1.0);
+	sysmodel->evaluateJacobian();
+	sysmodel->getJacobian().printMatrix();
+
+
 
 
 	return 0;

Examples/SparseTest/SparseTest.cpp

@@ -10,9 +10,6 @@
 
 #include <SparseMatrix/COO_Matrix.hpp>
 
-
-
-
 int main(int argc, char const *argv[])
 {
 	std::vector<double> val{0.1, 0.2, 0.3, 0.4};
@@ -35,21 +32,55 @@ int main(int argc, char const *argv[])
 	}
 
 	std::cout << "A:\n"; 
-	A.printMatrix(true);
+	A.printMatrix();
 
 	std::vector<double> val2{0.5, 0.6, 0.7, 0.8, 1.0};
 	std::vector<size_t> x2{0,2,0,2,1};
 	std::vector<size_t> y2{3,3,2,2,3};
 	COO_Matrix<double, size_t> B = COO_Matrix<double, size_t>(x2,y2,val2,m,n);
 
 	std::cout << "B:\n";
-	B.printMatrix(true);
+	B.printMatrix();
 
-	A.AXPY(2.0, &B);
+	A.AXPY(2.0, B);
 
 	std::cout << "A + 2B:\n";
-	A.printMatrix(true);
-
+	A.printMatrix();
+
+	std::vector<size_t> r;
+	std::vector<size_t> c;
+	std::vector<double> v;
+	std::tie(r,c,v) = A.getDataToCSR();
+
+	for (size_t i = 0; i < r.size() - 1; i++)
+	{
+		std::cout << r[i] << std::endl;
+		size_t rdiff = r[i+1] - r[i];
+		for (size_t j = 0; j < rdiff; j++)
+		{
+			std::cout << c[j + r[i]] << ", " << v[j + r[i]] << std::endl;
+		}
+	}
+	std::cout << r[r.size()-1] << std::endl;
+
+	//Basic Verification test
+	std::vector<size_t> rtest = {0, 2, 4, 7, 8};
+	std::vector<size_t> ctest = {2,3,2,3,1,2,3,2};
+	std::vector<double> valtest = {1.4, 1.0, 0.4, 2.2, 0.1, 1.6, 1.2, 0.3};
+
+	assert(rtest.size() == r.size());
+	assert(ctest.size() == c.size());
+	assert(valtest.size() == v.size());
+
+	int failval = 0;
+	for (size_t i = 0; i < rtest.size(); i++) if (r[i] != rtest[i]) failval--; 
+	for (size_t i = 0; i < ctest.size(); i++)
+	{
+		double vdiff = v[i] - valtest[i];
+		if (c[i] != ctest[i] || -1e-14 > vdiff  || vdiff > 1e-14) failval--;
+	}
 
-	return 0;
+	std::cout << failval << std::endl;
+
+	return failval;
 }

PowerElectronicsModel.hpp

@@ -23,8 +23,8 @@ class PowerElectronicsModel : public ModelEvaluatorImpl<ScalarT, IdxT>
     // using ModelEvaluatorImpl<ScalarT, IdxT>::size_quad_;
     // using ModelEvaluatorImpl<ScalarT, IdxT>::size_opt_;
     using ModelEvaluatorImpl<ScalarT, IdxT>::nnz_;
-    // using ModelEvaluatorImpl<ScalarT, IdxT>::time_;
-    // using ModelEvaluatorImpl<ScalarT, IdxT>::alpha_;
+    using ModelEvaluatorImpl<ScalarT, IdxT>::time_;
+    using ModelEvaluatorImpl<ScalarT, IdxT>::alpha_;
     using ModelEvaluatorImpl<ScalarT, IdxT>::y_;
     using ModelEvaluatorImpl<ScalarT, IdxT>::yp_;
     // using ModelEvaluatorImpl<ScalarT, IdxT>::yB_;
@@ -98,6 +98,11 @@ class PowerElectronicsModel : public ModelEvaluatorImpl<ScalarT, IdxT>
         return 0;
     }
 
+    /**
+     * @brief Set intial y and y' of each component
+     * 
+     * @return int 
+     */
     int initialize()
     {
 
@@ -106,7 +111,18 @@ class PowerElectronicsModel : public ModelEvaluatorImpl<ScalarT, IdxT>
         {
             component->initialize();
         }
+        this->distributeVectors();
+
+        return 0;
+    }
 
+    /**
+     * @brief Distribute y and y' to each component based of node connection graph
+     * 
+     * @return int 
+     */
+    int distributeVectors()
+    {
         for(const auto& component : components_)
         {
             for(IdxT j=0; j<component->size(); ++j)
@@ -123,23 +139,19 @@ class PowerElectronicsModel : public ModelEvaluatorImpl<ScalarT, IdxT>
         return 0;
     }
 
+    /**
+     * @brief Evaluate Residuals at each component then collect them
+     * 
+     * @return int 
+     */
     int evaluateResidual()
     {
         for (IdxT i = 0; i < this->f_.size(); i++)
         {
             f_[i] = 0;
         }
 
-        // Update variables
-        for(const auto& component : components_)
-        {
-            for(IdxT j=0; j<component->size(); ++j)
-            {
-                component->y()[j]  = y_[component->getNodeConnection(j)];
-				component->yp()[j] = yp_[component->getNodeConnection(j)];
-            }
-            component->evaluateResidual();
-        }
+        this->distributeVectors();
 
         // Update system residual vector
 
@@ -161,6 +173,25 @@ class PowerElectronicsModel : public ModelEvaluatorImpl<ScalarT, IdxT>
      */
     int evaluateJacobian()
 	{
+        this->J_.zeroMatrix();
+        this->distributeVectors();
+
+        //Evaluate component jacs
+        for(const auto& component : components_)
+        {
+            component->evaluateJacobian();
+            std::vector<IdxT> r;
+            std::vector<IdxT> c;
+            std::vector<ScalarT> v;
+	        std::tie(r, c, v) = component->getJacobian().getEntrieCopies();
+            for (IdxT i = 0; i < static_cast<IdxT>(r.size()); i++)
+            {
+                r[i] = component->getNodeConnection(r[i]);
+                c[i] = component->getNodeConnection(c[i]);
+            }
+            this->J_.AXPY(1.0, r, c, v);
+        }
+
 		return 0;
 	}
 
@@ -205,8 +236,18 @@ class PowerElectronicsModel : public ModelEvaluatorImpl<ScalarT, IdxT>
         return 0;
     }
 
+    /**
+     * @brief Distribute time and time scaling for each component
+     * 
+     * @param t 
+     * @param a 
+     */
     void updateTime(ScalarT t, ScalarT a)
     {
+        for(const auto& component : components_)
+        {
+            component->updateTime(t, a);
+        }
     }
 
     void addComponent(component_type* component)