diff --git a/CHANGELOG.md b/CHANGELOG.md index e0e9cce6cba2786da73cfcec11a43652577e286d..4cc3bc8fdb6dee02512a92116363d100f071e53f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,7 @@ v5.0.0 - Add input module for Semtex field files (!777) - Fixed interppoints module (!760) - Move StreamFunction utility to a FieldConvert module (!809) +- Extend wss module to compressible flows (!810) - Allow explicitly setting bool options of FieldConvert modules as false (!811) - Enable output to multiple files (!844) - Allow using xml file without expansion tag in FieldConvert (!849) diff --git a/library/FieldUtils/ProcessModules/ProcessCreateExp.cpp b/library/FieldUtils/ProcessModules/ProcessCreateExp.cpp index 99f590ff3f0e9fde03d663ed04ceba98445c71c5..9624a0a2f8f7bca4e816415f3e1fe7b31dd73e97 100644 --- a/library/FieldUtils/ProcessModules/ProcessCreateExp.cpp +++ b/library/FieldUtils/ProcessModules/ProcessCreateExp.cpp @@ -193,13 +193,12 @@ void ProcessCreateExp::Process(po::variables_map &vm) m_f->m_session->LoadParameter("Strip_Z", nstrips, 1); - vector vars; + vector vars = m_f->m_session->GetVariables(); if (vm.count("useSessionVariables")) { - m_f->m_variables = m_f->m_session->GetVariables(); + m_f->m_variables = vars; } nfields = m_f->m_variables.size(); - vars = m_f->m_variables; m_f->m_exp.resize(nfields * nstrips); diff --git a/library/FieldUtils/ProcessModules/ProcessWSS.cpp b/library/FieldUtils/ProcessModules/ProcessWSS.cpp index 9c5890ad47cbdb8c5bf14f80ac579dee4b6c406a..0db43c26d336849cbf4a53ff0c88e809bdb03dd1 100644 --- a/library/FieldUtils/ProcessModules/ProcessWSS.cpp +++ b/library/FieldUtils/ProcessModules/ProcessWSS.cpp @@ -35,13 +35,14 @@ #include #include -using namespace std; #include "ProcessWSS.h" #include #include +using namespace std; + namespace Nektar { namespace FieldUtils @@ -64,17 +65,12 @@ void ProcessWSS::Process(po::variables_map &vm) { ProcessBoundaryExtract::Process(vm); - NekDouble kinvis = m_f->m_session->GetParameter("Kinvis"); - int i, j; int nfields = m_f->m_variables.size(); - int spacedim = m_f->m_graph->GetSpaceDimension(); - if ((m_f->m_numHomogeneousDir) == 1 || (m_f->m_numHomogeneousDir) == 2) - { - spacedim += m_f->m_numHomogeneousDir; - } + int expdim = m_f->m_graph->GetSpaceDimension(); + m_spacedim = expdim + m_f->m_numHomogeneousDir; - if (spacedim == 2) + if (m_spacedim == 2) { m_f->m_variables.push_back("Shear_x"); m_f->m_variables.push_back("Shear_y"); @@ -87,43 +83,42 @@ void ProcessWSS::Process(po::variables_map &vm) m_f->m_variables.push_back("Shear_z"); m_f->m_variables.push_back("Shear_mag"); } + if (m_f->m_exp[0]->GetNumElmts() == 0) { return; } - ASSERTL0(m_f->m_variables[0] == "u", - "Implicit assumption that input is in incompressible format of " - "(u,v,p) or (u,v,w,p)"); - - if (spacedim == 1) + if (m_spacedim == 1) { ASSERTL0(false, "Error: wss for a 1D problem cannot " "be computed"); } - int newfields = spacedim + 1; - int nshear = spacedim + 1; - int nstress = spacedim * spacedim; - int ngrad = spacedim * spacedim; + // Load viscosity coefficients + NekDouble kinvis, lambda; + GetViscosity( kinvis, lambda); + + // Declare arrays + int nshear = m_spacedim + 1; + int nstress = m_spacedim * m_spacedim; + int ngrad = m_spacedim * m_spacedim; - Array > velocity(nfields), grad(ngrad), - fgrad(ngrad); + Array > velocity(nfields); + Array > grad(ngrad); Array > stress(nstress), fstress(nstress); Array > fshear(nshear); - Array BndExp(newfields); - Array BndElmtExp(spacedim); + Array BndExp(nshear); + Array BndElmtExp(nfields); - m_f->m_exp.resize(nfields + newfields); - for (i = 0; i < newfields; ++i) + // Resize m_exp + m_f->m_exp.resize(nfields + nshear); + for (i = 0; i < nshear; ++i) { - m_f->m_exp[nfields + i] = - m_f->AppendExpList(m_f->m_numHomogeneousDir); + m_f->m_exp[nfields + i] = m_f->AppendExpList(m_f->m_numHomogeneousDir); } - - // Create map of boundary ids for partitioned domains SpatialDomains::BoundaryConditions bcs(m_f->m_session, m_f->m_exp[0]->GetGraph()); @@ -135,6 +130,7 @@ void ProcessWSS::Process(po::variables_map &vm) { BndRegionMap[breg_it.first] = cnt++; } + // Loop over boundaries to Write for (int b = 0; b < m_f->m_bndRegionsToWrite.size(); ++b) { @@ -143,12 +139,12 @@ void ProcessWSS::Process(po::variables_map &vm) int bnd = BndRegionMap[m_f->m_bndRegionsToWrite[b]]; // Get expansion list for boundary and for elements containing this // bnd - for (i = 0; i < newfields; i++) + for (i = 0; i < nshear; i++) { BndExp[i] = m_f->m_exp[nfields + i]->UpdateBndCondExpansion(bnd); } - for (i = 0; i < spacedim; i++) + for (i = 0; i < nfields; i++) { m_f->m_exp[i]->GetBndElmtExpansion(bnd, BndElmtExp[i]); } @@ -171,56 +167,72 @@ void ProcessWSS::Process(po::variables_map &vm) stress[i] = Array(nqe); } + Array div(nqe, 0.0); + // initialise arrays in the boundary for (i = 0; i < nstress; ++i) { fstress[i] = Array(nqb); } - for (i = 0; i < ngrad; ++i) - { - fgrad[i] = Array(nqb); - } - for (i = 0; i < nshear; ++i) { fshear[i] = Array(nqb, 0.0); } // Extract Velocities - for (i = 0; i < spacedim; ++i) - { - velocity[i] = BndElmtExp[i]->GetPhys(); - } + GetVelocity( BndElmtExp, velocity); - // Compute gradients (velocity correction scheme method) - for (i = 0; i < spacedim; ++i) + // Compute gradients + for (i = 0; i < m_spacedim; ++i) { - if (spacedim == 2) + if (m_spacedim == 2) { BndElmtExp[i]->PhysDeriv(velocity[i], - grad[i * spacedim + 0], - grad[i * spacedim + 1]); + grad[i * m_spacedim + 0], + grad[i * m_spacedim + 1]); } else { - BndElmtExp[i]->PhysDeriv( - velocity[i], grad[i * spacedim + 0], - grad[i * spacedim + 1], grad[i * spacedim + 2]); + BndElmtExp[i]->PhysDeriv(velocity[i], + grad[i * m_spacedim + 0], + grad[i * m_spacedim + 1], + grad[i * m_spacedim + 2]); } + // Add contribution to div(u) + Vmath::Vadd(nqe, grad[i * m_spacedim + i], 1, div, 1, div, 1); } - // Compute stress component terms tau_ij = mu*(u_i,j + u_j,i) - for (i = 0; i < spacedim; ++i) + // Velocity divergence scaled by lambda * mu + Vmath::Smul(nqe, lambda*kinvis, div, 1, div, 1); + + // Compute stress component terms + // tau_ij = mu*(u_i,j + u_j,i) + mu*lambda*delta_ij*div(u) + for (i = 0; i < m_spacedim; ++i) { - for (j = 0; j < spacedim; ++j) + for (j = i; j < m_spacedim; ++j) { - Vmath::Vadd(nqe, grad[i * spacedim + j], 1, - grad[j * spacedim + i], 1, - stress[i * spacedim + j], 1); - - Vmath::Smul(nqe, kinvis, stress[i * spacedim + j], 1, - stress[i * spacedim + j], 1); + Vmath::Vadd(nqe, grad[i * m_spacedim + j], 1, + grad[j * m_spacedim + i], 1, + stress[i * m_spacedim + j], 1); + + Vmath::Smul(nqe, kinvis, + stress[i * m_spacedim + j], 1, + stress[i * m_spacedim + j], 1); + + if (i == j) + { + // Add divergence term to diagonal + Vmath::Vadd(nqe, stress[i * m_spacedim + j], 1, + div, 1, + stress[i * m_spacedim + j], 1); + } + else + { + // Copy to make symmetric + Vmath::Vcopy(nqe, stress[i * m_spacedim + j], 1, + stress[j * m_spacedim + i], 1); + } } } @@ -234,32 +246,33 @@ void ProcessWSS::Process(po::variables_map &vm) Array > normals; m_f->m_exp[0]->GetBoundaryNormals(bnd, normals); // Reverse normals, to get correct orientation for the body - for (i = 0; i < spacedim; ++i) + for (i = 0; i < m_spacedim; ++i) { Vmath::Neg(nqb, normals[i], 1); } // calculate wss, and update coeffs in the boundary expansion // S = tau_ij * n_j - for (i = 0; i < spacedim; ++i) + for (i = 0; i < m_spacedim; ++i) { - for (j = 0; j < spacedim; ++j) + for (j = 0; j < m_spacedim; ++j) { - Vmath::Vvtvp(nqb, normals[j], 1, fstress[i * spacedim + j], - 1, fshear[i], 1, fshear[i], 1); + Vmath::Vvtvp(nqb, normals[j], 1, + fstress[i * m_spacedim + j], 1, + fshear[i], 1, fshear[i], 1); } } // T = S - (S.n)n - for (i = 0; i < spacedim; ++i) + for (i = 0; i < m_spacedim; ++i) { Vmath::Vvtvp(nqb, normals[i], 1, fshear[i], 1, fshear[nshear - 1], 1, fshear[nshear - 1], 1); } - Vmath::Smul(nqb, -1.0, fshear[nshear - 1], 1, fshear[nshear - 1], - 1); + Vmath::Smul(nqb, -1.0, fshear[nshear - 1], 1, + fshear[nshear - 1], 1); - for (i = 0; i < spacedim; i++) + for (i = 0; i < m_spacedim; i++) { Vmath::Vvtvp(nqb, normals[i], 1, fshear[nshear - 1], 1, fshear[i], 1, fshear[i], 1); @@ -269,7 +282,7 @@ void ProcessWSS::Process(po::variables_map &vm) // Tw Vmath::Zero(nqb, fshear[nshear - 1], 1); - for (i = 0; i < spacedim; ++i) + for (i = 0; i < m_spacedim; ++i) { Vmath::Vvtvp(nqb, fshear[i], 1, fshear[i], 1, fshear[nshear - 1], 1, fshear[nshear - 1], 1); @@ -280,5 +293,64 @@ void ProcessWSS::Process(po::variables_map &vm) } } } + +void ProcessWSS::GetViscosity(NekDouble &kinvis, NekDouble &lambda) +{ + if(boost::iequals(m_f->m_variables[0], "u")) + { + // IncNavierStokesSolver + kinvis = m_f->m_session->GetParameter("Kinvis"); + lambda = 0; + } + else if(boost::iequals(m_f->m_variables[0], "rho") && + boost::iequals(m_f->m_variables[1], "rhou")) + { + // CompressibleFlowSolver + m_f->m_session->LoadParameter ("mu", kinvis, 1.78e-05); + m_f->m_session->LoadParameter ("lambda", lambda, -2.0/3.0); + } + else + { + // Unknown + ASSERTL0(false, "Invalid variables for WSS"); + } +} + +void ProcessWSS::GetVelocity( + const Array exp, + Array > &vel) +{ + int npoints = exp[0]->GetNpoints(); + if(boost::iequals(m_f->m_variables[0], "u")) + { + // IncNavierStokesSolver + for (int i = 0; i < m_spacedim; ++i) + { + vel[i] = Array(npoints); + Vmath::Vcopy(npoints, + exp[i]->GetPhys(), 1, + vel[i], 1); + } + } + else if(boost::iequals(m_f->m_variables[0], "rho") && + boost::iequals(m_f->m_variables[1], "rhou")) + { + // CompressibleFlowSolver + for (int i = 0; i < m_spacedim; ++i) + { + vel[i] = Array(npoints); + Vmath::Vdiv(npoints, + exp[i + 1]->GetPhys(), 1, + exp[0]->GetPhys(), 1, + vel[i], 1); + } + } + else + { + // Unknown + ASSERTL0(false, "Could not identify velocity for WSS"); + } +} + } } diff --git a/library/FieldUtils/ProcessModules/ProcessWSS.h b/library/FieldUtils/ProcessModules/ProcessWSS.h index f2ecd6b975a72aa4c91563602ef1d4054be3bd4a..eb34281a57fb93ccab1546128ca6a4a3d2b695d3 100644 --- a/library/FieldUtils/ProcessModules/ProcessWSS.h +++ b/library/FieldUtils/ProcessModules/ProcessWSS.h @@ -73,6 +73,15 @@ public: return "Calculating wall shear stress"; } +protected: + void GetViscosity(NekDouble &kinvis, NekDouble &lambda); + + void GetVelocity(const Array exp, + Array > &vel); + +private: + int m_spacedim; + }; } }