Commit 27d8ae65 authored by Gianmarco Mengaldo's avatar Gianmarco Mengaldo
Browse files

Regression tests added for both ADRSolver (diffusion operator) and...

Regression tests added for both ADRSolver (diffusion operator) and CompressibleFlowSolver (NS). Solved a bug in LDG (Vn instead of m_traceNormals into the Upwind method).
parent 13aa1f16
......@@ -61,6 +61,7 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
//cout<<setprecision(16);
int i, j, k;
int nDim = fields[0]->GetCoordim(0);
int nPts = fields[0]->GetTotPoints();
......@@ -168,6 +169,13 @@ namespace Nektar
}
fields[0]->GetTrace()->GetNormals(m_traceNormals);
// Get the normal velocity Vn
for(i = 0; i < nDim; ++i)
{
Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
Vn, 1, Vn, 1);
}
// Get the sign of (v \cdot n), v = an arbitrary vector
// Evaluate upwind flux:
// uflux = \hat{u} \phi \cdot u = u^{(+,-)} n
......@@ -186,7 +194,7 @@ namespace Nektar
// edge::eForward, if V*n<0 <=> V*n_F<0, pick uflux = uBwd
// edge::eBackward, if V*n<0 <=> V*n_B>=0, pick uflux = uBwd
fields[i]->GetTrace()->Upwind(m_traceNormals[j],
fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
Fwd, Bwd, fluxtemp);
// Imposing weak boundary condition with flux
......@@ -315,6 +323,13 @@ namespace Nektar
}
fields[0]->GetTrace()->GetNormals(m_traceNormals);
// Get the normal velocity Vn
for(i = 0; i < nDim; ++i)
{
Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
Vn, 1, Vn, 1);
}
// Evaulate upwind flux:
// qflux = \hat{q} \cdot u = q \cdot n - C_(11)*(u^+ - u^-)
for (i = 0; i < nvariables; ++i)
......@@ -337,7 +352,7 @@ namespace Nektar
// edge::eBackward, if V*n<0 <=> V*n_B>=0, pick
// qflux = qFwd = q+
fields[i]->GetTrace()->Upwind(m_traceNormals[j],
fields[i]->GetTrace()->Upwind(/*m_traceNormals[j]*/Vn,
qBwd, qFwd,
qfluxtemp);
......
......@@ -84,7 +84,7 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
cout<<setprecision(16);
//cout<<setprecision(16);
int i, j;
int nDim = fields[0]->GetCoordim(0);
int nScalars = inarray.num_elements();
......
......@@ -35,6 +35,7 @@
#include <SolverUtils/Diffusion/DiffusionLFR.h>
#include <LibUtilities/Polylib/Polylib.h>
#include <boost/math/special_functions/gamma.hpp>
#include <iostream>
#include <iomanip>
......@@ -768,9 +769,9 @@ namespace Nektar
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
//cout<<setprecision(16);
int i, j, n, z;
int nLocalSolutionPts, phys_offset;
Array<TwoD, const NekDouble> gmat;
Array<OneD, const NekDouble> jac;
Array<OneD, NekDouble> auxArray1, auxArray2, auxArray3;
......@@ -802,10 +803,10 @@ namespace Nektar
iuFluxO1[i][j] = Array<OneD, NekDouble>(nTracePts, 0.0);
}
}
// Compute interface numerical fluxes for inarray in physical space
v_NumFluxforScalar(fields, inarray, iuFluxO1);
switch(nDim)
{
// 1D problems
......
......@@ -35,6 +35,7 @@
#include <SolverUtils/Diffusion/DiffusionLFRNS.h>
#include <LibUtilities/Polylib/Polylib.h>
#include <boost/math/special_functions/gamma.hpp>
#include <iostream>
#include <iomanip>
......@@ -141,13 +142,13 @@ namespace Nektar
// Extract the Q factors at each edge point
pFields[0]->GetExp(n)->GetEdgeQFactors(
0, auxArray1 = m_Q2D_e0[n]);
0, auxArray1 = m_Q2D_e0[n]);
pFields[0]->GetExp(n)->GetEdgeQFactors(
1, auxArray1 = m_Q2D_e1[n]);
1, auxArray1 = m_Q2D_e1[n]);
pFields[0]->GetExp(n)->GetEdgeQFactors(
2, auxArray1 = m_Q2D_e2[n]);
2, auxArray1 = m_Q2D_e2[n]);
pFields[0]->GetExp(n)->GetEdgeQFactors(
3, auxArray1 = m_Q2D_e3[n]);
3, auxArray1 = m_Q2D_e3[n]);
}
break;
}
......@@ -775,7 +776,7 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > &outarray)
{
cout<<setprecision(16);
int i, j, n;
int i, j, n, z;
int nLocalSolutionPts, phys_offset;
Array<TwoD, const NekDouble> gmat;
......@@ -1001,6 +1002,10 @@ namespace Nektar
divFC(nConvectiveFields);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
tmp(nConvectiveFields);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
tmp1(nConvectiveFields);
Array<OneD, Array<OneD, Array<OneD, NekDouble> > >
tmp2(nConvectiveFields);
derivativesO1[0] = Array<OneD, Array<OneD, NekDouble> >(
nScalars);
......@@ -1021,8 +1026,9 @@ namespace Nektar
nDim);
BderivativesO1[i] = Array<OneD, Array<OneD, NekDouble> >(
nDim);
tmp[i] = Array<OneD, Array<OneD, NekDouble> >(
nDim);
tmp[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
tmp1[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
tmp2[i] = Array<OneD, Array<OneD, NekDouble> >(nDim);
for (j = 0; j < nDim; ++j)
{
......@@ -1038,7 +1044,147 @@ namespace Nektar
Array<OneD, NekDouble>(nPts, 0.0);
tmp[i][j] =
Array<OneD, NekDouble>(nPts, 0.0);
tmp1[i][j] =
Array<OneD, NekDouble>(nPts, 0.0);
tmp2[i][j] =
Array<OneD, NekDouble>(nPts, 0.0);
// Computing the physical first-order discountinuous
// derivatives
for (n = 0; n < nElements; n++)
{
// Discontinuous flux
nLocalSolutionPts = fields[0]->GetExp(n)->
GetTotPoints();
phys_offset = fields[0]->GetPhys_Offset(n);
jac = fields[0]->GetExp(n)->GetGeom2D()->
GetJac();
gmat = fields[0]->GetExp(n)->GetGeom2D()->
GetGmat();
// Temporary vectors
Array<OneD, NekDouble> u1_hat(
nLocalSolutionPts, 0.0);
Array<OneD, NekDouble> u2_hat(
nLocalSolutionPts, 0.0);
switch (j)
{
case 0:
{
if (fields[0]->GetExp(n)->GetGeom2D()->
GetGtype() ==
SpatialDomains::eDeformed)
{
for (z = 0; z < nLocalSolutionPts;
z++)
{
u1_hat[z] =
(inarray[i][z+phys_offset]
* gmat[0][z]) * jac[z];
u2_hat[z] =
(inarray[i][z+phys_offset]
* gmat[1][z]) * jac[z];
}
}
else
{
for (z = 0; z < nLocalSolutionPts;
z++)
{
u1_hat[z] =
(inarray[i][z+phys_offset]
* gmat[0][0]) * jac[0];
u2_hat[z] =
(inarray[i][z+phys_offset]
* gmat[1][0])*jac[0];
}
}
break;
}
case 1:
{
if (fields[0]->GetExp(n)->GetGeom2D()->
GetGtype() ==
SpatialDomains::eDeformed)
{
for (z = 0; z < nLocalSolutionPts;
z++)
{
u1_hat[z] =
(inarray[i][z+phys_offset]
* gmat[2][z]) * jac[z];
u2_hat[z] =
(inarray[i][z+phys_offset]
* gmat[3][z]) * jac[z];
}
}
else
{
for (z = 0; z < nLocalSolutionPts;
z++)
{
u1_hat[z] =
(inarray[i][z+phys_offset]
* gmat[2][0]) * jac[0];
u2_hat[z] =
(inarray[i][z+phys_offset]
* gmat[3][0])*jac[0];
}
}
break;
}
}
fields[i]->GetExp(n)->StdPhysDeriv(0,
auxArray1 = u1_hat,
auxArray2 = tmp1[i][j] + phys_offset);
fields[i]->GetExp(n)->StdPhysDeriv(1,
auxArray1 = u2_hat,
auxArray2 = tmp2[i][j] + phys_offset);
}
Vmath::Vadd(nPts,
auxArray1 = tmp1[i][j], 1,
auxArray2 = tmp2[i][j], 1,
DinarrayO1[i][j], 1);
// Multiply by the metric terms
for (n = 0; n < nElements; ++n)
{
nLocalSolutionPts = fields[0]->
GetExp(n)->GetTotPoints();
phys_offset = fields[0]->GetPhys_Offset(n);
jac = fields[0]->GetExp(n)->GetGeom2D()->GetJac();
gmat = fields[0]->GetExp(n)->GetGeom2D()->GetGmat();
if (fields[0]->GetExp(n)->GetGeom2D()->
GetGtype() == SpatialDomains::eDeformed)
{
for (z = 0; z < nLocalSolutionPts; z++)
{
DinarrayO1[i][j][phys_offset + z] =
DinarrayO1[i][j][phys_offset + z]/jac[z];
}
}
else
{
Vmath::Smul(nLocalSolutionPts, 1/jac[0],
auxArray1 = DinarrayO1[i][j]
+ phys_offset, 1,
auxArray2 = DinarrayO1[i][j]
+ phys_offset, 1);
}
}
/*
// Computing the physical first-order discountinuous
// derivatives
for (n = 0; n < nElements; n++)
......@@ -1049,6 +1195,7 @@ namespace Nektar
auxArray1 = inarray[i] + phys_offset,
auxArray2 = DinarrayO1[i][j] + phys_offset);
}
*/
// Computing the standard first-order correction
// derivatives
......
......@@ -74,8 +74,21 @@ IF( NEKTAR_SOLVER_ADR )
ADD_NEKTAR_TEST(Advection_m12_Order2)
ADD_NEKTAR_TEST(Advection_m14_DG_Order4)
ADD_NEKTAR_TEST(Advection_m14_Order4)
ADD_NEKTAR_TEST(ExDiffusion_m3)
ADD_NEKTAR_TEST(ExDiffusion_m8)
ADD_NEKTAR_TEST(ExDiffusion_1D_LDG)
ADD_NEKTAR_TEST(ExDiffusion_1D_LFRDG)
ADD_NEKTAR_TEST(ExDiffusion_1D_LFRHU)
ADD_NEKTAR_TEST(ExDiffusion_1D_LFRSD)
ADD_NEKTAR_TEST(ExDiffusion_2D_LDG_hybrid_m3)
ADD_NEKTAR_TEST(ExDiffusion_2D_LDG_hybrid_m8)
ADD_NEKTAR_TEST(ExDiffusion_2D_LDG_regular_Neumann)
ADD_NEKTAR_TEST(ExDiffusion_2D_LFRDG_regular_Neumann)
ADD_NEKTAR_TEST(ExDiffusion_2D_LFRHU_regular_Neumann)
ADD_NEKTAR_TEST(ExDiffusion_2D_LFRSD_regular_Neumann)
ADD_NEKTAR_TEST(ExDiffusion_2D_LDG_deformed)
ADD_NEKTAR_TEST(ExDiffusion_2D_LFRDG_deformed)
ADD_NEKTAR_TEST(ExDiffusion_2D_LFRHU_deformed)
ADD_NEKTAR_TEST(ExDiffusion_2D_LFRSD_deformed)
ADD_NEKTAR_TEST(Helmholtz1D_8modes_DG)
ADD_NEKTAR_TEST(Helmholtz1D_8modes)
ADD_NEKTAR_TEST(Helmholtz1D_8nodes)
......
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>2D unsteady DG explicit diffusion, order 4, P=8</description>
<description>1D unsteady LDG explicit diffusion, order 3, P=2</description>
<executable>ADRSolver</executable>
<parameters>ExDiffusion_m8.xml</parameters>
<parameters>ExDiffusion_1D_LDG.xml</parameters>
<files>
<file description="Session File">ExDiffusion_m8.xml</file>
<file description="Session File">ExDiffusion_1D_LDG.xml</file>
</files>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">0.000337147</value>
<value variable="u" tolerance="1e-12">0.00196469</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">0.00279079</value>
<value variable="u" tolerance="1e-12">0.00925869</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<NEKTAR>
<GEOMETRY DIM="1" SPACE="1">
<VERTEX>
<V ID="0"> -1.0 0.0 0.0</V>
<V ID="1"> -0.8 0.0 0.0</V>
<V ID="2"> -0.6 0.0 0.0</V>
<V ID="3"> -0.4 0.0 0.0</V>
<V ID="4"> -0.2 0.0 0.0</V>
<V ID="5"> 0.0 0.0 0.0</V>
<V ID="6"> 0.2 0.0 0.0</V>
<V ID="7"> 0.4 0.0 0.0</V>
<V ID="8"> 0.6 0.0 0.0</V>
<V ID="9"> 0.8 0.0 0.0</V>
<V ID="10"> 1.0 0.0 0.0</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>
</ELEMENT>
<COMPOSITE>
<C ID="0"> S[0-9] </C>
<C ID="1"> V[0] </C>
<C ID="2"> V[10] </C>
</COMPOSITE>
<DOMAIN> C[0] </DOMAIN>
</GEOMETRY>
<EXPANSIONS>
<E COMPOSITE="C[0]" FIELDS="u" TYPE="GLL_LAGRANGE_SEM" NUMMODES="3"/>
</EXPANSIONS>
<CONDITIONS>
<PARAMETERS>
<P> TimeStep = 0.000001 </P>
<P> NumSteps = 100 </P>
<P> FinTime = TimeStep * NumSteps </P>
<P> IO_CheckSteps = 100000 </P>
<P> IO_InfoSteps = 100000 </P>
</PARAMETERS>
<SOLVERINFO>
<I PROPERTY="EQTYPE" VALUE="UnsteadyDiffusion" />
<I PROPERTY="Projection" VALUE="DisContinuous" />
<I PROPERTY="DiffusionType" VALUE="LDG" />
<I PROPERTY="DiffusionAdvancement" VALUE="Explicit" />
<I PROPERTY="TimeIntegrationMethod" VALUE="ClassicalRungeKutta4"/>
</SOLVERINFO>
<VARIABLES>
<V ID="0"> u </V>
</VARIABLES>
<BOUNDARYREGIONS>
<B ID="0"> C[1] </B>
<B ID="1"> C[2] </B>
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" USERDEFINEDTYPE="TimeDependent"
VALUE="exp(-2.0*PI*PI*t)*sin(PI*x)" />
</REGION>
<REGION REF="1">
<D VAR="u" USERDEFINEDTYPE="TimeDependent"
VALUE="exp(-2.0*PI*PI*t)*sin(PI*x)" />
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="InitialConditions">
<E VAR="u" VALUE="sin(PI*x)" />
</FUNCTION>
<FUNCTION NAME="ExactSolution">
<E VAR="u" VALUE="exp(-2.0*PI*PI*FinTime)*sin(PI*x)" />
</FUNCTION>
</CONDITIONS>
</NEKTAR>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>1D unsteady LFRDG explicit diffusion, order 3, P=2</description>
<executable>ADRSolver</executable>
<parameters>ExDiffusion_1D_LFRDG.xml</parameters>
<files>
<file description="Session File">ExDiffusion_1D_LFRDG.xml</file>
</files>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">0.00295</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">0.0137468</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<NEKTAR>
<GEOMETRY DIM="1" SPACE="1">
<VERTEX>
<V ID="0"> -1.0 0.0 0.0</V>
<V ID="1"> -0.8 0.0 0.0</V>
<V ID="2"> -0.6 0.0 0.0</V>
<V ID="3"> -0.4 0.0 0.0</V>
<V ID="4"> -0.2 0.0 0.0</V>
<V ID="5"> 0.0 0.0 0.0</V>
<V ID="6"> 0.2 0.0 0.0</V>
<V ID="7"> 0.4 0.0 0.0</V>
<V ID="8"> 0.6 0.0 0.0</V>
<V ID="9"> 0.8 0.0 0.0</V>
<V ID="10"> 1.0 0.0 0.0</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>
</ELEMENT>
<COMPOSITE>
<C ID="0"> S[0-9] </C>
<C ID="1"> V[0] </C>
<C ID="2"> V[10] </C>
</COMPOSITE>
<DOMAIN> C[0] </DOMAIN>
</GEOMETRY>
<EXPANSIONS>
<E COMPOSITE="C[0]" FIELDS="u" TYPE="GLL_LAGRANGE_SEM" NUMMODES="3"/>
</EXPANSIONS>
<CONDITIONS>
<PARAMETERS>
<P> TimeStep = 0.000001 </P>
<P> NumSteps = 100 </P>
<P> FinTime = TimeStep * NumSteps </P>
<P> IO_CheckSteps = 100000 </P>
<P> IO_InfoSteps = 100000 </P>
</PARAMETERS>
<SOLVERINFO>
<I PROPERTY="EQTYPE" VALUE="UnsteadyDiffusion" />
<I PROPERTY="Projection" VALUE="DisContinuous" />
<I PROPERTY="DiffusionType" VALUE="LFRDG" />
<I PROPERTY="DiffusionAdvancement" VALUE="Explicit" />
<I PROPERTY="TimeIntegrationMethod" VALUE="ClassicalRungeKutta4"/>
</SOLVERINFO>
<VARIABLES>
<V ID="0"> u </V>
</VARIABLES>
<BOUNDARYREGIONS>
<B ID="0"> C[1] </B>
<B ID="1"> C[2] </B>
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<D VAR="u" USERDEFINEDTYPE="TimeDependent"
VALUE="exp(-2.0*PI*PI*t)*sin(PI*x)" />
</REGION>
<REGION REF="1">
<D VAR="u" USERDEFINEDTYPE="TimeDependent"
VALUE="exp(-2.0*PI*PI*t)*sin(PI*x)" />
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="InitialConditions">
<E VAR="u" VALUE="sin(PI*x)" />
</FUNCTION>
<FUNCTION NAME="ExactSolution">
<E VAR="u" VALUE="exp(-2.0*PI*PI*FinTime)*sin(PI*x)" />
</FUNCTION>
</CONDITIONS>
</NEKTAR>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>1D unsteady LFRHU explicit diffusion, order 3, P=2</description>
<executable>ADRSolver</executable>
<parameters>ExDiffusion_1D_LFRHU.xml</parameters>
<files>
<file description="Session File">ExDiffusion_1D_LFRHU.xml</file>
</files>
<metrics>
<metric type="L2" id="1">
<value variable="u" tolerance="1e-12">0.00196469</value>
</metric>
<metric type="Linf" id="2">
<value variable="u" tolerance="1e-12">0.00925869</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8" ?>
<NEKTAR>