Commit 480c850d authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'feature/Euler1D' into 'master'

Fixes for Euler 1D solvers



See merge request !565
parents 93da83c0 a5df2009
......@@ -31,6 +31,7 @@ v4.3.3
**CompressibleFlowSolver**:
- Fix issue with residual output (!647)
- Issues with 1D Euler solver fixed (!565)
**Packaging**:
- Fix NekMesh dependencies for DEB package (!650)
......
......@@ -895,11 +895,6 @@ namespace Nektar
Array<OneD, NekDouble> &Fwd,
Array<OneD, NekDouble> &Bwd)
{
// Expansion casts
LocalRegions::Expansion1DSharedPtr exp1D;
LocalRegions::Expansion1DSharedPtr exp1DFirst;
LocalRegions::Expansion1DSharedPtr exp1DLast;
// Counter variables
int n, v;
......@@ -911,10 +906,7 @@ namespace Nektar
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
// Basis shared pointer
LibUtilities::BasisSharedPtr Basis;
// Set forward and backard state to zero
Vmath::Zero(Fwd.num_elements(), Fwd, 1);
Vmath::Zero(Bwd.num_elements(), Bwd, 1);
......@@ -924,12 +916,8 @@ namespace Nektar
// Loop on the elements
for (cnt = n = 0; n < nElements; ++n)
{
exp1D = (*m_exp)[n]->as<LocalRegions::Expansion1D>();
// Set the offset of each element
phys_offset = GetPhys_Offset(n);
Basis = (*m_exp)[n]->GetBasis(0);
for(v = 0; v < 2; ++v, ++cnt)
{
......@@ -1052,24 +1040,19 @@ namespace Nektar
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
// Basis shared pointer
LibUtilities::BasisSharedPtr Basis;
vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
for (n = 0; n < GetExpSize(); ++n)
{
// Basis definition on each element
Basis = (*m_exp)[n]->GetBasis(0);
// Number of coefficients on each element
int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
offset = GetCoeff_Offset(n);
// Implementation for every points except Gauss points
if (Basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
LibUtilities::eGauss_Lagrange)
{
t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
if(negatedFluxNormal[2*n])
......
......@@ -181,9 +181,13 @@ namespace Nektar
switch (normals.num_elements())
{
case 1:
// do nothing
{ // do nothing
const int nq = inarray[0].num_elements();
const int vx = (int)vecLocs[i][0];
Vmath::Vmul (nq, inarray [vx], 1, normals [0], 1,
outarray[vx], 1);
break;
}
case 2:
{
const int nq = inarray[0].num_elements();
......@@ -268,9 +272,13 @@ namespace Nektar
switch (normals.num_elements())
{
case 1:
// do nothing
{ // do nothing
const int nq = normals[0].num_elements();
const int vx = (int)vecLocs[i][0];
Vmath::Vmul (nq, inarray [vx], 1, normals [0], 1,
outarray[vx], 1);
break;
}
case 2:
{
const int nq = normals[0].num_elements();
......
......@@ -33,7 +33,9 @@
//
///////////////////////////////////////////////////////////////////////////////
#include <StdRegions/StdSegExp.h>
#include <StdRegions/StdSegExp.h>
#include <LibUtilities/Foundations/InterpCoeff.h>
using namespace std;
......@@ -238,30 +240,66 @@ namespace Nektar
NekMatrix<NekDouble> B(nquad,m_ncoeffs,m_base[0]->GetBdata(),eWrapper);
out = B * in;
#endif //NEKTAR_USING_DIRECT_BLAS_CALLS
#endif //NEKTAR_USING_DIRECT_BLAS_CALLS
}
}
/** \brief Forward transform from physical quadrature space stored in
* \a inarray and evaluate the expansion coefficients and store in
* \a outarray
*
* Perform a forward transform using a Galerkin projection by taking the
* inner product of the physical points and multiplying by the inverse of
* the mass matrix using the Solve method of the standard matrix
* container holding the local mass matrix, i.e. \f$ {\bf \hat{u}} =
* {\bf M}^{-1} {\bf I} \f$ where \f$ {\bf I}[p] = \int^1_{-1}
* \phi_p(\xi_1) u(\xi_1) d\xi_1 \f$
void StdSegExp::v_ReduceOrderCoeffs(
int numMin,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray)
{
int n_coeffs = inarray.num_elements();
Array<OneD, NekDouble> coeff(n_coeffs);
Array<OneD, NekDouble> coeff_tmp(n_coeffs,0.0);
Array<OneD, NekDouble> tmp;
Array<OneD, NekDouble> tmp2;
int nmodes0 = m_base[0]->GetNumModes();
int numMax = nmodes0;
Vmath::Vcopy(n_coeffs,inarray,1,coeff_tmp,1);
const LibUtilities::PointsKey Pkey0(
nmodes0, LibUtilities::eGaussLobattoLegendre);
LibUtilities::BasisKey b0(m_base[0]->GetBasisType(),nmodes0,Pkey0);
LibUtilities::BasisKey bortho0(LibUtilities::eOrtho_A,nmodes0,Pkey0);
LibUtilities::InterpCoeff1D(
b0, coeff_tmp, bortho0, coeff);
Vmath::Zero(n_coeffs,coeff_tmp,1);
Vmath::Vcopy(numMin,
tmp = coeff,1,
tmp2 = coeff_tmp,1);
LibUtilities::InterpCoeff1D(
bortho0, coeff_tmp, b0, outarray);
}
/**
* \brief Forward transform from physical quadrature space stored in \a
* inarray and evaluate the expansion coefficients and store in \a
* outarray
*
* This function stores the expansion coefficients calculated by the
* transformation in the coefficient space array \a outarray
* Perform a forward transform using a Galerkin projection by taking the
* inner product of the physical points and multiplying by the inverse
* of the mass matrix using the Solve method of the standard matrix
* container holding the local mass matrix, i.e. \f$ {\bf \hat{u}} =
* {\bf M}^{-1} {\bf I} \f$ where \f$ {\bf I}[p] = \int^1_{-1}
* \phi_p(\xi_1) u(\xi_1) d\xi_1 \f$
*
* \param inarray: array of physical quadrature points to be transformed
* This function stores the expansion coefficients calculated by the
* transformation in the coefficient space array \a outarray
*
* \param outarray: the coeffficients of the expansion
*/
* \param inarray: array of physical quadrature points to be transformed
* \param outarray: the coeffficients of the expansion
*/
void StdSegExp::v_FwdTrans(const Array<OneD, const NekDouble>& inarray,
Array<OneD, NekDouble> &outarray)
{
......
......@@ -215,6 +215,15 @@ namespace Nektar
STD_REGIONS_EXPORT virtual DNekMatSharedPtr v_CreateStdMatrix(
const StdMatrixKey &mkey);
//---------------------------------------
// Operator evaluation functions
//---------------------------------------
STD_REGIONS_EXPORT virtual void v_ReduceOrderCoeffs(
int numMin,
const Array<OneD, const NekDouble> &inarray,
Array<OneD, NekDouble> &outarray);
private:
};
......
......@@ -34,6 +34,7 @@ IF( NEKTAR_SOLVER_COMPRESSIBLE_FLOW )
ADD_NEKTAR_TEST (CylinderSubsonicMix)
ADD_NEKTAR_TEST (CylinderSubsonic_P3)
ADD_NEKTAR_TEST_LENGTHY(CylinderSubsonic_P8)
ADD_NEKTAR_TEST (Euler1D)
ADD_NEKTAR_TEST (IsentropicVortex16_P3)
ADD_NEKTAR_TEST (IsentropicVortex_FRDG_SEM)
ADD_NEKTAR_TEST_LENGTHY(IsentropicVortex_FRSD_SEM)
......
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>Euler 1D P=3, WeakDG, MODIFIED</description>
<executable>CompressibleFlowSolver</executable>
<parameters>Euler1D.xml</parameters>
<files>
<file description="Session File">Euler1D.xml</file>
</files>
<metrics>
<metric type="L2" id="1">
<value variable="rho" tolerance="1e-12">1.98838e-06</value>
<value variable="rhou" tolerance="1e-12">0.00067684</value>
<value variable="E" tolerance="1e-12">0.575708</value>
</metric>
<metric type="Linf" id="2">
<value variable="rho" tolerance="1e-12">1.98524e-05</value>
<value variable="rhou" tolerance="1e-12">0.00675771</value>
<value variable="E" tolerance="1e-12">5.74799</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<NEKTAR>
<GEOMETRY DIM="1" SPACE="1">
<VERTEX>
<V ID="0"> 0.000e+00 0.000e+00 0.000e+00</V>
<V ID="1"> 0.020e+00 0.000e+00 0.000e+00</V>
<V ID="2"> 0.040e+00 0.000e+00 0.000e+00</V>
<V ID="3"> 0.060e+00 0.000e+00 0.000e+00</V>
<V ID="4"> 0.080e+00 0.000e+00 0.000e+00</V>
<V ID="5"> 0.100e+00 0.000e+00 0.000e+00</V>
<V ID="6"> 0.120e+00 0.000e+00 0.000e+00</V>
<V ID="7"> 0.140e+00 0.000e+00 0.000e+00</V>
<V ID="8"> 0.160e+00 0.000e+00 0.000e+00</V>
<V ID="9"> 0.180e+00 0.000e+00 0.000e+00</V>
<V ID="10"> 0.200e+00 0.000e+00 0.000e+00</V>
<V ID="11"> 0.220e+00 0.000e+00 0.000e+00</V>
<V ID="12"> 0.240e+00 0.000e+00 0.000e+00</V>
<V ID="13"> 0.260e+00 0.000e+00 0.000e+00</V>
<V ID="14"> 0.280e+00 0.000e+00 0.000e+00</V>
<V ID="15"> 0.300e+00 0.000e+00 0.000e+00</V>
<V ID="16"> 0.320e+00 0.000e+00 0.000e+00</V>
<V ID="17"> 0.340e+00 0.000e+00 0.000e+00</V>
<V ID="18"> 0.360e+00 0.000e+00 0.000e+00</V>
<V ID="19"> 0.380e+00 0.000e+00 0.000e+00</V>
<V ID="20"> 0.400e+00 0.000e+00 0.000e+00</V>
</VERTEX>
<ELEMENT>
<S ID="0"> 0 1 </S>
<S ID="1"> 1 2 </S>
<S ID="2"> 2 3 </S>
<S ID="3"> 3 4 </S>
<S ID="4"> 4 5 </S>
<S ID="5"> 5 6 </S>
<S ID="6"> 6 7 </S>
<S ID="7"> 7 8 </S>
<S ID="8"> 8 9 </S>
<S ID="9"> 9 10 </S>
<S ID="10"> 10 11 </S>
<S ID="11"> 11 12 </S>
<S ID="12"> 12 13 </S>
<S ID="13"> 13 14 </S>
<S ID="14"> 14 15 </S>
<S ID="15"> 15 16 </S>
<S ID="16"> 16 17 </S>
<S ID="17"> 17 18 </S>
<S ID="18"> 18 19 </S>
<S ID="19"> 19 20 </S>
</ELEMENT>
<COMPOSITE>
<C ID="0"> S[0-19] </C>
<C ID="1"> V[0] </C>
<C ID="2"> V[20] </C>
</COMPOSITE>
<DOMAIN>
<D ID="0"> C[0] </D>
</DOMAIN>
</GEOMETRY>
<EXPANSIONS>
<E COMPOSITE="C[0]" NUMMODES="3" FIELDS="rho,rhou,E" TYPE="MODIFIED" />
</EXPANSIONS>
<CONDITIONS>
<PARAMETERS>
<P> TimeStep = 1e-6 </P>
<P> FinTime = 40e-6 </P>
<P> NumSteps = FinTime/TimeStep </P>
<P> IO_CheckSteps = 100 </P>
<P> IO_InfoSteps = 10 </P>
<P> Gamma = 1.4 </P>
<P> pInf = 101325 </P>
<P> rhoInf = 1.225 </P>
<P> u0 = 0.1 </P>
<P> uInf = u0 </P>
<P> uInfL = u0 + 0.01 </P>
<P> uInfR = u0 </P>
</PARAMETERS>
<SOLVERINFO>
<I PROPERTY="EQTYPE" VALUE="EulerCFE" />
<I PROPERTY="Projection" VALUE="DisContinuous" />
<I PROPERTY="AdvectionType" VALUE="WeakDG" />
<I PROPERTY="DiffusionType" VALUE="LDGNS" />
<I PROPERTY="TimeIntegrationMethod" VALUE="ClassicalRungeKutta4" />
<I PROPERTY="UpwindType" VALUE="ExactToro"/>
<I PROPERTY="ProblemType" VALUE="General"/>
<I PROPERTY="ViscosityType" VALUE="Constant"/>
</SOLVERINFO>
<VARIABLES>
<V ID="0"> rho </V>
<V ID="1"> rhou </V>
<V ID="2"> E </V>
</VARIABLES>
<BOUNDARYREGIONS>
<B ID="0"> C[1] </B>
<B ID="1"> C[2] </B>
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="rho" VALUE="rhoInf" />
<D VAR="rhou" VALUE="rhoInf*uInfL" />
<D VAR="E" VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInfL*uInfL)" />
</REGION>
<REGION REF="1">
<D VAR="rho" VALUE="rhoInf" />
<D VAR="rhou" VALUE="rhoInf*uInfR" />
<D VAR="E" VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInfR*uInfR)" />
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="InitialConditions">
<E VAR="rho" DOMAIN="0" VALUE="rhoInf" />
<E VAR="rhou" DOMAIN="0" VALUE="rhoInf*u0" />
<E VAR="E" DOMAIN="0" VALUE="pInf/(Gamma-1)+0.5*rhoInf*(u0*u0)" />
</FUNCTION>
<FUNCTION NAME="ExactSolution">
<E VAR="rho" DOMAIN="0" VALUE="rhoInf" />
<E VAR="rhou" DOMAIN="0" VALUE="rhoInf*uInf" />
<E VAR="E" DOMAIN="0" VALUE="pInf/(Gamma-1)+0.5*rhoInf*(uInf*uInf)" />
</FUNCTION>
</CONDITIONS>
</NEKTAR>
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment