Merge branch 'flux_jacobian_matrix' into develop

Adds full matrix flux jacobians
Allows main diagonal ability to hold full matrices instead of just scalars

fluxJacobian.hpp

@@ -24,6 +24,7 @@
 #include <memory>          // unique_ptr
 #include "vector3d.hpp"
 #include "uncoupledScalar.hpp"
+#include "matrix.hpp"
 
 using std::vector;
 using std::string;
@@ -40,26 +41,21 @@ class sutherland;
 class turbModel;
 class input;
 class genArray;
-class squareMatrix;
 
 // This class holds the flux jacobians for the flow and turbulence equations.
 // In the LU-SGS method the jacobians are scalars.
 
 class fluxJacobian {
-  double flowJacobian_;
-  double turbJacobian_;
+  squareMatrix flowJacobian_;
+  squareMatrix turbJacobian_;
 
  public:
   // constructors
-  fluxJacobian(const double &flow, const double &turb) :
-      flowJacobian_(flow), turbJacobian_(turb) {}
+  fluxJacobian(const double &flow, const double &turb);
+  fluxJacobian(const int &flowSize, const int &turbSize);
   fluxJacobian() : fluxJacobian(0.0, 0.0) {}
   explicit fluxJacobian(const uncoupledScalar &specRad) :
       fluxJacobian(specRad.FlowVariable(), specRad.TurbVariable()) {}
-  fluxJacobian(const primVars &, const unitVec3dMag<double> &,
-               const unitVec3dMag<double> &, const idealGas &,
-               const sutherland &, const double &,
-               const unique_ptr<turbModel> &, const input &, const bool &);
 
   // move constructor and assignment operator
   fluxJacobian(fluxJacobian&&) noexcept = default;
@@ -70,27 +66,38 @@ class fluxJacobian {
   fluxJacobian& operator=(const fluxJacobian&) = default;
 
   // member functions
-  double FlowJacobian() const { return flowJacobian_; }
-  double TurbulenceJacobian() const { return turbJacobian_; }
-
-  void AddInviscidJacobian(const primVars &, const unitVec3dMag<double> &,
-                           const unitVec3dMag<double> &, const idealGas &,
-                           const unique_ptr<turbModel> &, const bool &);
-  void AddViscousJacobian(const primVars &, const unitVec3dMag<double> &,
-                          const unitVec3dMag<double> &, const idealGas &,
-                          const sutherland &, const double &,
-                          const unique_ptr<turbModel> &, const bool &);
-  void AddTurbSourceJacobian(const primVars &, const sutherland &,
-                             const double &, const unique_ptr<turbModel> &);
+  squareMatrix FlowJacobian() const { return flowJacobian_; }
+  squareMatrix TurbulenceJacobian() const { return turbJacobian_; }
+
+  void AddToFlowJacobian(const squareMatrix &jac) {flowJacobian_ += jac;}
+  void AddToTurbJacobian(const squareMatrix &jac) {turbJacobian_ += jac;}
+  void SubtractFromFlowJacobian(const squareMatrix &jac) {flowJacobian_ -= jac;}
+  void SubtractFromTurbJacobian(const squareMatrix &jac) {turbJacobian_ -= jac;}
+
+  void RusanovFluxJacobian(const primVars &, const primVars &,
+                           const idealGas &, const vector3d<double> &,
+                           const bool &, const input &);
+  void InvFluxJacobian(const primVars &, const idealGas &,
+                       const vector3d<double> &, const input &);
+  void ApproxRoeFluxJacobian(const primVars &, const primVars &,
+                             const idealGas &, const vector3d<double> &,
+                             const bool &, const input &);
+  void DelPrimativeDelConservative(const primVars &, const idealGas &,
+                                   const input &);
+
+  void ApproxTSLJacobian(const primVars &, const idealGas &,
+                         const sutherland &,
+                         const vector3d<double> &, const double &,
+                         const unique_ptr<turbModel> &, const input &);
 
   void Zero() {
-    flowJacobian_ = 0.0;
-    turbJacobian_ = 0.0;
+    flowJacobian_.Zero();
+    turbJacobian_.Zero();
   }
 
   genArray ArrayMult(genArray) const;
-
-  fluxJacobian Inverse(const bool &) const;
+  bool IsScalar() const;
+  void Inverse(const bool &);
 
   inline fluxJacobian & operator+=(const fluxJacobian &);
   inline fluxJacobian & operator-=(const fluxJacobian &);
@@ -225,14 +232,6 @@ inline const fluxJacobian operator/(const double &lhs, fluxJacobian rhs) {
 
 ostream &operator<<(ostream &os, const fluxJacobian &);
 
-squareMatrix RusanovFluxJacobian(const primVars &, const primVars &,
-                                 const idealGas &, const vector3d<double> &,
-                                 const bool &);
-squareMatrix InvFluxJacobian(const primVars &, const idealGas &,
-                             const vector3d<double> &);
-squareMatrix ApproxRoeFluxJacobian(const primVars &, const primVars &,
-                                   const idealGas &, const vector3d<double> &,
-                                   const bool &);
 genArray RusanovOffDiagonal(const primVars &, const genArray &,
                             const unitVec3dMag<double> &,
                             const unitVec3dMag<double> &,
@@ -247,11 +246,5 @@ genArray RoeOffDiagonal(const primVars &, const primVars &, const genArray &,
                         const unique_ptr<turbModel> &, const input &,
                         const bool &);
 
-squareMatrix DelPrimativeDelConservative(const primVars &, const idealGas &);
-
-squareMatrix ApproxTSLJacobian(const primVars &, const idealGas &,
-                               const sutherland &,
-                               const vector3d<double> &, const double &,
-                               const unique_ptr<turbModel> &);
 
 #endif

input.cpp

@@ -549,14 +549,23 @@ void input::CalcCFL(const int &i) {
   }
 }
 
+// member function to determine number of turbulence equations
+int input::NumTurbEquations() const {
+  auto numEqns = 0;
+  if (this->IsTurbulent()) {
+    numEqns = 2;
+  }
+  return numEqns;
+}
+
 // member function to determine number of equations to solver for
 int input::NumEquations() const {
   auto numEqns = 0;
   if ((equationSet_ == "euler") ||
       (equationSet_ == "navierStokes")) {
-    numEqns = 5;
+    numEqns = this->NumFlowEquations();
   } else if (equationSet_ == "rans") {
-    numEqns = 7;
+    numEqns = this->NumFlowEquations() + this->NumTurbEquations();
   } else {
     cerr << "ERROR: Equations set is not recognized. Cannot determine number "
             "of equations!" << endl;

input.hpp

@@ -23,6 +23,7 @@
 #include <memory>                   // unique_ptr
 #include "vector3d.hpp"
 #include "boundaryConditions.hpp"
+#include "macros.hpp"
 
 using std::vector;
 using std::string;
@@ -161,14 +162,16 @@ class input {
 
   int NumVars() const {return vars_.size();}
   int NumEquations() const;
+  int NumFlowEquations() const {return NUMFLOWVARS;}
+  int NumTurbEquations() const;
 
   void ReadInput(const int &);
 
   bool IsImplicit() const;
   bool IsViscous() const;
   bool IsTurbulent() const;
   bool IsBlockMatrix() const;
-  
+
   string OrderOfAccuracy() const;
 
   double FarfieldTurbIntensity() const {return farfieldTurbInten_;}

kdtree.cpp

@@ -15,8 +15,7 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>. */
 
 #include <iostream>     // cout
-#include <cstdlib>      // exit()
-#include <algorithm>    // sort
+#include <algorithm>    // nth_element
 #include <limits>       // numeric_limits
 #include <vector>       // vector
 #include "kdtree.hpp"
@@ -88,36 +87,16 @@ void kdtree::BuildKdtree(const int &start, const int &end,
   // top portion of vector, and all other elements in bottom
   // portion. This is necessary because median calculation does
   // not sort the data
-
   // start is at median, so use next index down
-  auto pTop = start + 1;
-  // end is just past the last value, so use index up
-  auto pBot = end - 1;
-
-  while (pTop <= pBot) {
-    // if inside this conditional, point found belongs in top
-    // portion of vector, so just increment index and move on
-    if (nodes_[pTop][dim] <= median) {
-      pTop++;
-    } else {
-      // if inside this conditional, point from top index
-      // belongs on bottom
-
-      // while loop is used to find a point from bottom index
-      // that belongs on top
-      while (nodes_[pBot][dim] > median && pTop < pBot) {
-        pBot--;
-      }
-
-      // swap points so point from top belonging on bottom,
-      // and point from bottom belonging on top are in correct place
-      std::swap(nodes_[pTop], nodes_[pBot]);
-      pBot--;
-    }
-  }
+  const auto splitIndex = std::partition(nodes_.begin() + start + 1,
+                                         nodes_.begin() + end,
+                                         [&median, &dim]
+                                         (const vector3d<double> &p1)
+                                         {return p1[dim] <= median;});
 
   // record split index as index of right branch
-  right_[start] = pTop;
+  // std::partition returns iterator to 2nd partition which is right side
+  right_[start] = splitIndex - nodes_.begin();
 
   // recursively build left and right branches
   this->BuildKdtree(start + 1, right_[start], depth + 1);  // left

macros.hpp

@@ -20,12 +20,11 @@
 
 #define NUMVARS 7
 #define NUMFLOWVARS 5
-#define NUMTURBVARS 2
 #define EPS 1.0e-30
 #define ROOTP 0
 #define DEFAULTWALLDIST 1.0e10
 #define MAJORVERSION 0
-#define MINORVERSION 1
+#define MINORVERSION 2
 #define PATCHNUMBER 0
 
 #endif

matrix.cpp

@@ -173,13 +173,13 @@ void squareMatrix::Identity() {
 }
 
 // member function to do matrix/vector multplication
-genArray squareMatrix::VecMult(const genArray &vec) const {
+genArray squareMatrix::ArrayMult(const genArray &vec, const int pos) const {
   // vec -- vector to multiply with
-  
+
   genArray product(0.0);
   for (auto rr = 0; rr < size_; rr++) {
     for (auto cc = 0; cc < size_; cc++) {
-      product[rr] += (*this)(rr, cc) * vec[cc];
+      product[rr] += (*this)(rr, cc) * vec[pos + cc];
     }
   }
   return product;

matrix.hpp

@@ -61,8 +61,8 @@ class squareMatrix {
   void Zero();
   void Identity();
   squareMatrix MatMult(const squareMatrix &) const;
-  genArray VecMult(const genArray &) const;
-  
+  genArray ArrayMult(const genArray &, const int = 0) const;
+
   // operator overloads
   double & operator()(const int &r, const int &c) {
     return data_[this->GetLoc(r, c)];