...

Commits (44)
 ... ... @@ -444,8 +444,7 @@ Under the two following sections it is possible to define the initial conditions \section{Examples} \subsection{Shock capturing} Compressible flows can be characterised by abrupt changes in density within the flow domain often referred to as shocks. These discontinuities lead to numerical instabilities (Gibbs phenomena). This problem is prevented by locally adding a diffusion term to the equations to damp the numerical oscillations. An artificial diffusion term is introduced locally to the Euler equations to deal with flow discontinuity and the consequential numerical oscillations. Two models are implemented, a non-smooth and a smooth artificial viscosity model. Compressible flows can be characterised by abrupt changes in density within the flow domain often referred to as shocks. These discontinuities can lead to numerical instabilities (Gibbs phenomena). This problem is prevented by locally adding a diffusion term to the equations to damp the numerical oscillations. \subsubsection{Non-smooth artificial viscosity model} For the non-smooth artificial viscosity model the added artificial viscosity is constant in each element and discontinuous between the elements. The Euler system is augmented by an added Laplacian term on right hand side of equation \ref{eq:euler} \cite{persson2006sub}. ... ... @@ -480,9 +479,9 @@ where $s_0 = s_\kappa - 4.25\;log_{10}(p)$. To enable the non-smooth viscosity model, the following line has to be added to the \inltt{SOLVERINFO} section: \begin{lstlisting}[style=XmlStyle] \end{lstlisting} The diffusivity and the sensor can be controlled by the following parameters: \begin{lstlisting}[style=XmlStyle] ... ... @@ -502,45 +501,6 @@ The diffusivity and the sensor can be controlled by the following parameters: \end{center} \end{figure} \subsubsection{Smooth artificial viscosity model} For the smooth artificial viscosity model an extra PDE for the artificial viscosity is appended to the Euler system \label{eq:eulerplusvis}\begin{split} \frac{\partial \epsilon}{\partial t} &= \nabla\cdot \left(\nabla \epsilon\right) + \frac{1}{\tau}\left(\frac{h}{p}\lambda_{max}S_\kappa - \epsilon\right)\ \ \ \ \ \ \textrm{on}\ \ \ \ \ \ \ \Omega\\ \frac{\partial \epsilon}{\partial n} &= 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \textrm{on}\ \ \ \ \ \ \ \Gamma \end{split} where $S_\kappa$ is a normalised sensor value and serves as a forcing term for the artificial viscosity. A smooth artificial viscosity distribution is obtained.\\ \\ To enable the smooth viscosity model, the following line has to be added to the \inltt{SOLVERINFO} section: \begin{lstlisting}[style=XmlStyle] \end{lstlisting} Furthermore, the extra viscosity variable \inltt{eps} has to be added to the variable list: \begin{lstlisting}[style=XmlStyle] rho rhou rhov E eps \end{lstlisting} A similar addition has to be made for the boundary conditions and initial conditions. The tests that have been run started with a uniform homogeneous boundary condition and initial condition. The following parameters can be set in the xml session file: \begin{lstlisting}[style=XmlStyle]

Skappa = -1.3

Kappa = 0.2

mu0 = 1.0

FH = 3

FL = 0.01*FH

C1 = 0.03

C2 = 5/3*C1

\end{lstlisting} where for now \inltt{FH} and \inltt{FL} are used to tune which range of the sensor is used as a forcing term and \inltt{C1} and \inltt{C2} are fixed constants which can be played around with to make the model more diffusive or not. However these constants are generally fixed. \subsection{Variable polynomial order} A sensor based $p$-adaptive algorithm is implemented to optimise the computational cost and accuracy. The DG scheme allows one to use different polynomial orders since the fluxes over the elements are determined using a Riemann solver and there is now further coupling between the elements. Furthermore, the initial $p$-adaptive algorithm uses the same sensor as the shock capturing algorithm to identify the smoothness of the local solution so it rather straightforward to implement both algorithms at the same time.\\ ... ...
 ... ... @@ -1258,6 +1258,12 @@ namespace Nektar id = m_traceMap->GetBndCondCoeffsToGlobalCoeffsMap(i); BndRhs[id] += m_bndCondExpansions[i]->GetCoeff(0); } else if (m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::ePeriodic) { ASSERTL0(false, "HDG implementation does not support " "periodic boundary conditions at present."); } } //---------------------------------- ... ...
 ... ... @@ -1957,6 +1957,12 @@ namespace Nektar BndRhs[id] += m_bndCondExpansions[i]->GetCoeffs()[j]; } } else if (m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::ePeriodic) { ASSERTL0(false, "HDG implementation does not support " "periodic boundary conditions at present."); } } //---------------------------------- ... ...
 ... ... @@ -2332,6 +2332,12 @@ using namespace std; BndRhs[id] += m_bndCondExpansions[i]->GetCoeffs()[j]; } } else if (m_bndConditions[i]->GetBoundaryConditionType() == SpatialDomains::ePeriodic) { ASSERTL0(false, "HDG implementation does not support " "periodic boundary conditions at present."); } } //---------------------------------- ... ...
 ... ... @@ -130,11 +130,6 @@ namespace Nektar int numConvFields = nConvectiveFields; if (m_shockCaptureType == "Smooth") { numConvFields = nConvectiveFields - 1; } for (j = 0; j < nDim; ++j) { for (i = 0; i < numConvFields; ++i) ... ...
 ... ... @@ -91,7 +91,7 @@ namespace Nektar m_session->MatchSolverInfo("REACTIONADVANCEMENT", "Explicit", m_explicitReaction, true); m_session->LoadParameter("CheckNanSteps", m_nanSteps, 1); m_session->LoadParameter("CheckAbortSteps", m_abortSteps, 1); // Steady state tolerance m_session->LoadParameter("SteadyStateTol", m_steadyStateTol, 0.0); // Frequency for checking steady state ... ... @@ -263,20 +263,28 @@ namespace Nektar NekDouble elapsed = 0.0; NekDouble totFilterTime = 0.0; while (step < m_steps || m_time < m_fintime - NekConstants::kNekZeroTol) Array abortFlags(2, 0); string abortFile = "abort"; if (m_session->DefinesSolverInfo("CheckAbortFile")) { abortFile = m_session->GetSolverInfo("CheckAbortFile"); } while ((step < m_steps || m_time < m_fintime - NekConstants::kNekZeroTol) && abortFlags[1] == 0) { if (m_cflSafetyFactor) { m_timestep = GetTimeStep(fields); // Ensure that the final timestep finishes at the final // time, or at a prescribed IO_CheckTime. if (m_time + m_timestep > m_fintime && m_fintime > 0.0) { m_timestep = m_fintime - m_time; } else if (m_checktime && else if (m_checktime && m_time + m_timestep - lastCheckTime >= m_checktime) { lastCheckTime += m_checktime; ... ... @@ -284,7 +292,7 @@ namespace Nektar doCheckTime = true; } } // Perform any solver-specific pre-integration steps timer.Start(); if (v_PreIntegrate(step)) ... ... @@ -358,21 +366,34 @@ namespace Nektar } } // search for NaN and quit if found if (m_nanSteps && !((step+1) % m_nanSteps) ) // test for abort conditions (nan, or abort file) if (m_abortSteps && !((step+1) % m_abortSteps) ) { int nanFound = 0; abortFlags[0] = 0; for (i = 0; i < nvariables; ++i) { if (Vmath::Nnan(fields[i].num_elements(), fields[i], 1) > 0) { nanFound = 1; abortFlags[0] = 1; } } m_session->GetComm()->AllReduce(nanFound, //rank zero looks for abort file and deltes it //if it exists. The communicates the abort if(m_session->GetComm()->GetRank() == 0) { if(boost::filesystem::exists(abortFile)) { boost::filesystem::remove(abortFile); abortFlags[1] = 1; } } m_session->GetComm()->AllReduce(abortFlags, LibUtilities::ReduceMax); ASSERTL0 (!nanFound, ASSERTL0 (!abortFlags[0], "NaN found during time integration."); } ... ... @@ -458,7 +479,7 @@ namespace Nektar ++step; ++stepCounter; } // Print out summary statistics if (m_session->GetComm()->GetRank() == 0) { ... ... @@ -473,7 +494,7 @@ namespace Nektar cout << "Time-integration : " << intTime << "s" << endl; } } // If homogeneous, transform back into physical space if necessary. if(m_HomogeneousType != eNotHomogeneous) { ... ... @@ -502,7 +523,7 @@ namespace Nektar { x.second->Finalise(m_fields, m_time); } // Print for 1D problems if(m_spacedim == 1) { ... ...
 ... ... @@ -60,6 +60,8 @@ public: protected: /// Number of time steps between outputting status information. int m_infosteps; /// Number of steps between checks for abort conditions. int m_abortSteps; /// Number of time steps between outputting filters information. int m_filtersInfosteps; int m_nanSteps; ... ...
 ... ... @@ -209,14 +209,14 @@ C[1-8] ... ...
 3D unsteady DG advection, hexahedra, order 1, P=12,periodic bcs 3D unsteady DG advection, hexahedra, order 1, P=10,periodic bcs ADRSolver --use-scotch Advection3D_m12_DG_hex_VarP.xml --use-scotch Advection3D_m10_DG_hex_VarP.xml 3 Advection3D_m12_DG_hex_VarP.xml Advection3D_m10_DG_hex_VarP.xml 3.99686e-06 5.67026e-06 1.50105e-05 1.94197e-05 \ No newline at end of file
 ... ... @@ -51,7 +51,7 @@ C[0] ... ...
 ... ... @@ -2,16 +2,16 @@ 3D unsteady DG advection, tetrahedra, order 4, P=Variable ADRSolver Advection3D_m12_DG_prism_VarP.xml Advection3D_m10_DG_prism_VarP.xml Advection3D_m12_DG_prism_VarP.xml Advection3D_m10_DG_prism_VarP.xml 1.74811e-05 1.5189e-07 0.00323074 7.42306e-06 ... ...
 ... ... @@ -498,15 +498,15 @@ C[0-8] ... ...
 ... ... @@ -2,16 +2,16 @@ 3D unsteady DG advection, prisms, order 4, P=Variable, HDF output ADRSolver --io-format Hdf5 Advection3D_m12_DG_prism_VarP.xml --io-format Hdf5 Advection3D_m10_DG_prism_VarP.xml Advection3D_m12_DG_prism_VarP.xml Advection3D_m10_DG_prism_VarP.xml 1.74811e-05 1.5189e-07 0.00323074 7.42306e-06 ... ...
 ... ... @@ -2,17 +2,17 @@ 3D unsteady DG advection, tetrahedra, order 4, P=Variable ADRSolver --use-scotch Advection3D_m12_DG_prism_VarP.xml --use-scotch Advection3D_m10_DG_prism_VarP.xml 3 Advection3D_m12_DG_prism_VarP.xml Advection3D_m10_DG_prism_VarP.xml 1.74811e-05 1.5189e-07 0.00323074 7.42306e-06 ... ...
 ... ... @@ -2,16 +2,16 @@ 3D unsteady DG advection, hexahedra, order 4, P=12 ADRSolver Advection3D_m12_DG_hex.xml Advection3D_m8_DG_hex.xml Advection3D_m12_DG_hex.xml Advection3D_m8_DG_hex.xml 3.1094e-8 1.08584e-05 1.3031e-6 2.52297e-05 ... ...
 ... ... @@ -126,7 +126,7 @@ C[0] ... ...
 ... ... @@ -126,7 +126,7 @@ C[0] ... ...
 ... ... @@ -2,17 +2,17 @@ 3D unsteady DG advection, hexahedra, order 4, P=12, par(3) ADRSolver --use-scotch Advection3D_m12_DG_hex.xml --use-scotch Advection3D_m8_DG_hex.xml 3 Advection3D_m12_DG_hex.xml Advection3D_m8_DG_hex.xml 3.1094e-8 1.08584e-05 1.3031e-6 2.52297e-05 ... ...
 ... ... @@ -9,7 +9,7 @@ 0.000129406 0.000129406 0.0461794 ... ...
 ... ... @@ -8,10 +8,10 @@ 0.0337347 0.0371179 0.205011 0.115599
 ... ... @@ -98,7 +98,7 @@ C[1] ... ...
 ... ... @@ -8,10 +8,10 @@ 0.0116705 0.0264016 0.177522 0.0861253
 ... ... @@ -540,7 +540,7 @@ C[1] ... ...
 /////////////////////////////////////////////////////////////////////////////// // // File: SmoothShockCapture.h // // For more information, please see: http://www.nektar.info // // The MIT License // // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA), // Department of Aeronautics, Imperial College London (UK), and Scientific // Computing and Imaging Institute, University of Utah (USA). // // Permission is hereby granted, free of charge, to any person obtaining a // copy of this software and associated documentation files (the "Software"), // to deal in the Software without restriction, including without limitation // the rights to use, copy, modify, merge, publish, distribute, sublicense, // and/or sell copies of the Software, and to permit persons to whom the // Software is furnished to do so, subject to the following conditions: // // The above copyright notice and this permission notice shall be included // in all copies or substantial portions of the Software. // // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER // DEALINGS IN THE SOFTWARE. // // Description: Smooth artificial diffusion for shock capture // /////////////////////////////////////////////////////////////////////////////// #ifndef NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_ARTIFICIALDIFFUSION_SMOOTH #define NEKTAR_SOLVERS_COMPRESSIBLEFLOWSOLVER_ARTIFICIALDIFFUSION_SMOOTH #include "ArtificialDiffusion.h" namespace Nektar { /** * @brief Smooth artificial diffusion for shock capture for compressible flow * problems. */ class SmoothShockCapture : public ArtificialDiffusion { public: friend class MemoryManager; /// Creates an instance of this class static ArtificialDiffusionSharedPtr create( const LibUtilities::SessionReaderSharedPtr& pSession, const Array& pFields, const int spacedim) { ArtificialDiffusionSharedPtr p = MemoryManager:: AllocateSharedPtr(pSession, pFields, spacedim); return p; } ///Name of the class static std::string className; protected: virtual void v_DoArtificialDiffusion( const Array > &inarray, Array > &outarray); virtual void v_GetArtificialViscosity( const Array > &physfield, Array &mu); private: SmoothShockCapture(const LibUtilities::SessionReaderSharedPtr& pSession, const Array& pFields, const int spacedim); virtual ~SmoothShockCapture(void){}; void GetForcingTerm( const Array > &inarray, Array > outarrayForcing); /// Parameters NekDouble m_FacL; NekDouble m_FacH; NekDouble m_hFactor; NekDouble m_C1; NekDouble m_C2; NekDouble m_mu0; int m_offset; }; } #endif
 ... ... @@ -10,7 +10,6 @@ IF( NEKTAR_SOLVER_COMPRESSIBLE_FLOW ) SET(CompressibleFlowSolverSource ./ArtificialDiffusion/ArtificialDiffusion.cpp ./ArtificialDiffusion/NonSmoothShockCapture.cpp ./ArtificialDiffusion/SmoothShockCapture.cpp ./BoundaryConditions/CFSBndCond.cpp ./BoundaryConditions/ExtrapOrder0BC.cpp ./BoundaryConditions/IsentropicVortexBC.cpp ... ...
 ... ... @@ -49,27 +49,29 @@

TimeStep = 0.00001

NumSteps = 0

FinTime = 0.38

IO_CheckSteps = 10000

IO_InfoSteps = 100

Gamma = 1.4

IO_CheckTime = 100

FinTime = 1

IO_CheckTime = 0.1

IO_InfoSteps = 10

Gamma = 1.4

pInf = 101325

rhoInf = 1.225

uInf = 0.1

vInf = 0.0

CFL = 0.5

... ... @@ -92,13 +94,6 @@
... ...
 ... ... @@ -247,7 +247,7 @@ ... ... @@ -263,8 +263,8 @@

TimeStep = 0.002

NumSteps = 70

TimeStep = 0.007

NumSteps = 20

IO_CheckSteps = 1000

IO_InfoSteps = 20

Re = 7500

... ...
 ... ... @@ -253,13 +253,13 @@

TimeStep = 0.002

NumSteps = 70

TimeStep = 0.007

NumSteps = 20

IO_CheckSteps = 1000

IO_InfoSteps = 20

Re = 7500

Kinvis = 1.0/Re

kdim = 16

kdim = 16

nvec =2

evtol =1e-6

... ...
 ... ... @@ -14,7 +14,7 @@

TimeStep = 0.01

NumSteps = 1.0/TimeStep

NumSteps = 0.5/TimeStep

IO_CheckSteps = 25.0/TimeStep

IO_InfoSteps = 0.5/TimeStep

IO_CFLSteps = IO_InfoSteps

... ...
 3D channel flow, Hex elements, par(2), P=8 3D channel flow, Hex elements, par(2), P=3 IncNavierStokesSolver --use-scotch Hex_channel_m8_par.xml --use-scotch Hex_channel_m3_par.xml 2 Hex_channel_m8_par.xml Hex_channel_m3_par.xml ... ...
 ... ... @@ -146,7 +146,7 @@ C[0] ... ...
 3D channel flow, Hexahedral elements, P=8, Successive RHS(5) 3D channel flow, Hexahedral elements, P=3, Successive RHS(5) IncNavierStokesSolver Hex_channel_m8_srhs.xml Hex_channel_m3_srhs.xml Hex_channel_m8_srhs.xml Hex_channel_m3_srhs.xml ... ...
 ... ... @@ -146,7 +146,7 @@ C[0]