FEAT: Allow solver_conduction to use sources w linear dependence on value

- especially useful for e- temp, where there are sources with & without dependence on temperature

> This does not yet address the BC's. But that is doable now that we know how to use optional arguments w functions.

Author: abukowski21

include/solvers.h

@@ -40,7 +40,8 @@ arma_vec solver_conduction(
 			   arma_vec dx,
 			   precision_t dt,
 			   int64_t nGCs,
-			   bool return_diff);
+			   bool return_diff = false,
+			   arma_vec source2 = arma_vec());
 
 arma_cube solver_chemistry(arma_cube density,
 			   arma_cube source,

src/solver_conduction.cpp

@@ -25,7 +25,9 @@ arma_vec solver_conduction(arma_vec value,
                            arma_vec dx,
                            precision_t dt,
                            int64_t nGCs,
-                           bool return_diff) {
+                           bool return_diff, // (optional) False by default (return new `value`)
+                           arma_vec source2 // (optional) Sources dependent on `value`
+                           ) {
 
   int64_t nPts = value.n_elem;
 
@@ -47,14 +49,20 @@ arma_vec solver_conduction(arma_vec value,
   arma_vec conduction(nPts);
   conduction.zeros();
 
+  // If source2 is not given, set it to zero:
+  if (source2.n_elem == 0){
+    source2.set_size(source.n_elem);
+    source2.zeros();
+  }
+
   int64_t i;
 
   for (i = nGCs; i < nPts - nGCs; i++)
     dl(i) = di(i + 1) - di(i - 1) * r(i) * r(i) - di(i) * (1.0 - r(i) * r(i));
 
   arma_vec a = di / du22 % r - dl / du12 % r % r;
   arma_vec c = di / du22 + dl / du12;
-  arma_vec b = -1.0 / m - di / du22 % (1.0 + r) - dl / du12 % (1.0 - r % r);
+  arma_vec b = -1.0 / m - di / du22 % (1.0 + r) - dl / du12 % (1.0 - r % r) + source2 % front * dt;
   arma_vec d = -1.0 * (value / m + source % front * dt);
 
   // Lower BCs (fixed value):