From 1f3019c0f402c99c2c842c3a3cc7464dd2d5749b Mon Sep 17 00:00:00 2001 From: Michael Date: Sun, 24 Apr 2016 21:34:06 -0700 Subject: [PATCH 1/5] Changed fluxJacobian class to use squareMatrix instead of double to allow for block matrix methods --- fluxJacobian.cpp | 58 ++++++++++++++++++++++++++++++++---------------- fluxJacobian.hpp | 22 +++++++++--------- matrix.cpp | 6 ++--- matrix.hpp | 4 ++-- procBlock.cpp | 27 +++++++++++----------- 5 files changed, 68 insertions(+), 49 deletions(-) diff --git a/fluxJacobian.cpp b/fluxJacobian.cpp index 3eccb19..96933fb 100644 --- a/fluxJacobian.cpp +++ b/fluxJacobian.cpp @@ -23,7 +23,6 @@ #include "input.hpp" // input #include "primVars.hpp" // primVars #include "genArray.hpp" // genArray -#include "matrix.hpp" // squareMatrix #include "inviscidFlux.hpp" // ConvectiveFluxUpdate using std::cout; @@ -34,6 +33,21 @@ using std::string; using std:: unique_ptr; // constructor +// if constructed with two doubles, create scalar squareMatrix +fluxJacobian::fluxJacobian(const double &flow, const double &turb) { + flowJacobian_ = squareMatrix(1); + flowJacobian_ += flow; + + turbJacobian_ = squareMatrix(1); + turbJacobian_ += turb; +} + +// if constructed with two intss, create scalar squareMatrix with given size +fluxJacobian::fluxJacobian(const int &flowSize, const int &turbSize) { + flowJacobian_ = squareMatrix(flowSize); + turbJacobian_ = squareMatrix(turbSize); +} + fluxJacobian::fluxJacobian(const primVars &state, const unitVec3dMag &fAreaL, const unitVec3dMag &fAreaR, @@ -53,8 +67,8 @@ fluxJacobian::fluxJacobian(const primVars &state, // implicit matrix // initialize jacobians - flowJacobian_ = 0.0; - turbJacobian_ = 0.0; + flowJacobian_ = squareMatrix(1); + turbJacobian_ = squareMatrix(1); const auto isTurbulent = inp.IsTurbulent(); @@ -140,27 +154,33 @@ void fluxJacobian::AddTurbSourceJacobian(const primVars &state, // member function to multiply the flux jacobians with a genArray genArray fluxJacobian::ArrayMult(genArray arr) const { - arr[0] *= flowJacobian_; - arr[1] *= flowJacobian_; - arr[2] *= flowJacobian_; - arr[3] *= flowJacobian_; - arr[4] *= flowJacobian_; - - arr[5] *= turbJacobian_; - arr[6] *= turbJacobian_; - + if (this->IsScalar()) { + arr[0] *= flowJacobian_(0, 0); + arr[1] *= flowJacobian_(0, 0); + arr[2] *= flowJacobian_(0, 0); + arr[3] *= flowJacobian_(0, 0); + arr[4] *= flowJacobian_(0, 0); + + arr[5] *= turbJacobian_(0, 0); + arr[6] *= turbJacobian_(0, 0); + } else { + arr = flowJacobian_.ArrayMult(arr); + arr = turbJacobian_.ArrayMult(arr, NUMFLOWVARS); + } return arr; } +bool fluxJacobian::IsScalar() const { + return (flowJacobian_.Size() > 1) ? false : true; +} + // function to take the inverse of a flux jacobian -fluxJacobian fluxJacobian::Inverse(const bool &isTurbulent) const { - auto inv = *this; - inv.flowJacobian_ = 1.0 / inv.flowJacobian_; +void fluxJacobian::Inverse(const bool &isTurbulent) { + flowJacobian_.Inverse(); if (isTurbulent) { - inv.turbJacobian_ = 1.0 / inv.turbJacobian_; + turbJacobian_.Inverse(); } - return inv; } // non-member functions @@ -408,8 +428,8 @@ genArray RoeOffDiagonal(const primVars &left, const primVars &right, 0.5 * specRad.ArrayMult(update) : fAreaL.Mag() * ((newFlux - oldFlux).ConvertToGenArray()) - 0.5 * specRad.ArrayMult(update); - // fAreaL.Mag() * invFluxJac.VecMult(update) + 0.5 * specRadArr * update : - // fAreaL.Mag() * invFluxJac.VecMult(update) - 0.5 * specRadArr * update; + // fAreaL.Mag() * invFluxJac.ArrayMult(update) + 0.5 * specRadArr * update : + // fAreaL.Mag() * invFluxJac.ArrayMult(update) - 0.5 * specRadArr * update; } // change of variable matrix going frim primative to conservative variables diff --git a/fluxJacobian.hpp b/fluxJacobian.hpp index 6e0a0f1..18bc689 100644 --- a/fluxJacobian.hpp +++ b/fluxJacobian.hpp @@ -24,6 +24,7 @@ #include // unique_ptr #include "vector3d.hpp" #include "uncoupledScalar.hpp" +#include "matrix.hpp" using std::vector; using std::string; @@ -40,19 +41,18 @@ 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()) {} @@ -70,8 +70,8 @@ class fluxJacobian { fluxJacobian& operator=(const fluxJacobian&) = default; // member functions - double FlowJacobian() const { return flowJacobian_; } - double TurbulenceJacobian() const { return turbJacobian_; } + squareMatrix FlowJacobian() const { return flowJacobian_; } + squareMatrix TurbulenceJacobian() const { return turbJacobian_; } void AddInviscidJacobian(const primVars &, const unitVec3dMag &, const unitVec3dMag &, const idealGas &, @@ -84,13 +84,13 @@ class fluxJacobian { const double &, const unique_ptr &); 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 &); diff --git a/matrix.cpp b/matrix.cpp index 6e097ac..c859601 100644 --- a/matrix.cpp +++ b/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; diff --git a/matrix.hpp b/matrix.hpp index 481c4ab..3b8fcf6 100644 --- a/matrix.hpp +++ b/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)]; diff --git a/procBlock.cpp b/procBlock.cpp index 2f39386..e7b1738 100644 --- a/procBlock.cpp +++ b/procBlock.cpp @@ -932,8 +932,7 @@ void procBlock::InvertDiagonal(multiArray3d &mainDiagonal, mainDiagonal(ip, jp, kp) += diagVolTime; // invert main diagonal - mainDiagonal(ip, jp, kp) = - mainDiagonal(ip, jp, kp).Inverse(inp.IsTurbulent()); + mainDiagonal(ip, jp, kp).Inverse(inp.IsTurbulent()); } } } @@ -1067,14 +1066,14 @@ For viscous simulations, the viscous contribution to the spectral radius K is used, and everything else remains the same. */ void procBlock::LUSGS_Forward(const vector> &reorder, - multiArray3d &x, - const multiArray3d &solTimeMmN, - const multiArray3d &solDeltaNm1, - const idealGas &eqnState, const input &inp, - const sutherland &suth, - const unique_ptr &turb, - const multiArray3d &aInv, - const int &sweep) const { + multiArray3d &x, + const multiArray3d &solTimeMmN, + const multiArray3d &solDeltaNm1, + const idealGas &eqnState, const input &inp, + const sutherland &suth, + const unique_ptr &turb, + const multiArray3d &aInv, + const int &sweep) const { // reorder -- order of cells to visit (this should be ordered in hyperplanes) // x -- correction - added to solution at time n to get to time n+1 (assumed // to be zero to start) @@ -1191,10 +1190,10 @@ void procBlock::LUSGS_Forward(const vector> &reorder, // normal at lower boundaries needs to be reversed, so add instead // of subtract L x(ig, jg, kg) = aInv(ip, jp, kp).ArrayMult(-thetaInv * - residual_(ip, jp, kp) - - solDeltaNm1(ip, jp, kp) - - solTimeMmN(ip, jp, kp) + - L - U); + residual_(ip, jp, kp) - + solDeltaNm1(ip, jp, kp) - + solTimeMmN(ip, jp, kp) + + L - U); } // end forward sweep } From 9d87048fa5b94f3cde13d9e1caa0ae6317cfa53e Mon Sep 17 00:00:00 2001 From: Michael Date: Mon, 25 Apr 2016 23:30:36 -0700 Subject: [PATCH 2/5] Changing flux jacobian functions to work with squareMatrix --- fluxJacobian.cpp | 461 +++++++++++++++++++++-------------------------- fluxJacobian.hpp | 42 ++--- input.cpp | 13 +- input.hpp | 5 +- macros.hpp | 3 +- procBlock.cpp | 51 +++--- 6 files changed, 261 insertions(+), 314 deletions(-) diff --git a/fluxJacobian.cpp b/fluxJacobian.cpp index 96933fb..9fdb19e 100644 --- a/fluxJacobian.cpp +++ b/fluxJacobian.cpp @@ -48,110 +48,8 @@ fluxJacobian::fluxJacobian(const int &flowSize, const int &turbSize) { turbJacobian_ = squareMatrix(turbSize); } -fluxJacobian::fluxJacobian(const primVars &state, - const unitVec3dMag &fAreaL, - const unitVec3dMag &fAreaR, - const idealGas &eos, const sutherland &suth, - const double &vol, - const unique_ptr &turb, - const input &inp, const bool &mainDiagonal) { - // state -- primative variables - // fAreaL -- face area for left face - // fAreaR -- face area for right face - // eos -- equation of state - // suth -- sutherland's law for viscosity - // vol -- cell volume - // turb -- turbulence model - // inp -- input variables - // mainDiagonal -- flag to determine if jacobian is for main diagonal of - // implicit matrix - - // initialize jacobians - flowJacobian_ = squareMatrix(1); - turbJacobian_ = squareMatrix(1); - - const auto isTurbulent = inp.IsTurbulent(); - - this->AddInviscidJacobian(state, fAreaL, fAreaR, eos, turb, isTurbulent); - - // if viscous, add viscous jacobians - if (inp.IsViscous()) { - this->AddViscousJacobian(state, fAreaL, fAreaR, eos, suth, vol, - turb, isTurbulent); - } - - // source term is only added to main diagonal, - // and only for turbulence equations - if (mainDiagonal && isTurbulent) { - this->AddTurbSourceJacobian(state, suth, vol, turb); - } -} // member functions -// member function to add inviscid jacobians -void fluxJacobian::AddInviscidJacobian(const primVars &state, - const unitVec3dMag &fAreaL, - const unitVec3dMag &fAreaR, - const idealGas &eos, - const unique_ptr &turb, - const bool &isTurbulent) { - // state -- primative variables - // fAreaL -- face area for left face - // fAreaR -- face area for right face - // eos -- equation of state - // turb -- turbulence model - // isTurbulent -- flag to determine if simulation is turbulent - - flowJacobian_ += state.InvCellSpectralRadius(fAreaL, fAreaR, eos); - - if (isTurbulent) { - turbJacobian_ += turb->InviscidSpecRad(state, fAreaL, fAreaR); - } -} - -// member function to add viscous jacobians -void fluxJacobian::AddViscousJacobian(const primVars &state, - const unitVec3dMag &fAreaL, - const unitVec3dMag &fAreaR, - const idealGas &eos, - const sutherland &suth, - const double &vol, - const unique_ptr &turb, - const bool &isTurbulent) { - // state -- primative variables - // fAreaL -- face area for left face - // fAreaR -- face area for right face - // eos -- equation of state - // suth -- sutherland's law for viscosity - // vol -- cell volume - // turb -- turbulence model - // isTurbulent -- flag to determine if simulation is turbulent - - // factor of 2 because viscous spectral radius is not halved (Blazek 6.53) - flowJacobian_ += 2.0 * state.ViscCellSpectralRadius(fAreaL, fAreaR, eos, suth, - vol, turb); - - if (isTurbulent) { - // factor of 2 because viscous spectral radius is not halved (Blazek 6.53) - turbJacobian_ += 2.0 * turb->ViscSpecRad(state, fAreaL, fAreaR, eos, suth, - vol); - } -} - -// member function to add source jacobians -// this should only be used for the main diagonal -void fluxJacobian::AddTurbSourceJacobian(const primVars &state, - const sutherland &suth, - const double &vol, - const unique_ptr &turb) { - // state -- primative variables - // suth -- sutherland's law for viscosity - // vol -- cell volume - // turb -- turbulence model - - turbJacobian_ -= turb->SrcSpecRad(state, suth, vol); -} - // member function to multiply the flux jacobians with a genArray genArray fluxJacobian::ArrayMult(genArray arr) const { if (this->IsScalar()) { @@ -183,15 +81,6 @@ void fluxJacobian::Inverse(const bool &isTurbulent) { } } -// non-member functions -// ---------------------------------------------------------------------------- -// operator overload for << - allows use of cout, cerr, etc. -ostream &operator<<(ostream &os, fluxJacobian &jacobian) { - os << jacobian.FlowJacobian() << endl; - os << jacobian.TurbulenceJacobian() << endl; - return os; -} - /* Function to calculate Rusanov flux jacobian. The Rusanov flux is defined as shown below. @@ -206,15 +95,18 @@ jacobians. In the above equations the dissipation term L is held constant during differentiation. A represents the convective flux jacobian matrix. */ -squareMatrix RusanovFluxJacobian(const primVars &left, const primVars &right, - const idealGas &eos, - const vector3d &areaNorm, - const bool &positive) { +void fluxJacobian::RusanovFluxJacobian(const primVars &left, + const primVars &right, + const idealGas &eos, + const vector3d &areaNorm, + const bool &positive, + const input &inp) { // left -- primative variables from left side // right -- primative variables from right side // eos -- ideal gas equation of state // areaNorm -- face area vector // positive -- flag to determine whether to add or subtract dissipation + // inp -- input variables // dot product of velocities with unit area vector const auto specRad = std::max(fabs(left.Velocity().DotProd(areaNorm)) + @@ -223,26 +115,39 @@ squareMatrix RusanovFluxJacobian(const primVars &left, const primVars &right, right.SoS(eos)); // form dissipation matrix based on spectral radius - squareMatrix dissipation(5); - dissipation.Identity(); - dissipation *= specRad; + fluxJacobian dissipation(inp.NumFlowEquations(), inp.NumTurbEquations()); + dissipation.flowJacobian_.Identity(); + dissipation.flowJacobian_ *= specRad; + // begin jacobian calculation - const auto fluxJac = positive ? - InvFluxJacobian(left, eos, areaNorm) : - InvFluxJacobian(right, eos, areaNorm); + positive ? this->InvFluxJacobian(left, eos, areaNorm, inp) : + this->InvFluxJacobian(right, eos, areaNorm, inp); + + + // compute turbulent jacobian if necessary + if (inp.IsTurbulent()) { + dissipation.turbJacobian_.Identity(); + + const auto turbSpecRad = std::max(fabs(left.Velocity().DotProd(areaNorm)), + fabs(right.Velocity().DotProd(areaNorm))); - return positive ? 0.5 * (fluxJac + dissipation) : - 0.5 * (fluxJac - dissipation); + dissipation.turbJacobian_ *= turbSpecRad; + } + + positive ? (*this) += dissipation : (*this) -= dissipation; + (*this) *= 0.5; } // function to calculate inviscid flux jacobian -squareMatrix InvFluxJacobian(const primVars &state, - const idealGas &eqnState, - const vector3d &areaNorm) { +void fluxJacobian::InvFluxJacobian(const primVars &state, + const idealGas &eqnState, + const vector3d &areaNorm, + const input &inp) { // state -- primative variables from left side // eqnState -- ideal gas equation of state // areaNorm -- face area vector + // inp -- input variables const auto velNorm = state.Velocity().DotProd(areaNorm); const auto gammaMinusOne = eqnState.Gamma() - 1.0; @@ -251,45 +156,57 @@ squareMatrix InvFluxJacobian(const primVars &state, const auto a3 = eqnState.Gamma() - 2.0; // begin jacobian calculation - squareMatrix A(5); + flowJacobian_ = squareMatrix(inp.NumFlowEquations()); + turbJacobian_ = squareMatrix(inp.NumTurbEquations()); // calculate flux derivatives wrt left state // column zero - A(0, 0) = 0.0; - A(1, 0) = phi * areaNorm.X() - state.U() * velNorm; - A(2, 0) = phi * areaNorm.Y() - state.V() * velNorm; - A(3, 0) = phi * areaNorm.Z() - state.W() * velNorm; - A(4, 0) = velNorm * (phi - a1); + flowJacobian_(0, 0) = 0.0; + flowJacobian_(1, 0) = phi * areaNorm.X() - state.U() * velNorm; + flowJacobian_(2, 0) = phi * areaNorm.Y() - state.V() * velNorm; + flowJacobian_(3, 0) = phi * areaNorm.Z() - state.W() * velNorm; + flowJacobian_(4, 0) = velNorm * (phi - a1); // column one - A(0, 1) = areaNorm.X(); - A(1, 1) = velNorm - a3 * areaNorm.X() * state.U(); - A(2, 1) = state.V() * areaNorm.X() - gammaMinusOne * state.U() * areaNorm.Y(); - A(3, 1) = state.W() * areaNorm.X() - gammaMinusOne * state.U() * areaNorm.Z(); - A(4, 1) = a1 * areaNorm.X() - gammaMinusOne * state.U() * velNorm; + flowJacobian_(0, 1) = areaNorm.X(); + flowJacobian_(1, 1) = velNorm - a3 * areaNorm.X() * state.U(); + flowJacobian_(2, 1) = state.V() * areaNorm.X() - + gammaMinusOne * state.U() * areaNorm.Y(); + flowJacobian_(3, 1) = state.W() * areaNorm.X() - + gammaMinusOne * state.U() * areaNorm.Z(); + flowJacobian_(4, 1) = a1 * areaNorm.X() - gammaMinusOne * state.U() * velNorm; // column two - A(0, 2) = areaNorm.Y(); - A(1, 2) = state.U() * areaNorm.Y() - gammaMinusOne * state.V() * areaNorm.X(); - A(2, 2) = velNorm - a3 * areaNorm.Y() * state.V(); - A(3, 2) = state.W() * areaNorm.Y() - gammaMinusOne * state.V() * areaNorm.Z(); - A(4, 2) = a1 * areaNorm.Y() - gammaMinusOne * state.V() * velNorm; + flowJacobian_(0, 2) = areaNorm.Y(); + flowJacobian_(1, 2) = state.U() * areaNorm.Y() - + gammaMinusOne * state.V() * areaNorm.X(); + flowJacobian_(2, 2) = velNorm - a3 * areaNorm.Y() * state.V(); + flowJacobian_(3, 2) = state.W() * areaNorm.Y() - + gammaMinusOne * state.V() * areaNorm.Z(); + flowJacobian_(4, 2) = a1 * areaNorm.Y() - gammaMinusOne * state.V() * velNorm; // column three - A(0, 3) = areaNorm.Z(); - A(1, 3) = state.U() * areaNorm.Z() - gammaMinusOne * state.W() * areaNorm.X(); - A(2, 3) = state.V() * areaNorm.Z() - gammaMinusOne * state.W() * areaNorm.Y(); - A(3, 3) = velNorm - a3 * areaNorm.Z() * state.W(); - A(4, 3) = a1 * areaNorm.Z() - gammaMinusOne * state.W() * velNorm; + flowJacobian_(0, 3) = areaNorm.Z(); + flowJacobian_(1, 3) = state.U() * areaNorm.Z() - + gammaMinusOne * state.W() * areaNorm.X(); + flowJacobian_(2, 3) = state.V() * areaNorm.Z() - + gammaMinusOne * state.W() * areaNorm.Y(); + flowJacobian_(3, 3) = velNorm - a3 * areaNorm.Z() * state.W(); + flowJacobian_(4, 3) = a1 * areaNorm.Z() - gammaMinusOne * state.W() * velNorm; // column four - A(0, 4) = 0.0; - A(1, 4) = gammaMinusOne * areaNorm.X(); - A(2, 4) = gammaMinusOne * areaNorm.Y(); - A(3, 4) = gammaMinusOne * areaNorm.Z(); - A(4, 4) = eqnState.Gamma() * velNorm; + flowJacobian_(0, 4) = 0.0; + flowJacobian_(1, 4) = gammaMinusOne * areaNorm.X(); + flowJacobian_(2, 4) = gammaMinusOne * areaNorm.Y(); + flowJacobian_(3, 4) = gammaMinusOne * areaNorm.Z(); + flowJacobian_(4, 4) = eqnState.Gamma() * velNorm; + + // turbulent jacobian here + if (inp.IsTurbulent()) { + // DO STUFF HERE + - return A; + } } /* Function to calculate approximate Roe flux jacobian. The Roe flux is @@ -306,29 +223,155 @@ jacobians. In the above equations the Roe matrix Aroe is held constant during differentiation. A represents the convective flux jacobian matrix. */ -squareMatrix ApproxRoeFluxJacobian(const primVars &left, const primVars &right, - const idealGas &eos, - const vector3d &areaNorm, - const bool &positive) { +void fluxJacobian::ApproxRoeFluxJacobian(const primVars &left, + const primVars &right, + const idealGas &eos, + const vector3d &areaNorm, + const bool &positive, + const input &inp) { // left -- primative variables from left side // right -- primative variables from right side // eos -- ideal gas equation of state // areaNorm -- face unit area vector // positive -- flag to determine whether to add or subtract dissipation + // inp -- input variables // compute Roe averaged state const auto roeAvg = RoeAveragedState(left, right, eos); // compute Roe matrix - const auto Aroe = InvFluxJacobian(roeAvg, eos, areaNorm); + fluxJacobian roeMatrix; + roeMatrix.InvFluxJacobian(roeAvg, eos, areaNorm, inp); // compute convective flux jacobian - const auto fluxJac = positive ? InvFluxJacobian(left, eos, areaNorm) : - InvFluxJacobian(right, eos, areaNorm); + positive ? this->InvFluxJacobian(left, eos, areaNorm, inp) : + this->InvFluxJacobian(right, eos, areaNorm, inp); + + positive ? (*this) += roeMatrix : (*this) -= roeMatrix; + (*this) *= 0.5; +} + +// change of variable matrix going frim primative to conservative variables +// from Dwight +void fluxJacobian::DelPrimativeDelConservative(const primVars &state, + const idealGas &eos, + const input &inp) { + // state -- primative variables + // eos -- equation of state + // inp -- input variables + + const auto gammaMinusOne = eos.Gamma() - 1.0; + const auto invRho = 1.0 / state.Rho(); + + flowJacobian_ = squareMatrix(inp.NumFlowEquations()); + turbJacobian_ = squareMatrix(inp.NumTurbEquations()); + + // assign first column + flowJacobian_(0, 0) = 1.0; + flowJacobian_(1, 0) = -invRho * state.U(); + flowJacobian_(2, 0) = -invRho * state.V(); + flowJacobian_(3, 0) = -invRho * state.W(); + flowJacobian_(4, 0) = 0.5 * gammaMinusOne * + state.Velocity().DotProd(state.Velocity()); + + // assign second column + flowJacobian_(1, 1) = invRho; + flowJacobian_(4, 1) = -gammaMinusOne * state.U(); + + // assign third column + flowJacobian_(2, 2) = invRho; + flowJacobian_(4, 2) = -gammaMinusOne * state.V(); + + // assign fourth column + flowJacobian_(3, 3) = invRho; + flowJacobian_(4, 3) = -gammaMinusOne * state.W(); + + // assign fifth column + flowJacobian_(4, 4) = gammaMinusOne; + + // turbulent jacobian here + if (inp.IsTurbulent()) { + // DO STUFF HERE + - return positive ? 0.5 * (fluxJac + Aroe) : 0.5 * (fluxJac - Aroe); + } } +// approximate thin shear layer jacobian following implementation in Dwight. +// does not use any gradients +void fluxJacobian::ApproxTSLJacobian(const primVars &state, const idealGas &eos, + const sutherland &suth, + const vector3d &area, + const double &dist, + const unique_ptr &turb, + const input &inp) { + // state -- primative variables + // eos -- equation of state + // suth -- sutherland's law for viscosity + // area -- face area unit vector + // dist -- distance from cell center to cell center + // turb -- turbulence model + // inp -- input variables + + flowJacobian_ = squareMatrix(inp.NumFlowEquations()); + turbJacobian_ = squareMatrix(inp.NumTurbEquations()); + + const auto mu = suth.EffectiveViscosity(state.Temperature(eos)); + const auto mut = turb->EddyViscNoLim(state) * suth.NondimScaling(); + const auto velNorm = state.Velocity().DotProd(area); + + // assign first column + flowJacobian_(4, 0) = -eos.Conductivity(mu + mut) * state.Temperature(eos) / + ((mu + mut) * state.Rho()); + + // assign second column + flowJacobian_(1, 1) = (1.0 / 3.0) * area.X() * area.X() + 1.0; + flowJacobian_(1, 2) = (1.0 / 3.0) * area.X() * area.Y(); + flowJacobian_(1, 3) = (1.0 / 3.0) * area.X() * area.Z(); + flowJacobian_(1, 4) = (1.0 / 3.0) * area.X() * velNorm + state.U(); + + // assign third column + flowJacobian_(2, 1) = (1.0 / 3.0) * area.Y() * area.X(); + flowJacobian_(2, 2) = (1.0 / 3.0) * area.Y() * area.Y() + 1.0; + flowJacobian_(2, 3) = (1.0 / 3.0) * area.Y() * area.Z(); + flowJacobian_(2, 4) = (1.0 / 3.0) * area.Y() * velNorm + state.V(); + + // assign fourth column + flowJacobian_(3, 1) = (1.0 / 3.0) * area.Z() * area.X(); + flowJacobian_(3, 2) = (1.0 / 3.0) * area.Z() * area.Y(); + flowJacobian_(3, 3) = (1.0 / 3.0) * area.Z() * area.Z() + 1.0; + flowJacobian_(3, 4) = (1.0 / 3.0) * area.Z() * velNorm + state.W(); + + // assign fifth column + flowJacobian_(4, 4) = eos.Conductivity(mu + mut) / ((mu + mut) * state.Rho()); + + flowJacobian_ *= (mu + mut) / dist; + + fluxJacobian prim2Cons; + prim2Cons.DelPrimativeDelConservative(state, eos, inp); + + flowJacobian_ = flowJacobian_.MatMult(prim2Cons.flowJacobian_); + turbJacobian_ = turbJacobian_.MatMult(prim2Cons.turbJacobian_); + + // calculate turbulent jacobian if necessary + if (inp.IsTurbulent()) { + // DO STUFF HERE + + + } +} + + +// non-member functions +// ---------------------------------------------------------------------------- +// operator overload for << - allows use of cout, cerr, etc. +ostream &operator<<(ostream &os, fluxJacobian &jacobian) { + os << jacobian.FlowJacobian() << endl; + os << jacobian.TurbulenceJacobian() << endl; + return os; +} + + genArray RusanovOffDiagonal(const primVars &state, const genArray &update, const unitVec3dMag &fAreaL, const unitVec3dMag &fAreaR, @@ -403,119 +446,23 @@ genArray RoeOffDiagonal(const primVars &left, const primVars &right, // add contribution for viscous terms uncoupledScalar specRad(0.0, 0.0); - // squareMatrix jac(5); if (inp.IsViscous()) { const auto offState = positive ? left : right; specRad.AddToFlowVariable( offState.ViscCellSpectralRadius(fAreaL, fAreaR, eos, suth, vol, turb)); - // jac = ApproxTSLJacobian(offState, eos, suth, areaNorm, 1.0, turb); - if (inp.IsTurbulent()) { specRad.AddToTurbVariable(turb->ViscSpecRad(offState, fAreaL, fAreaR, eos, suth, vol)); } } - // const auto invFluxJac = RusanovFluxJacobian(left, right, eos, areaNorm, - // positive); - - // don't need 0.5 factor on roe flux because RoeFlux function already does it return positive ? fAreaL.Mag() * ((newFlux - oldFlux).ConvertToGenArray()) + 0.5 * specRad.ArrayMult(update) : fAreaL.Mag() * ((newFlux - oldFlux).ConvertToGenArray()) - 0.5 * specRad.ArrayMult(update); - // fAreaL.Mag() * invFluxJac.ArrayMult(update) + 0.5 * specRadArr * update : - // fAreaL.Mag() * invFluxJac.ArrayMult(update) - 0.5 * specRadArr * update; } -// change of variable matrix going frim primative to conservative variables -// from Dwight -squareMatrix DelPrimativeDelConservative(const primVars &state, - const idealGas &eos) { - // state -- primative variables - // eos -- equation of state - - const auto gammaMinusOne = eos.Gamma() - 1.0; - const auto invRho = 1.0 / state.Rho(); - - squareMatrix dPdC(5); - - // assign first column - dPdC(0, 0) = 1.0; - dPdC(1, 0) = -invRho * state.U(); - dPdC(2, 0) = -invRho * state.V(); - dPdC(3, 0) = -invRho * state.W(); - dPdC(4, 0) = 0.5 * gammaMinusOne * state.Velocity().DotProd(state.Velocity()); - - // assign second column - dPdC(1, 1) = invRho; - dPdC(4, 1) = -gammaMinusOne * state.U(); - - // assign third column - dPdC(2, 2) = invRho; - dPdC(4, 2) = -gammaMinusOne * state.V(); - - // assign fourth column - dPdC(3, 3) = invRho; - dPdC(4, 3) = -gammaMinusOne * state.W(); - - // assign fifth column - dPdC(4, 4) = gammaMinusOne; - - return dPdC; -} - -// approximate thin shear layer jacobian following implementation in Dwight. -// does not use any gradients -squareMatrix ApproxTSLJacobian(const primVars &state, const idealGas &eos, - const sutherland &suth, - const vector3d &area, const double &dist, - const unique_ptr &turb) { - // state -- primative variables - // eos -- equation of state - // suth -- sutherland's law for viscosity - // area -- face area unit vector - // dist -- distance from cell center to cell center - // turb -- turbulence model - - squareMatrix jacobian(5); - - const auto mu = suth.EffectiveViscosity(state.Temperature(eos)); - const auto mut = turb->EddyViscNoLim(state) * suth.NondimScaling(); - const auto velNorm = state.Velocity().DotProd(area); - - // assign first column - jacobian(4, 0) = -eos.Conductivity(mu + mut) * state.Temperature(eos) / - ((mu + mut) * state.Rho()); - - // assign second column - jacobian(1, 1) = (1.0 / 3.0) * area.X() * area.X() + 1.0; - jacobian(1, 2) = (1.0 / 3.0) * area.X() * area.Y(); - jacobian(1, 3) = (1.0 / 3.0) * area.X() * area.Z(); - jacobian(1, 4) = (1.0 / 3.0) * area.X() * velNorm + state.U(); - - // assign third column - jacobian(2, 1) = (1.0 / 3.0) * area.Y() * area.X(); - jacobian(2, 2) = (1.0 / 3.0) * area.Y() * area.Y() + 1.0; - jacobian(2, 3) = (1.0 / 3.0) * area.Y() * area.Z(); - jacobian(2, 4) = (1.0 / 3.0) * area.Y() * velNorm + state.V(); - - // assign fourth column - jacobian(3, 1) = (1.0 / 3.0) * area.Z() * area.X(); - jacobian(3, 2) = (1.0 / 3.0) * area.Z() * area.Y(); - jacobian(3, 3) = (1.0 / 3.0) * area.Z() * area.Z() + 1.0; - jacobian(3, 4) = (1.0 / 3.0) * area.Z() * velNorm + state.W(); - - // assign fifth column - jacobian(4, 4) = eos.Conductivity(mu + mut) / ((mu + mut) * state.Rho()); - - jacobian *= (mu + mut) / dist; - - const auto delPrimDelCons = DelPrimativeDelConservative(state, eos); - - return jacobian.MatMult(delPrimDelCons); -} diff --git a/fluxJacobian.hpp b/fluxJacobian.hpp index 18bc689..c7887e3 100644 --- a/fluxJacobian.hpp +++ b/fluxJacobian.hpp @@ -56,10 +56,6 @@ class fluxJacobian { fluxJacobian() : fluxJacobian(0.0, 0.0) {} explicit fluxJacobian(const uncoupledScalar &specRad) : fluxJacobian(specRad.FlowVariable(), specRad.TurbVariable()) {} - fluxJacobian(const primVars &, const unitVec3dMag &, - const unitVec3dMag &, const idealGas &, - const sutherland &, const double &, - const unique_ptr &, const input &, const bool &); // move constructor and assignment operator fluxJacobian(fluxJacobian&&) noexcept = default; @@ -73,15 +69,21 @@ class fluxJacobian { squareMatrix FlowJacobian() const { return flowJacobian_; } squareMatrix TurbulenceJacobian() const { return turbJacobian_; } - void AddInviscidJacobian(const primVars &, const unitVec3dMag &, - const unitVec3dMag &, const idealGas &, - const unique_ptr &, const bool &); - void AddViscousJacobian(const primVars &, const unitVec3dMag &, - const unitVec3dMag &, const idealGas &, - const sutherland &, const double &, - const unique_ptr &, const bool &); - void AddTurbSourceJacobian(const primVars &, const sutherland &, - const double &, const unique_ptr &); + void RusanovFluxJacobian(const primVars &, const primVars &, + const idealGas &, const vector3d &, + const bool &, const input &); + void InvFluxJacobian(const primVars &, const idealGas &, + const vector3d &, const input &); + void ApproxRoeFluxJacobian(const primVars &, const primVars &, + const idealGas &, const vector3d &, + const bool &, const input &); + void DelPrimativeDelConservative(const primVars &, const idealGas &, + const input &); + + void ApproxTSLJacobian(const primVars &, const idealGas &, + const sutherland &, + const vector3d &, const double &, + const unique_ptr &, const input &); void Zero() { flowJacobian_.Zero(); @@ -225,14 +227,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 &, - const bool &); -squareMatrix InvFluxJacobian(const primVars &, const idealGas &, - const vector3d &); -squareMatrix ApproxRoeFluxJacobian(const primVars &, const primVars &, - const idealGas &, const vector3d &, - const bool &); genArray RusanovOffDiagonal(const primVars &, const genArray &, const unitVec3dMag &, const unitVec3dMag &, @@ -247,11 +241,5 @@ genArray RoeOffDiagonal(const primVars &, const primVars &, const genArray &, const unique_ptr &, const input &, const bool &); -squareMatrix DelPrimativeDelConservative(const primVars &, const idealGas &); - -squareMatrix ApproxTSLJacobian(const primVars &, const idealGas &, - const sutherland &, - const vector3d &, const double &, - const unique_ptr &); #endif diff --git a/input.cpp b/input.cpp index a866a8e..09aa41b 100644 --- a/input.cpp +++ b/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; diff --git a/input.hpp b/input.hpp index c9c3dcc..e568f83 100644 --- a/input.hpp +++ b/input.hpp @@ -23,6 +23,7 @@ #include // unique_ptr #include "vector3d.hpp" #include "boundaryConditions.hpp" +#include "macros.hpp" using std::vector; using std::string; @@ -161,6 +162,8 @@ class input { int NumVars() const {return vars_.size();} int NumEquations() const; + int NumFlowEquations() const {return NUMFLOWVARS;} + int NumTurbEquations() const; void ReadInput(const int &); @@ -168,7 +171,7 @@ class input { bool IsViscous() const; bool IsTurbulent() const; bool IsBlockMatrix() const; - + string OrderOfAccuracy() const; double FarfieldTurbIntensity() const {return farfieldTurbInten_;} diff --git a/macros.hpp b/macros.hpp index 96be0c9..4fbb0be 100644 --- a/macros.hpp +++ b/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 diff --git a/procBlock.cpp b/procBlock.cpp index e7b1738..f560e02 100644 --- a/procBlock.cpp +++ b/procBlock.cpp @@ -332,9 +332,10 @@ void procBlock::CalcInvFluxI(const idealGas &eqnState, const input &inp, // if using a block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { - mainDiagonal(ip, jp, kp).AddInviscidJacobian( - state_(ig, jg, kg), fAreaI_(ig, jg, kg), - fAreaI_(ig + 1, jg, kg), eqnState, turb, inp.IsTurbulent()); + // fluxJacobian fluxJac; + // fluxJac.RusanovFluxJacobian(state_(ig - 1, jg, kg), state_(ig, jg, kg), + // eqnState, fAreaI_(ig, jg, kg), true, inp); + // mainDiagonal(ip, jp, kp) += fluxJac; } } } @@ -450,10 +451,10 @@ void procBlock::CalcInvFluxJ(const idealGas &eqnState, const input &inp, // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { - mainDiagonal(ip, jp, kp).AddInviscidJacobian( - state_(ig, jg, kg), fAreaJ_(ig, jg, kg), - fAreaJ_(ig, jg + 1, kg), eqnState, turb, - inp.IsTurbulent()); + // mainDiagonal(ip, jp, kp).AddInviscidJacobian( + // state_(ig, jg, kg), fAreaJ_(ig, jg, kg), + // fAreaJ_(ig, jg + 1, kg), eqnState, turb, + // inp.IsTurbulent()); } } } @@ -570,10 +571,10 @@ void procBlock::CalcInvFluxK(const idealGas &eqnState, const input &inp, // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { - mainDiagonal(ip, jp, kp).AddInviscidJacobian( - state_(ig, jg, kg), fAreaK_(ig, jg, kg), - fAreaK_(ig, jg, kg + 1), eqnState, turb, - inp.IsTurbulent()); + // mainDiagonal(ip, jp, kp).AddInviscidJacobian( + // state_(ig, jg, kg), fAreaK_(ig, jg, kg), + // fAreaK_(ig, jg, kg + 1), eqnState, turb, + // inp.IsTurbulent()); } } } @@ -1944,10 +1945,10 @@ void procBlock::CalcViscFluxI(const sutherland &suth, const idealGas &eqnState, // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { - mainDiagonal(ip, jp, kp).AddViscousJacobian( - state_(ig, jg, kg), fAreaI_(ig, jg, kg), - fAreaI_(ig + 1, jg, kg), eqnState, suth, - vol_(ig, jg, kg), turb, inp.IsTurbulent()); + // mainDiagonal(ip, jp, kp).AddViscousJacobian( + // state_(ig, jg, kg), fAreaI_(ig, jg, kg), + // fAreaI_(ig + 1, jg, kg), eqnState, suth, + // vol_(ig, jg, kg), turb, inp.IsTurbulent()); } } } @@ -2108,10 +2109,10 @@ void procBlock::CalcViscFluxJ(const sutherland &suth, const idealGas &eqnState, // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { - mainDiagonal(ip, jp, kp).AddViscousJacobian( - state_(ig, jg, kg), fAreaJ_(ig, jg, kg), - fAreaJ_(ig, jg + 1, kg), eqnState, suth, - vol_(ig, jg, kg), turb, inp.IsTurbulent()); + // mainDiagonal(ip, jp, kp).AddViscousJacobian( + // state_(ig, jg, kg), fAreaJ_(ig, jg, kg), + // fAreaJ_(ig, jg + 1, kg), eqnState, suth, + // vol_(ig, jg, kg), turb, inp.IsTurbulent()); } } } @@ -2271,10 +2272,10 @@ void procBlock::CalcViscFluxK(const sutherland &suth, const idealGas &eqnState, // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { - mainDiagonal(ip, jp, kp).AddViscousJacobian( - state_(ig, jg, kg), fAreaK_(ig, jg, kg), - fAreaK_(ig, jg, kg + 1), eqnState, suth, - vol_(ig, jg, kg), turb, inp.IsTurbulent()); + // mainDiagonal(ip, jp, kp).AddViscousJacobian( + // state_(ig, jg, kg), fAreaK_(ig, jg, kg), + // fAreaK_(ig, jg, kg + 1), eqnState, suth, + // vol_(ig, jg, kg), turb, inp.IsTurbulent()); } } } @@ -6634,8 +6635,8 @@ void procBlock::CalcSrcTerms(const gradients &grads, const sutherland &suth, // add contribution of source spectral radius to flux jacobian if (inp.IsBlockMatrix()) { - mainDiagonal(ip, jp, kp).AddTurbSourceJacobian( - state_(ig, jg, kg), suth, vol_(ig, jg, kg), turb); + // mainDiagonal(ip, jp, kp).AddTurbSourceJacobian( + // state_(ig, jg, kg), suth, vol_(ig, jg, kg), turb); } } } From 8a8c8ba90063e0ae42d47fda3b433f4e85b0ee41 Mon Sep 17 00:00:00 2001 From: Michael Date: Tue, 26 Apr 2016 22:45:56 -0700 Subject: [PATCH 3/5] Changed main diagonal to accumulate differently from spectral radius. Still need to add in turbulence jacobians for full matrix functions. --- fluxJacobian.hpp | 5 ++ procBlock.cpp | 157 ++++++++++++++++++++++++++--------------------- 2 files changed, 93 insertions(+), 69 deletions(-) diff --git a/fluxJacobian.hpp b/fluxJacobian.hpp index c7887e3..97e9ad6 100644 --- a/fluxJacobian.hpp +++ b/fluxJacobian.hpp @@ -69,6 +69,11 @@ class fluxJacobian { 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 &, const bool &, const input &); diff --git a/procBlock.cpp b/procBlock.cpp index f560e02..0adab53 100644 --- a/procBlock.cpp +++ b/procBlock.cpp @@ -320,15 +320,15 @@ void procBlock::CalcInvFluxI(const idealGas &eqnState, const input &inp, ip, jp, kp); // calculate component of wave speed. This is done on a cell by cell // basis, so only at the upper faces - specRadius_(ip, jp, kp).AddToFlowVariable( - state_(ig, jg, kg).InvCellSpectralRadius(fAreaI_(ig, jg, kg), - fAreaI_(ig + 1, jg, kg), - eqnState)); - if (inp.IsTurbulent()) { - specRadius_(ip, jp, kp).AddToTurbVariable( - turb->InviscidSpecRad(state_(ig, jg, kg), fAreaI_(ig, jg, kg), - fAreaI_(ig + 1, jg, kg))); - } + const auto invSpecRad = state_(ig, jg, kg).InvCellSpectralRadius( + fAreaI_(ig, jg, kg), fAreaI_(ig + 1, jg, kg), eqnState); + + const auto turbInvSpecRad = inp.IsTurbulent() ? + turb->InviscidSpecRad(state_(ig, jg, kg), fAreaI_(ig, jg, kg), + fAreaI_(ig + 1, jg, kg)): 0.0; + + const uncoupledScalar specRad(invSpecRad, turbInvSpecRad); + specRadius_(ip, jp, kp) += specRad; // if using a block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { @@ -336,6 +336,8 @@ void procBlock::CalcInvFluxI(const idealGas &eqnState, const input &inp, // fluxJac.RusanovFluxJacobian(state_(ig - 1, jg, kg), state_(ig, jg, kg), // eqnState, fAreaI_(ig, jg, kg), true, inp); // mainDiagonal(ip, jp, kp) += fluxJac; + } else { + mainDiagonal(ip, jp, kp) += fluxJacobian(specRad); } } } @@ -439,15 +441,15 @@ void procBlock::CalcInvFluxJ(const idealGas &eqnState, const input &inp, // calculate component of wave speed. This is done on a cell by cell // basis, so only at the upper faces - specRadius_(ip, jp, kp).AddToFlowVariable( - state_(ig, jg, kg).InvCellSpectralRadius(fAreaJ_(ig, jg, kg), - fAreaJ_(ig, jg + 1, kg), - eqnState)); - if (inp.IsTurbulent()) { - specRadius_(ip, jp, kp).AddToTurbVariable( - turb->InviscidSpecRad(state_(ig, jg, kg), fAreaJ_(ig, jg, kg), - fAreaJ_(ig, jg + 1, kg))); - } + const auto invSpecRad = state_(ig, jg, kg).InvCellSpectralRadius( + fAreaJ_(ig, jg, kg), fAreaJ_(ig, jg + 1, kg), eqnState); + + const auto turbInvSpecRad = inp.IsTurbulent() ? + turb->InviscidSpecRad(state_(ig, jg, kg), fAreaJ_(ig, jg, kg), + fAreaJ_(ig, jg + 1, kg)): 0.0; + + const uncoupledScalar specRad(invSpecRad, turbInvSpecRad); + specRadius_ += specRad; // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { @@ -455,6 +457,8 @@ void procBlock::CalcInvFluxJ(const idealGas &eqnState, const input &inp, // state_(ig, jg, kg), fAreaJ_(ig, jg, kg), // fAreaJ_(ig, jg + 1, kg), eqnState, turb, // inp.IsTurbulent()); + } else { + mainDiagonal(ip, jp, kp) += fluxJacobian(specRad); } } } @@ -559,15 +563,15 @@ void procBlock::CalcInvFluxK(const idealGas &eqnState, const input &inp, // calculate component of wave speed. This is done on a cell by cell // basis, so only at the upper faces - specRadius_(ip, jp, kp).AddToFlowVariable( - state_(ig, jg, kg).InvCellSpectralRadius(fAreaK_(ig, jg, kg), - fAreaK_(ig, jg, kg + 1), - eqnState)); - if (inp.IsTurbulent()) { - specRadius_(ip, jp, kp).AddToTurbVariable( - turb->InviscidSpecRad(state_(ig, jg, kg), fAreaK_(ig, jg, kg), - fAreaK_(ig + 1, jg, kg + 1))); - } + const auto invSpecRad = state_(ig, jg, kg).InvCellSpectralRadius( + fAreaK_(ig, jg, kg), fAreaK_(ig, jg, kg + 1), eqnState); + + const auto turbInvSpecRad = inp.IsTurbulent() ? + turb->InviscidSpecRad(state_(ig, jg, kg), fAreaK_(ig, jg, kg), + fAreaK_(ig + 1, jg, kg + 1)) : 0.0; + + const uncoupledScalar specRad(invSpecRad, turbInvSpecRad); + specRadius_(ip, jp, kp) += specRad; // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { @@ -575,6 +579,8 @@ void procBlock::CalcInvFluxK(const idealGas &eqnState, const input &inp, // state_(ig, jg, kg), fAreaK_(ig, jg, kg), // fAreaK_(ig, jg, kg + 1), eqnState, turb, // inp.IsTurbulent()); + } else { + mainDiagonal(ip, jp, kp) += fluxJacobian(specRad); } } } @@ -922,12 +928,6 @@ void procBlock::InvertDiagonal(multiArray3d &mainDiagonal, inp.DualTimeCFL(); } - // if not using a block matrix on diagonal, construct from - // spectral radii - if (!inp.IsBlockMatrix()) { - mainDiagonal(ip, jp, kp) = fluxJacobian(specRadius_(ip, jp, kp)); - } - // add volume and time term mainDiagonal(ip, jp, kp) *= inp.MatrixRelaxation(); mainDiagonal(ip, jp, kp) += diagVolTime; @@ -1877,6 +1877,8 @@ void procBlock::CalcViscFluxI(const sutherland &suth, const idealGas &eqnState, // mainDiagonal -- main diagonal of LHS used to store flux jacobians for // implicit solver + const auto viscCoeff = inp.ViscousCFLCoefficient(); + // loop over all physical i-faces for (auto kp = 0, kg = numGhosts_; kg < fAreaI_.NumK() - numGhosts_; kg++, kp++) { @@ -1929,19 +1931,19 @@ void procBlock::CalcViscFluxI(const sutherland &suth, const idealGas &eqnState, // calculate component of wave speed. This is done on a cell by cell // basis, so only at the upper faces - // factor 2 because visc spectral radius is not halved (Blazek 6.53) - specRadius_(ip, jp, kp).AddToFlowVariable( - 2.0 * state_(ig, jg, kg).ViscCellSpectralRadius( + const auto viscSpecRad = + state_(ig, jg, kg).ViscCellSpectralRadius( fAreaI_(ig, jg, kg), fAreaI_(ig + 1, jg, kg), eqnState, suth, - vol_(ig, jg, kg), turb)); + vol_(ig, jg, kg), turb); + + const auto turbViscSpecRad = inp.IsTurbulent() ? + turb->ViscSpecRad(state_(ig, jg, kg), fAreaI_(ig, jg, kg), + fAreaI_(ig + 1, jg, kg), eqnState, suth, + vol_(ig, jg, kg)) : 0.0; + + const uncoupledScalar specRad(viscSpecRad, turbViscSpecRad); + specRadius_(ip, jp, kp) += specRad * viscCoeff; - if (inp.IsTurbulent()) { - // factor 2 because visc spectral radius is not halved (Blazek 6.53) - specRadius_(ip, jp, kp).AddToTurbVariable( - 2.0 * turb->ViscSpecRad(state_(ig, jg, kg), fAreaI_(ig, jg, kg), - fAreaI_(ig + 1, jg, kg), eqnState, suth, - vol_(ig, jg, kg))); - } // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { @@ -1949,6 +1951,9 @@ void procBlock::CalcViscFluxI(const sutherland &suth, const idealGas &eqnState, // state_(ig, jg, kg), fAreaI_(ig, jg, kg), // fAreaI_(ig + 1, jg, kg), eqnState, suth, // vol_(ig, jg, kg), turb, inp.IsTurbulent()); + } else { + // factor 2 because visc spectral radius is not halved (Blazek 6.53) + mainDiagonal(ip, jp, kp) += fluxJacobian(2.0 * specRad); } } } @@ -2041,6 +2046,8 @@ void procBlock::CalcViscFluxJ(const sutherland &suth, const idealGas &eqnState, // mainDiagonal -- main diagonal of LHS used to store flux jacobians for // implicit solver + const auto viscCoeff = inp.ViscousCFLCoefficient(); + // loop over all physical j-faces for (auto kp = 0, kg = numGhosts_; kg < fAreaJ_.NumK() - numGhosts_; kg++, kp++) { @@ -2093,19 +2100,19 @@ void procBlock::CalcViscFluxJ(const sutherland &suth, const idealGas &eqnState, // calculate component of wave speed. This is done on a cell by cell // basis, so only at the upper faces - // factor 2 because visc spectral radius is not halved (Blazek 6.53) - specRadius_(ip, jp, kp).AddToFlowVariable( - 2.0 * state_(ig, jg, kg).ViscCellSpectralRadius( + const auto viscSpecRad = + state_(ig, jg, kg).ViscCellSpectralRadius( fAreaJ_(ig, jg, kg), fAreaJ_(ig, jg + 1, kg), eqnState, suth, - vol_(ig, jg, kg), turb)); + vol_(ig, jg, kg), turb); + + const auto turbViscSpecRad = inp.IsTurbulent() ? + turb->ViscSpecRad(state_(ig, jg, kg), fAreaJ_(ig, jg, kg), + fAreaJ_(ig, jg + 1, kg), eqnState, suth, + vol_(ig, jg, kg)) : 0.0; + + const uncoupledScalar specRad(viscSpecRad, turbViscSpecRad); + specRadius_(ip, jp, kp) += specRad * viscCoeff; - if (inp.IsTurbulent()) { - // factor 2 because visc spectral radius is not halved (Blazek 6.53) - specRadius_(ip, jp, kp).AddToTurbVariable( - 2.0 * turb->ViscSpecRad(state_(ig, jg, kg), fAreaJ_(ig, jg, kg), - fAreaJ_(ig, jg + 1, kg), eqnState, suth, - vol_(ig, jg, kg))); - } // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { @@ -2113,6 +2120,9 @@ void procBlock::CalcViscFluxJ(const sutherland &suth, const idealGas &eqnState, // state_(ig, jg, kg), fAreaJ_(ig, jg, kg), // fAreaJ_(ig, jg + 1, kg), eqnState, suth, // vol_(ig, jg, kg), turb, inp.IsTurbulent()); + } else { + // factor 2 because visc spectral radius is not halved (Blazek 6.53) + mainDiagonal(ip, jp, kp) += fluxJacobian(2.0 * specRad); } } } @@ -2204,6 +2214,8 @@ void procBlock::CalcViscFluxK(const sutherland &suth, const idealGas &eqnState, // mainDiagonal -- main diagonal of LHS used to store flux jacobians for // implicit solver + const auto viscCoeff = inp.ViscousCFLCoefficient(); + // loop over all physical k-faces for (auto kp = 0, kg = numGhosts_; kg < fAreaK_.NumK() - numGhosts_; kg++, kp++) { @@ -2256,19 +2268,19 @@ void procBlock::CalcViscFluxK(const sutherland &suth, const idealGas &eqnState, // calculate component of wave speed. This is done on a cell by cell // basis, so only at the upper faces - // factor 2 because visc spectral radius is not halved (Blazek 6.53) - specRadius_(ip, jp, kp).AddToFlowVariable( - 2.0 * state_(ig, jg, kg).ViscCellSpectralRadius( + const auto viscSpecRad = + state_(ig, jg, kg).ViscCellSpectralRadius( fAreaK_(ig, jg, kg), fAreaK_(ig, jg, kg + 1), eqnState, suth, - vol_(ig, jg, kg), turb)); - - if (inp.IsTurbulent()) { - // factor 2 because visc spectral radius is not halved (Blazek 6.53) - specRadius_(ip, jp, kp).AddToTurbVariable( - 2.0 * turb->ViscSpecRad(state_(ig, jg, kg), fAreaK_(ig, jg, kg), - fAreaK_(ig, jg, kg + 1), eqnState, suth, - vol_(ig, jg, kg))); - } + vol_(ig, jg, kg), turb); + + const auto turbViscSpecRad = inp.IsTurbulent() ? + turb->ViscSpecRad(state_(ig, jg, kg), fAreaK_(ig, jg, kg), + fAreaK_(ig, jg, kg + 1), eqnState, suth, + vol_(ig, jg, kg)) : 0.0; + + const uncoupledScalar specRad(viscSpecRad, turbViscSpecRad); + specRadius_(ip, jp, kp) += specRad * viscCoeff; + // if using block matrix on main diagonal, calculate flux jacobian if (inp.IsBlockMatrix()) { @@ -2276,6 +2288,9 @@ void procBlock::CalcViscFluxK(const sutherland &suth, const idealGas &eqnState, // state_(ig, jg, kg), fAreaK_(ig, jg, kg), // fAreaK_(ig, jg, kg + 1), eqnState, suth, // vol_(ig, jg, kg), turb, inp.IsTurbulent()); + } else { + // factor 2 because visc spectral radius is not halved (Blazek 6.53) + mainDiagonal(ip, jp, kp) += fluxJacobian(2.0 * specRad); } } } @@ -6630,13 +6645,17 @@ void procBlock::CalcSrcTerms(const gradients &grads, const sutherland &suth, ip, jp, kp); // add source spectral radius for turbulence equations - specRadius_(ip, jp, kp).SubtractFromTurbVariable( - turb->SrcSpecRad(state_(ig, jg, kg), suth, vol_(ig, jg, kg))); + const auto srcSpecRad = turb->SrcSpecRad(state_(ig, jg, kg), suth, + vol_(ig, jg, kg)); + specRadius_(ip, jp, kp).SubtractFromTurbVariable(srcSpecRad); // add contribution of source spectral radius to flux jacobian if (inp.IsBlockMatrix()) { // mainDiagonal(ip, jp, kp).AddTurbSourceJacobian( // state_(ig, jg, kg), suth, vol_(ig, jg, kg), turb); + } else { + const uncoupledScalar srcJac(0.0, srcSpecRad); + mainDiagonal(ip, jp, kp) -= fluxJacobian(srcJac); } } } From 4dfaa224ad51e22c26a3dbe9dcccee435c953232 Mon Sep 17 00:00:00 2001 From: Michael Date: Wed, 27 Apr 2016 23:50:38 -0700 Subject: [PATCH 4/5] Fixed full matrix fluxJacobian calculation for turbulent flow. Changed K-Omega Wilcox model to use limited eddy viscosity in viscous flux calculation as it should. --- fluxJacobian.cpp | 17 ++++++----- turbulence.cpp | 51 ++++++++++++-------------------- turbulence.hpp | 75 ++++++++++++++++++++++-------------------------- 3 files changed, 61 insertions(+), 82 deletions(-) diff --git a/fluxJacobian.cpp b/fluxJacobian.cpp index 9fdb19e..2b5ba57 100644 --- a/fluxJacobian.cpp +++ b/fluxJacobian.cpp @@ -203,9 +203,8 @@ void fluxJacobian::InvFluxJacobian(const primVars &state, // turbulent jacobian here if (inp.IsTurbulent()) { - // DO STUFF HERE - - + turbJacobian_(0, 0) = velNorm; + turbJacobian_(1, 1) = velNorm; } } @@ -291,9 +290,8 @@ void fluxJacobian::DelPrimativeDelConservative(const primVars &state, // turbulent jacobian here if (inp.IsTurbulent()) { - // DO STUFF HERE - - + turbJacobian_(0, 0) = invRho; + turbJacobian_(1, 1) = invRho; } } @@ -351,13 +349,14 @@ void fluxJacobian::ApproxTSLJacobian(const primVars &state, const idealGas &eos, prim2Cons.DelPrimativeDelConservative(state, eos, inp); flowJacobian_ = flowJacobian_.MatMult(prim2Cons.flowJacobian_); - turbJacobian_ = turbJacobian_.MatMult(prim2Cons.turbJacobian_); // calculate turbulent jacobian if necessary if (inp.IsTurbulent()) { - // DO STUFF HERE - + turbJacobian_(0, 0) = mu + turb->SigmaK() * mut; + turbJacobian_(1, 1) = mu + turb->SigmaW() * mut; + turbJacobian_ *= 1.0 / dist; + turbJacobian_ = turbJacobian_.MatMult(prim2Cons.turbJacobian_); } } diff --git a/turbulence.cpp b/turbulence.cpp index 196e887..fe705c1 100644 --- a/turbulence.cpp +++ b/turbulence.cpp @@ -131,6 +131,23 @@ double turbModel::SpectralRadius(const primVars &state, return specRad; } +// member function to calculate inviscid spectral radius +// df_dq = [vel (dot) area 0 +// 0 vel (dot) area] +double turbModel::InviscidSpecRad(const primVars &state, + const unitVec3dMag &fAreaL, + const unitVec3dMag &fAreaR) const { + // state -- primative variables + // fAreaL -- face area for left face + // fAreaR -- face area for right face + + auto normAvg = (0.5 * (fAreaL.UnitVector() + + fAreaR.UnitVector())).Normalize(); + auto fMag = 0.5 * (fAreaL.Mag() + fAreaR.Mag()); + return state.Velocity().DotProd(normAvg) * fMag; +} + + // ------------------------------------------------------------------------- // member functions for the turbNone class @@ -325,8 +342,8 @@ double turbKWWilcox::EddyViscAndMolecDiffCoeff(const primVars &state, sigmaK = sigmaStar_; sigmaW = sigma_; - // return eddy viscosity without limiter, scaled for nondimensional equations - return this->EddyViscNoLim(state) * suth.NondimScaling(); + // return eddy viscosity, scaled for nondimensional equations + return this->EddyVisc(state, velGrad, suth, 0.0) * suth.NondimScaling(); } /* member function to calculate the spectral radius of the source jacobian @@ -349,21 +366,6 @@ double turbKWWilcox::SrcSpecRad(const primVars &state, return -2.0 * betaStar_ * state.Omega() * vol * suth.InvNondimScaling(); } -// member function to calculate inviscid spectral radius -// df_dq = [vel (dot) area 0 -// 0 vel (dot) area] -double turbKWWilcox::InviscidSpecRad(const primVars &state, - const unitVec3dMag &fAreaL, - const unitVec3dMag &fAreaR) const { - // state -- primative variables - // fAreaL -- face area for left face - // fAreaR -- face area for right face - auto normAvg = (0.5 * (fAreaL.UnitVector() + - fAreaR.UnitVector())).Normalize(); - auto fMag = 0.5 * (fAreaL.Mag() + fAreaR.Mag()); - return state.Velocity().DotProd(normAvg) * fMag; -} - // member function to calculate viscous spectral radius // dfv_dq = [ (area / vol) * (nu + sigmaStar * nut) 0 // 0 (area / vol) * (nu + sigma * nut)] @@ -592,21 +594,6 @@ double turbKWSst::SrcSpecRad(const primVars &state, return -2.0 * betaStar_ * state.Omega() * vol * suth.InvNondimScaling(); } -// member function to calculate inviscid spectral radius -// df_dq = [vel (dot) area 0 -// 0 vel (dot) area] -double turbKWSst::InviscidSpecRad(const primVars &state, - const unitVec3dMag &fAreaL, - const unitVec3dMag &fAreaR) const { - // state -- primative variables - // fAreaL -- face area for left face - // fAreaR -- face area for right face - - auto normAvg = (0.5 * (fAreaL.UnitVector() + - fAreaR.UnitVector())).Normalize(); - auto fMag = 0.5 * (fAreaL.Mag() + fAreaR.Mag()); - return state.Velocity().DotProd(normAvg) * fMag; -} // member function to calculate viscous spectral radius // dfv_dq = [ (area / vol) * (nu + sigmaStar * nut) 0 diff --git a/turbulence.hpp b/turbulence.hpp index 7b18d95..8e8c15b 100644 --- a/turbulence.hpp +++ b/turbulence.hpp @@ -57,6 +57,24 @@ class turbModel { virtual double TurbPrandtlNumber() const {return 0.9;} virtual double TkeMin() const {return 1.0e-20;} virtual double OmegaMin() const {return 1.0e-20;} + virtual double SigmaK() const {return 0.0;} + virtual double SigmaW() const {return 0.0;} + virtual double EddyVisc(const primVars &state, + const tensor &vGrad, + const sutherland &suth, + const double &f2) const {return 0.0;} + virtual double WallBeta() const {return 1.0;} + virtual double SrcSpecRad(const primVars &state, + const sutherland &suth, + const double &vol) const {return 0.0;} + virtual double InviscidSpecRad(const primVars &state, + const unitVec3dMag &fAreaL, + const unitVec3dMag &fAreaR) const; + virtual double ViscSpecRad(const primVars &state, + const unitVec3dMag &fAreaL, + const unitVec3dMag &fAreaR, + const idealGas &eos, const sutherland &suth, + const double &vol) const {return 0.0;} tensor BoussinesqReynoldsStress(const primVars &state, const tensor &velGrad, const sutherland &suth, @@ -85,9 +103,6 @@ class turbModel { const sutherland &suth, const idealGas &eos, const double &wallDist, const double &vol, double &ksrc, double &wsrc) const = 0; - virtual double EddyVisc(const primVars &state, - const tensor &vGrad, - const sutherland &suth, const double &f2) const = 0; virtual double EddyViscAndMolecDiffCoeff(const primVars &state, const tensor &vGrad, const vector3d &kGrad, @@ -97,20 +112,6 @@ class turbModel { const double &wallDist, double &sigmaK, double &sigmaW) const = 0; - virtual double WallBeta() const = 0; - virtual double SrcSpecRad(const primVars &state, - const sutherland &suth, - const double &vol) const = 0; - virtual double InviscidSpecRad(const primVars &state, - const unitVec3dMag &fAreaL, - const unitVec3dMag &fAreaR) - const = 0; - virtual double ViscSpecRad(const primVars &state, - const unitVec3dMag &fAreaL, - const unitVec3dMag &fAreaR, - const idealGas &eos, const sutherland &suth, - const double &vol) const = 0; - virtual void Print() const = 0; // destructor @@ -138,23 +139,6 @@ class turbNone : public turbModel { const sutherland &suth, const idealGas &eos, const double &wallDist, const double &vol, double &ksrc, double &wsrc) const override; - double SrcSpecRad(const primVars &state, - const sutherland &suth, - const double &vol) const override {return 0.0;} - double InviscidSpecRad(const primVars &state, - const unitVec3dMag &fAreaL, - const unitVec3dMag &fAreaR) const override - {return 0.0;} - double ViscSpecRad(const primVars &state, - const unitVec3dMag &fAreaL, - const unitVec3dMag &fAreaR, - const idealGas &eos, const sutherland &suth, - const double &vol) const override {return 0.0;} - - double EddyVisc(const primVars &state, const tensor &vGrad, - const sutherland &suth, - const double &f2) const override {return 0.0;} - double EddyViscNoLim(const primVars &state) const override {return 0.0;} double EddyViscAndMolecDiffCoeff(const primVars &state, const tensor &velGrad, const vector3d &kGrad, @@ -163,7 +147,13 @@ class turbNone : public turbModel { const idealGas &eos, const double &wallDist, double &sigmaK, double &sigmaW) const override {return 0.0;} - double WallBeta() const override {return 1.0;} + double EddyViscNoLim(const primVars &state) const override {return 0.0;} + double InviscidSpecRad(const primVars &state, + const unitVec3dMag &fAreaL, + const unitVec3dMag &fAreaR) const override { + return 0.0; + } + double TkeMin() const override {return 0.0;} double OmegaMin() const override {return 0.0;} @@ -220,9 +210,6 @@ class turbKWWilcox : public turbModel { double SrcSpecRad(const primVars &, const sutherland &, const double &) const override; - double InviscidSpecRad(const primVars &, - const unitVec3dMag &, - const unitVec3dMag &) const override; double ViscSpecRad(const primVars &, const unitVec3dMag &, const unitVec3dMag &, @@ -240,6 +227,9 @@ class turbKWWilcox : public turbModel { double Beta0() const {return beta0_;} double CLim() const {return clim_;} + double SigmaK() const override {return this->SigmaStar();} + double SigmaW() const override {return this->Sigma();} + void Print() const override; // destructor @@ -300,9 +290,6 @@ class turbKWSst : public turbModel { double SrcSpecRad(const primVars &, const sutherland &, const double &) const override; - double InviscidSpecRad(const primVars &, - const unitVec3dMag &, - const unitVec3dMag &) const override; double ViscSpecRad(const primVars &, const unitVec3dMag &, const unitVec3dMag &, @@ -325,6 +312,12 @@ class turbKWSst : public turbModel { double A1() const {return a1_;} double BetaStar() const {return betaStar_;} + // use coefficients from 1 because they are smaller + // this is used for TSL flux jacobian, so smaller will help increase + // diagonal dominance of implicit operator + double SigmaK() const override {return this->SigmaK1();} + double SigmaW() const override {return this->SigmaW1();} + void Print() const override; // destructor From 1a850121e9dd43b3b764b8af1ecad95450493710 Mon Sep 17 00:00:00 2001 From: Michael Date: Wed, 27 Apr 2016 23:52:29 -0700 Subject: [PATCH 5/5] Simplified building of k-d tree by using std::partition instead of home-cooked partitioning code --- kdtree.cpp | 37 ++++++++----------------------------- 1 file changed, 8 insertions(+), 29 deletions(-) diff --git a/kdtree.cpp b/kdtree.cpp index a7299eb..e343159 100644 --- a/kdtree.cpp +++ b/kdtree.cpp @@ -15,8 +15,7 @@ along with this program. If not, see . */ #include // cout -#include // exit() -#include // sort +#include // nth_element #include // numeric_limits #include // 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 &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