Commit 6cb5ebd5 by Giacomo Castiglioni Committed by Dave Moxey

### Fix/ldg penalty

parent 36b04d7e
 ... ... @@ -14,3 +14,6 @@ ThirdParty *.pdf # Kdevelop project files *.kdev4 # Sublime project files *.sublime-project *.sublime-workspace
 ... ... @@ -132,6 +132,7 @@ v5.0.0 - Fix surface extraction, added regression test (!994) - Fix 2D meshing running out of memory due to missing else (!1012) - Add support for .msh v4.1 file input (!1054) - Added penalty term to LDG and LDGNS, slight generalization of LDG (!1080) **FieldConvert**: - Add input module for Semtex field files (!777) ... ...
 ... ... @@ -11,14 +11,14 @@ The ADRSolver is designed to solve partial differential equations of the form: \alpha \dfrac{\partial u}{\partial t} + \lambda u + \nu \nabla u + \epsilon \nabla \cdot (D \nabla u) = f in either discontinuous or continuous projections of the solution field. For a full list of the equations which are supported, and the capabilities of each equation, in either discontinuous or continuous projections of the solution field. For a full list of the equations which are supported, and the capabilities of each equation, see the table below. \begin{table}[h!] \begin{center} \tiny \renewcommand\arraystretch{2.2} \renewcommand\arraystretch{2.2} \begin{tabular}{llll} \toprule \textbf{Equation to solve} & \textbf{EquationType} & \textbf{Dimensions} & ... ... @@ -26,22 +26,22 @@ see the table below. \midrule $u = f$ & \inltt{Projection} & All & Continuous/Discontinuous \\ $\nabla^2 u = 0$ & $\nabla^2 u = 0$ & \inltt{Laplace} & All & Continuous/Discontinuous \\ $\nabla^2 u = f$ & $\nabla^2 u = f$ & \inltt{Poisson} & All & Continuous/Discontinuous \\ $\nabla^2 u + \lambda u = f$ & $\nabla^2 u + \lambda u = f$ & \inltt{Helmholtz} & All & Continuous/Discontinuous \\ $\epsilon \nabla^2 u + \mathbf{V}\nabla u = f$ & $\epsilon \nabla^2 u + \mathbf{V}\nabla u = f$ & \inltt{SteadyAdvectionDiffusion} & 2D only & Continuous/Discontinuous \\ $\epsilon \nabla^2 u + \lambda u = f$ & $\epsilon \nabla^2 u + \lambda u = f$ & \inltt{SteadyDiffusionReaction} & 2D only & Continuous/Discontinuous \\ $\epsilon \nabla^2 u + \mathbf{V}\nabla u + \lambda u = f$ & \inltt{SteadyAdvectionDiffusionReaction} & 2D only & $\epsilon \nabla^2 u + \mathbf{V}\nabla u + \lambda u = f$ & \inltt{SteadyAdvectionDiffusionReaction} & 2D only & Continuous/Discontinuous \\ $\dfrac{\partial u}{\partial t} + \mathbf{V}\nabla u = f$ & \inltt{UnsteadyAdvection} & All & Continuous/Discontinuous \\ $\dfrac{\partial u}{\partial t} = \epsilon \nabla^2 u$ & $\dfrac{\partial u}{\partial t} = \epsilon \nabla^2 u$ & \inltt{UnsteadyDiffusion} & All & Continuous/Discontinuous \\ $\dfrac{\partial u}{\partial t} = \epsilon \nabla^2 u + R(u)$ & \inltt{UnsteadyReactionDiffusion} & All & Continuous \\ ... ... @@ -64,7 +64,7 @@ ADRSolver session.xml \section{Session file configuration} The type of equation which is to be solved is specified through the EquationType The type of equation which is to be solved is specified through the EquationType SOLVERINFO option in the session file. This can be set as in table \ref{t:ADR1}. At present, the Steady non-symmetric solvers cannot be used in parallel. \\ ... ... @@ -76,7 +76,7 @@ The solver info are listed below: \item \textbf{TimeIntegrationMethod}: The following types of time integration methods have been tested with each solver: \begin{center} \footnotesize \renewcommand\arraystretch{1.2} \renewcommand\arraystretch{1.2} \begin{tabular}{lcccc} \toprule \textbf{EqType} & \textbf{Explicit} & \textbf{Diagonally Implicit} & ... ... @@ -116,7 +116,7 @@ The solver info are listed below: \end{itemize} \item \textbf{DiffusionType}: \begin{itemize} \item \inltt{LDG}. \item \inltt{LDG} (The penalty term is proportional to an optional parameter \inltt{LDGc11} which is by default set to one; proportionality to polynomial order can be manually imposed by setting the parameter \inltt{LDGc11} equal to $p^2$). \end{itemize} \item \textbf{UpwindType}: \begin{itemize} ... ... @@ -129,13 +129,13 @@ The solver info are listed below: The following parameters can be specified in the \inltt{PARAMETERS} section of the session file: \begin{itemize} \item \inltt{epsilon}: sets the diffusion coefficient $\epsilon$.\\ \item \inltt{epsilon}: sets the diffusion coefficient $\epsilon$.\\ \textit{Can be used} in: SteadyDiffusionReaction, SteadyAdvectionDiffusionReaction, UnsteadyDiffusion, UnsteadyAdvectionDiffusion. \\ \textit{Default value}: 0. \item \inltt{d00}, \inltt{d11}, \inltt{d22}: sets the diagonal entries of the diffusion tensor $D$. \\ \textit{Can be used in}: UnsteadyDiffusion \\ \textit{Default value}: All set to 1 (i.e. identity matrix). \textit{Default value}: All set to 1 (i.e. identity matrix). \item \inltt{lambda}: sets the reaction coefficient $\lambda$. \\ \textit{Can be used in}: SteadyDiffusionReaction, Helmholtz, SteadyAdvectionDiffusionReaction\\ \textit{Default value}: 0. ... ... @@ -215,7 +215,7 @@ and then define the actual advection function as Two boundary regions are defined, one at each end of the domain, and periodicity is enforced \begin{lstlisting}[style=XMLStyle] \begin{lstlisting}[style=XMLStyle] C[1] C[2] ... ... @@ -268,7 +268,7 @@ We consider the elliptic partial differential equation: where $\nabla^2$ is the Laplacian and $\lambda$ is a real positive constant. \subsubsection{Input file} \subsubsection{Input file} The input for this example is given in the example file \inlsh{Helmholtz2D\_modal.xml} ... ... @@ -300,7 +300,7 @@ basis with $7$ modes (maximum polynomial order is $6$). \end{lstlisting} Only one parameter is needed for this problem. In this example $\lambda = 1$ and the Continuous Galerkin Method is used as projection scheme to solve the the Continuous Galerkin Method is used as projection scheme to solve the Helmholtz equation, so we need to specify the following parameters and solver information. \begin{lstlisting}[style=XMLStyle] ... ... @@ -331,7 +331,7 @@ region is then assigned an appropriate boundary condition. ... ... @@ -343,7 +343,7 @@ region is then assigned an appropriate boundary condition. \end{lstlisting} We know that for $f = -(\lambda + 2 \pi^2)sin(\pi x)cos(\pi y)$, the exact We know that for $f = -(\lambda + 2 \pi^2)sin(\pi x)cos(\pi y)$, the exact solution of the two-dimensional Helmholtz equation is $u = sin(\pi x)cos(\pi y)$. These functions are defined specified to initialise the problem and verify the correct solution is obtained by evaluating the $L_2$ and $L_{inf}$ errors. ... ... @@ -363,7 +363,7 @@ the correct solution is obtained by evaluating the $L_2$ and $L_{inf}$ errors. ADRSolver Test_Helmholtz2D_modal.xml \end{lstlisting} This execution should print out a summary of input file, the $L_2$ and This execution should print out a summary of input file, the $L_2$ and $L_{inf}$ errors and the time spent on the calculation. \subsubsection{Post-processing} ... ... @@ -417,9 +417,9 @@ Schmidt number (the relative thickness of the momentum to mass transfer boundary layer) is sufficiently large. The analytical solution for the non-dimensional mass transfer at the wall is given by: \begin{align*} S h(z) = \dfrac{2^{4/3}(Pe R/z)^{1/3}}{g^{1/3}\Gamma(4/3)} , S h(z) = \dfrac{2^{4/3}(Pe R/z)^{1/3}}{g^{1/3}\Gamma(4/3)} , \end{align*} where $z$ is the streamwise coordinate, $R$ the pipe radius, $\Gamma(4/3)$ an incomplete where $z$ is the streamwise coordinate, $R$ the pipe radius, $\Gamma(4/3)$ an incomplete Gamma function and $Pe$ the Peclet number given by: \begin{align*} Pe = \dfrac{2 U R}{\epsilon} ... ... @@ -538,7 +538,7 @@ flow), hence we can define this through an analytical function as follows: \end{lstlisting} We assume that the initial domain concentration is uniform everywhere and the same as the inlet. This is defined by, same as the inlet. This is defined by, \begin{lstlisting}[style=XMLStyle] ... ...
 ... ... @@ -272,7 +272,7 @@ Note that only \inltt{WeakDG} is fully supported, the other operators work only \item \inltt{DiffusionType} is the diffusion operator we want to use for the Navier-Stokes equations: \begin{itemize} \item \inltt{LDGNS} (LDG); \item \inltt{LDGNS} (LDG with primitive variables. The penalty term is inversely proportional to the element size, proportional to the local viscosity for the momentum equations and to the thermal conductivity for the energy equation, and proportional to an optional parameter \inltt{LDGNSc11} which is by default set to one; proportionality to polynomial order can be manually imposed by setting the parameter \inltt{LDGNSc11} equal to $p^2$); \item \inltt{LFRDGNS} (Flux-Reconstruction recovering nodal DG scheme); \item \inltt{LFRSDNS} (Flux-Reconstruction recovering a spectral difference (SD) scheme); \item \inltt{LFRHUNS} (Flux-Reconstruction recovering Huynh (G2) scheme); ... ...
 ... ... @@ -43,16 +43,16 @@ namespace Nektar static DiffusionFactory instance; return instance; } void Diffusion::InitObject( const LibUtilities::SessionReaderSharedPtr pSession, Array pFields) { v_InitObject(pSession, pFields); } void Diffusion::Diffuse( const int nConvectiveFields, const std::size_t nConvectiveFields, const Array &fields, const Array > &inarray, Array > &outarray, ... ...
 ... ... @@ -52,23 +52,28 @@ namespace Nektar namespace SolverUtils { typedef std::function >&, Array >&, Array >&)> DiffusionFluxVecCB; const Array > &, const Array > >&, Array > >&)> DiffusionFluxVecCB; typedef std::function >&, Array > >&, Array > >&)> Array > >&)> DiffusionFluxVecCBNS; typedef std::function >&, const Array >&, Array >&)> DiffusionFluxPenaltyNS; typedef std::function >&, Array&)> DiffusionArtificialDiffusion; class Diffusion { public: ... ... @@ -79,34 +84,33 @@ namespace Nektar SOLVER_UTILS_EXPORT void InitObject( LibUtilities::SessionReaderSharedPtr pSession, Array pFields); SOLVER_UTILS_EXPORT void Diffuse( const int nConvectiveFields, const std::size_t nConvectiveFields, const Array &fields, const Array > &inarray, Array > &outarray, const Array > &pFwd = NullNekDoubleArrayofArray, const Array > &pBwd = NullNekDoubleArrayofArray); SOLVER_UTILS_EXPORT void FluxVec( Array > > &fluxvector); template template void SetFluxVector(FuncPointerT func, ObjectPointerT obj) { m_fluxVector = std::bind( func, obj, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5); std::placeholders::_3); } void SetFluxVectorVec(DiffusionFluxVecCB fluxVector) void SetFluxVector(DiffusionFluxVecCB fluxVector) { m_fluxVector = fluxVector; } template template void SetFluxVectorNS(FuncPointerT func, ObjectPointerT obj) { m_fluxVectorNS = std::bind( ... ... @@ -114,16 +118,22 @@ namespace Nektar std::placeholders::_3); } void SetFluxVectorNS(DiffusionFluxVecCBNS fluxVector) { m_fluxVectorNS = fluxVector; } template void SetArtificialDiffusionVector(FuncPointerT func, ObjectPointerT obj) void SetFluxPenaltyNS(FuncPointerT func, ObjectPointerT obj) { m_ArtificialDiffusionVector = std::bind( func, obj, std::placeholders::_1, std::placeholders::_2); m_fluxPenaltyNS = std::bind( func, obj, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3); } void SetFluxVectorNS(DiffusionFluxVecCBNS fluxVector) void SetFluxPenaltyNS(DiffusionFluxPenaltyNS flux) { m_fluxVectorNS = fluxVector; m_fluxPenaltyNS = flux; } inline void SetHomoDerivs(Array > &deriv) ... ... @@ -135,11 +145,11 @@ namespace Nektar { return v_GetFluxTensor(); } protected: DiffusionFluxVecCB m_fluxVector; DiffusionFluxVecCBNS m_fluxVectorNS; DiffusionArtificialDiffusion m_ArtificialDiffusionVector; DiffusionFluxPenaltyNS m_fluxPenaltyNS; virtual void v_InitObject( LibUtilities::SessionReaderSharedPtr pSession, ... ... @@ -147,9 +157,9 @@ namespace Nektar { boost::ignore_unused(pSession, pFields); }; virtual void v_Diffuse( const int nConvectiveFields, const std::size_t nConvectiveFields, const Array &fields, const Array > &inarray, Array > &outarray, ... ... @@ -161,17 +171,17 @@ namespace Nektar { boost::ignore_unused(deriv); } virtual Array > > &v_GetFluxTensor() { static Array > > tmp; return tmp; } }; }; /// A shared pointer to an EquationSystem object typedef std::shared_ptr DiffusionSharedPtr; /// Datatype of the NekFactory used to instantiate classes derived /// from the Diffusion class. typedef LibUtilities::NekFactory DiffusionFactory; ... ...
 ... ... @@ -72,7 +72,7 @@ namespace Nektar }; /** * @brief Diffusion3DHomogeneous1D uses the 2D WeakDG approach * @brief Diffusion3DHomogeneous1D uses the 2D WeakDG approach * to compute the diffusion term looping on the planes in the z * direction and adding the flux in z direction at the end. */ ... ... @@ -80,12 +80,12 @@ namespace Nektar { // Strip trailing string "3DHomogeneous1D" to determine 2D diffusion // type, and create a diffusion object for the plane. string name = diffType.substr(0, diffType.length()-15); m_planeDiff = GetDiffusionFactory().CreateInstance(name, name); m_diffType = diffType.substr(0, diffType.length()-15); m_planeDiff = GetDiffusionFactory().CreateInstance(m_diffType, m_diffType); } /** * @brief Initiliase Diffusion3DHomogeneous1D objects and store * @brief Initiliase Diffusion3DHomogeneous1D objects and store * them before starting the time-stepping. * * @param pSession Pointer to session reader. ... ... @@ -114,7 +114,26 @@ namespace Nektar m_homoLen = pFields[0]->GetHomoLen(); m_trans = pFields[0]->GetTransposition(); m_planeCounter = 0; m_planeDiff->SetFluxVectorNS(m_fluxVectorNS); if (m_diffType == "LDG") { // Set viscous flux for LDG m_planeDiff->SetFluxVector(m_fluxVector); } else if (m_diffType == "LDGNS") { // Set viscous flux for LDGNS m_planeDiff->SetFluxVectorNS(m_fluxVectorNS); // Set penalty flux m_planeDiff->SetFluxPenaltyNS(m_fluxPenaltyNS); } else if (m_diffType == "LFRDGNS" || m_diffType == "LFRHUNS" || m_diffType == "LFRSDNS" ) { // Set viscous flux for FR cases m_planeDiff->SetFluxVectorNS(m_fluxVectorNS); } m_fieldsPlane = Array (nConvectiveFields); ... ... @@ -169,7 +188,7 @@ namespace Nektar * using an LDG interface flux and the the flux in the third direction. */ void Diffusion3DHomogeneous1D::v_Diffuse( const int nConvectiveFields, const std::size_t nConvectiveFields, const Array &fields, const Array > &inarray, Array > &outarray, ... ... @@ -181,22 +200,21 @@ namespace Nektar Array tmp(m_numPoints), tmp2; Array > viscHComp; const int nPointsTot = fields[0]->GetNpoints(); int i, j; NekDouble beta; if (m_fluxVectorNS) { viscHComp = Array >(nConvectiveFields); for (i = 0; i < nConvectiveFields - 1; ++i) for (int i = 0; i < nConvectiveFields - 1; ++i) { fields[0]->PhysDeriv(2, inarray[i], m_homoDerivStore[i]); viscHComp[i] = Array(m_numPoints); } } for (i = 0; i < m_numPlanes; ++i) for (int i = 0; i < m_numPlanes; ++i) { // Set up memory references for fields, inarray and outarray for // this plane. ... ... @@ -212,7 +230,7 @@ namespace Nektar m_outarrayPlane[j] = Array( m_numPointsPlane, tmp2 = outarray[j] + m_planePos[i]); } if (m_fluxVectorNS) { ... ... @@ -220,15 +238,43 @@ namespace Nektar } m_planeDiff->Diffuse(nConvectiveFields, m_fieldsPlane, m_inarrayPlane, m_outarrayPlane); if (m_diffType == "LDGNS") { // Store plane Fwd/Bwd traces std::size_t nTracePts = m_fieldsPlane[0]->GetTrace() ->GetTotPoints(); std::size_t nScalar = m_inarrayPlane.num_elements(); Array > Fwd(nScalar); Array > Bwd(nScalar); { for(std::size_t k = 0; k < nScalar; ++k) { Fwd[k] = Array(nTracePts, 0.0); Bwd[k] = Array(nTracePts, 0.0); m_fieldsPlane[k]->GetFwdBwdTracePhys( m_inarrayPlane[k], Fwd[k], Bwd[k]); } } m_planeDiff->Diffuse(nConvectiveFields, m_fieldsPlane, m_inarrayPlane, m_outarrayPlane, Fwd, Bwd); } else { m_planeDiff->Diffuse(nConvectiveFields, m_fieldsPlane, m_inarrayPlane, m_outarrayPlane); } if (m_fluxVectorNS) { Array > > &viscTensor = m_planeDiff->GetFluxTensor(); Array > > &viscTensor = m_planeDiff->GetFluxTensor(); // Extract H (viscTensor[2]) for (int j = 0; j < nConvectiveFields - 1; ++j) ... ... @@ -239,24 +285,25 @@ namespace Nektar } } } if (m_fluxVectorNS) { for (j = 0; j < nConvectiveFields - 1; ++j) for (int j = 0; j < nConvectiveFields - 1; ++j) { fields[j+1]->PhysDeriv(2, viscHComp[j], tmp); Vmath::Vadd(nPointsTot, outarray[j+1], 1, tmp, 1, outarray[j+1], 1); Vmath::Vadd(nPointsTot, outarray[j+1], 1, tmp, 1, outarray[j+1], 1); } } else { for (j = 0; j < nConvectiveFields; ++j) for (int j = 0; j < nConvectiveFields; ++j) { fields[j]->HomogeneousFwdTrans(inarray[j], tmp); for (i = 0; i < m_numPlanes; ++i) for (int i = 0; i < m_numPlanes; ++i) { beta = 2*M_PI*m_trans->GetK(i)/m_homoLen; beta *= beta; ... ...
 ... ... @@ -51,7 +51,7 @@ namespace Nektar new Diffusion3DHomogeneous1D(diffType)); } static std::string type[]; protected: Diffusion3DHomogeneous1D(std::string diffType); ... ... @@ -59,10 +59,10 @@ namespace Nektar std::string m_diffType; SolverUtils::DiffusionSharedPtr m_planeDiff; NekDouble m_homoLen; int m_numPoints; int m_numPointsPlane; int m_numPlanes; int m_planeCounter; std::size_t m_numPoints; std::size_t m_numPointsPlane; std::size_t m_numPlanes; std::size_t m_planeCounter; Array m_planes; Array m_planePos; Array > m_homoDerivStore; ... ... @@ -76,14 +76,14 @@ namespace Nektar LibUtilities::SessionReaderSharedPtr pSession, Array pFields); virtual void v_Diffuse( const int nConvective, const std::size_t nConvective, const Array &fields, const Array > &inarray, Array > &outarray, const Array > &pFwd = NullNekDoubleArrayofArray, const Array > &pBwd = NullNekDoubleArrayofArray); }; }; } } #endif
This diff is collapsed.
 ... ... @@ -51,59 +51,63 @@ namespace Nektar boost::ignore_unused(diffType); return DiffusionSharedPtr(new DiffusionLDG()); } static std::string type; protected: DiffusionLDG(); std::string m_shockCaptureType; std::string m_shockCaptureType; /// Coefficient of penalty term NekDouble m_C11; Array > m_traceNormals; LibUtilities::SessionReaderSharedPtr m_session; virtual void v_InitObject( LibUtilities::SessionReaderSharedPtr pSession, Array pFields); virtual void v_Diffuse( const int nConvective, const std::size_t nConvective, const Array &fields, const Array > &inarray, Array > &outarray, const Array > &pFwd = NullNekDoubleArrayofArray, const Array > &pBwd = NullNekDoubleArrayofArray); virtual void v_NumFluxforScalar( void NumFluxforScalar( const Array &fields, const Array > &ufield, Array > >&uflux, const Array > &pFwd = NullNekDoubleArrayofArray, const Array > &pBwd = NullNekDoubleArrayofArray); virtual void v_WeakPenaltyforScalar( void ApplyScalarBCs( const Array &fields, const int var, const std::size_t var, const Array &ufield, const Array &uplus, const Array &Fwd, const Array &Bwd, Array &penaltyflux); virtual void v_NumFluxforVector( void NumFluxforVector( const Array &fields, const Array > &ufield, Array > >&qfield, Array > &qflux); virtual void v_WeakPenaltyforVector( void ApplyVectorBCs( const Array &fields, const int var, const int dir, const std::size_t var, const std::size_t dir, const Array &qfield, const Array &qtemp, Array &penaltyflux, NekDouble C11); }; const Array &qFwd, const Array &qBwd, Array &penaltyflux); }; } } #endif
This source diff could not be displayed because it is too large. You can view the blob instead.