Commit daf91d19 authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Updates as part of substepping implementation and redefining...

Updates as part of substepping implementation and redefining eDiscontinuousGalerking to be eDiscontinuous


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4011 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent 1ae11094
......@@ -430,7 +430,7 @@ namespace Nektar
int npoints = m_fields[0]->GetNpoints();
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
switch(m_equationType)
{
......@@ -1012,7 +1012,8 @@ namespace Nektar
fields[i] = m_fields[i]->UpdateCoeffs();
}
if(m_projectionType==MultiRegions::eGalerkin)
if((m_projectionType==MultiRegions::eGalerkin)||
(m_projectionType==MultiRegions::eMixed_CG_Discontinuous))
{
// calculate the variable u* = Mu
// we are going to TimeIntegrate this new variable u*
......@@ -1105,7 +1106,8 @@ namespace Nektar
m_time += m_timestep;
if(m_projectionType==MultiRegions::eGalerkin)
if((m_projectionType==MultiRegions::eGalerkin)||
(m_projectionType==MultiRegions::eMixed_CG_Discontinuous))
{
// Project the solution u* onto the boundary conditions to
// obtain the actual solution
......
......@@ -87,7 +87,7 @@ namespace Nektar
switch (m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
int ncoeffs = inarray[0].num_elements();
Array<OneD, Array<OneD, NekDouble> > WeakAdv(nvariables);
......@@ -111,6 +111,7 @@ namespace Nektar
break;
}
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
// Calculate -V\cdot Grad(u);
for(j = 0; j < nvariables; ++j)
......@@ -136,7 +137,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
// Just copy over array
int npoints = GetNpoints();
......@@ -148,6 +149,7 @@ namespace Nektar
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
......
......@@ -89,17 +89,18 @@ namespace Nektar
switch (m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
{
dofs = ncoeffs;
break;
}
case MultiRegions::eGalerkin:
{
dofs = GetContNcoeffs();
UseContCoeffs = true;
break;
}
case MultiRegions::eDiscontinuous:
{
dofs = ncoeffs;
break;
}
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
dofs = GetContNcoeffs();
UseContCoeffs = true;
break;
}
}
cout << endl;
......@@ -127,50 +128,51 @@ namespace Nektar
switch (m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
WeakDGAdvection(inarray, WeakAdv,true,true,1);
m_fields[0]->MultiplyByElmtInvMass(WeakAdv[0],WeakAdv[0]);
m_fields[0]->BwdTrans(WeakAdv[0],outarray[0]);
Vmath::Neg(npoints,outarray[0],1);
break;
WeakDGAdvection(inarray, WeakAdv,true,true,1);
m_fields[0]->MultiplyByElmtInvMass(WeakAdv[0],WeakAdv[0]);
m_fields[0]->BwdTrans(WeakAdv[0],outarray[0]);
Vmath::Neg(npoints,outarray[0],1);
break;
}
case MultiRegions::eGalerkin:
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
// Calculate -V\cdot Grad(u);
// Calculate -V\cdot Grad(u);
for(i = 0; i < nvariables; ++i)
{
//Projection
m_fields[i]->FwdTrans(inarray[i],WeakAdv[i]);
m_fields[i]->BwdTrans_IterPerExp(WeakAdv[i],tmp[i]);
//Advection operator
AdvectionNonConservativeForm(m_velocity,tmp[i],outarray[i]);
Vmath::Neg(npoints,outarray[i],1);
//m_fields[i]->MultiplyByInvMassMatrix(WeakAdv[i],WeakAdv[i]);
//Projection
m_fields[i]->FwdTrans(outarray[i],WeakAdv[i]);
m_fields[i]->BwdTrans_IterPerExp(WeakAdv[i],outarray[i]);
m_fields[i]->FwdTrans(inarray[i],WeakAdv[i]);
m_fields[i]->BwdTrans_IterPerExp(WeakAdv[i],tmp[i]);
//Advection operator
AdvectionNonConservativeForm(m_velocity,tmp[i],outarray[i]);
Vmath::Neg(npoints,outarray[i],1);
//m_fields[i]->MultiplyByInvMassMatrix(WeakAdv[i],WeakAdv[i]);
//Projection
m_fields[i]->FwdTrans(outarray[i],WeakAdv[i]);
m_fields[i]->BwdTrans_IterPerExp(WeakAdv[i],outarray[i]);
}
break;
break;
}
}
/// The result is stored in outarray (is the j-th columns of the weak advection operator).
/// We now store it in MATRIX(j)
Vmath::Vcopy(npoints,&(outarray[0][0]),1,&(MATRIX[j]),npoints);
/// Set the j-th entry of inarray back to zero
inarray[0][j] = 0.0;
/// The result is stored in outarray (is the j-th columns of the weak advection operator).
/// We now store it in MATRIX(j)
Vmath::Vcopy(npoints,&(outarray[0][0]),1,&(MATRIX[j]),npoints);
/// Set the j-th entry of inarray back to zero
inarray[0][j] = 0.0;
}
////////////////////////////////////////////////////////////////////////////////
/// Calulating the eigenvalues of the weak advection operator stored in (MATRIX)
/// using Lapack routines
......
......@@ -86,7 +86,7 @@ namespace Nektar
break;
}
// Discontinuous field
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
string advName;
string riemName;
......@@ -214,7 +214,7 @@ namespace Nektar
switch(m_projectionType)
{
// Discontinuous projection
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
// Number of quadrature points
int nQuadraturePts = GetNpoints();
......@@ -226,9 +226,10 @@ namespace Nektar
}
break;
}
// Continuous projection
case MultiRegions::eGalerkin:
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
......
......@@ -86,7 +86,7 @@ namespace Nektar
switch (m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
int ncoeffs = GetNcoeffs();
Array<OneD, Array<OneD, NekDouble> > WeakAdv(nvariables);
......@@ -109,7 +109,8 @@ namespace Nektar
break;
}
case MultiRegions::eGalerkin:
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
// Calculate -V\cdot Grad(u);
for(i = 0; i < nvariables; ++i)
......
......@@ -98,7 +98,7 @@ namespace Nektar
switch (m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
int ncoeffs = inarray[0].num_elements();
Array<OneD, Array<OneD, NekDouble> > WeakAdv(nvariables);
......@@ -121,6 +121,7 @@ namespace Nektar
break;
}
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
ASSERTL0(false, "Explicit Galerkin diffusion not set up.");
}
......@@ -174,7 +175,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
// Just copy over array
int npoints = GetNpoints();
......@@ -186,6 +187,7 @@ namespace Nektar
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
......
......@@ -74,7 +74,7 @@ namespace Nektar
break;
}
// Discontinuous field
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
string advName;
string riemName;
......@@ -224,7 +224,7 @@ namespace Nektar
switch(m_projectionType)
{
// Discontinuous projection
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
// Number of quadrature points
int nQuadraturePts = GetNpoints();
......@@ -239,6 +239,7 @@ namespace Nektar
// Continuous projection
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
Array<OneD, NekDouble> coeffs(m_fields[0]->GetNcoeffs());
......
......@@ -87,7 +87,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
//-------------------------------------------------------
//inarray in physical space
......@@ -138,7 +138,8 @@ namespace Nektar
}
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
Array<OneD, Array<OneD, NekDouble> > physarray(nvariables);
Array<OneD, Array<OneD, NekDouble> > modarray(nvariables);
......@@ -196,14 +197,15 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
ConservativeToPrimitive(inarray,outarray);
SetBoundaryConditions(outarray,time);
PrimitiveToConservative(outarray,outarray);
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
ConservativeToPrimitive(inarray,outarray);
EquationSystem::SetBoundaryConditions(time);
......@@ -693,27 +695,27 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
{
GetSource(source,m_time);
//Source term solely for the p' equation (outarray[0])
m_fields[0]->IProductWRTBase(source,source);
Vmath::Vadd(ncoeffs,source,1,outarray[0],1,outarray[0],1);
}
case MultiRegions::eDiscontinuous:
{
GetSource(source,m_time);
//Source term solely for the p' equation (outarray[0])
m_fields[0]->IProductWRTBase(source,source);
Vmath::Vadd(ncoeffs,source,1,outarray[0],1,outarray[0],1);
}
break;
case MultiRegions::eGalerkin:
{
GetSource(source,m_time);
//Source term solely for the p' equation (outarray[0])
Vmath::Vadd(ncoeffs,source,1,outarray[0],1,outarray[0],1); }
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
GetSource(source,m_time);
//Source term solely for the p' equation (outarray[0])
Vmath::Vadd(ncoeffs,source,1,outarray[0],1,outarray[0],1); }
break;
default:
ASSERTL0(false,"Unknown projection scheme for the APE");
default:
ASSERTL0(false,"Unknown projection scheme for the APE");
break;
}
}
......
......@@ -93,8 +93,8 @@ namespace Nektar
"Invalid time integration type.");
// if discontinuous Galerkin determine numerical flux to use
if (m_projectionType == MultiRegions::eDiscontinuousGalerkin)
// if discontinuous determine numerical flux to use
if (m_projectionType == MultiRegions::eDiscontinuous)
{
ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
"No UPWINDTYPE defined in session.");
......
......@@ -120,7 +120,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
//-------------------------------------------------------
//inarray in physical space
......@@ -154,6 +154,7 @@ namespace Nektar
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
ASSERTL0(false,"Continouos scheme not implemented for EulerCFE");
break;
default:
......@@ -172,7 +173,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
// Just copy over array
int npoints = GetNpoints();
......@@ -185,6 +186,7 @@ namespace Nektar
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
ASSERTL0(false,"No Continuous Galerkin for Euler equations");
break;
......
......@@ -159,7 +159,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
// Just copy over array
int npoints = GetNpoints();
......@@ -172,6 +172,7 @@ namespace Nektar
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
ASSERTL0(false,"No Continuous Galerkin for Euler equations");
break;
......
......@@ -120,7 +120,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
//-------------------------------------------------------
//inarray in physical space
......@@ -154,6 +154,7 @@ namespace Nektar
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
ASSERTL0(false,"Continouos scheme not implemented for NavierStokesCFE");
break;
default:
......@@ -172,7 +173,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
// Just copy over array
int npoints = GetNpoints();
......@@ -185,6 +186,7 @@ namespace Nektar
}
break;
case MultiRegions::eGalerkin:
case MultiRegions::eMixed_CG_Discontinuous:
{
ASSERTL0(false,"No Continuous Galerkin for NavierStokes equations");
break;
......
......@@ -94,7 +94,7 @@ namespace Nektar
// if discontinuous Galerkin determine numerical flux to use
if (m_projectionType == MultiRegions::eDiscontinuousGalerkin)
if (m_projectionType == MultiRegions::eDiscontinuous)
{
ASSERTL0(m_session->DefinesSolverInfo("UPWINDTYPE"),
"No UPWINDTYPE defined in session.");
......@@ -118,7 +118,7 @@ namespace Nektar
{
ASSERTL0(false,"No Continuous Galerkin implemented for compressible flow solver.");
}
m_session->LoadParameter("Gamma",m_gamma,1.4);
m_session->LoadParameter("GasConstant",m_GasConstant,287.058);
m_session->LoadParameter("CFLParameter",m_cfl,0.0);
......
......@@ -830,7 +830,8 @@ namespace Nektar
fields[i] = m_fields[i]->UpdateCoeffs();
}
if(m_projectionType==MultiRegions::eGalerkin)
if((m_projectionType==MultiRegions::eGalerkin)||
(m_projectionType==MultiRegions::eMixed_CG_Discontinuous))
{
// calculate the variable u* = Mu
// we are going to TimeIntegrate this new variable u*
......@@ -960,7 +961,8 @@ namespace Nektar
m_time += m_timestep;
if(m_projectionType==MultiRegions::eGalerkin)
if((m_projectionType==MultiRegions::eGalerkin)||
(m_projectionType==MultiRegions::eMixed_CG_Discontinuous))
{
// Project the solution u* onto the boundary conditions to
// obtain the actual solution
......
......@@ -109,7 +109,7 @@ namespace Nektar
// Load parameter alpha.
m_session->LoadParameter("Alpha", m_alpha);
ASSERTL0(m_projectionType == MultiRegions::eDiscontinuousGalerkin,
ASSERTL0(m_projectionType == MultiRegions::eDiscontinuous,
"CG not implemented yet.");
// Set up storage arrays.
......@@ -193,7 +193,7 @@ namespace Nektar
switch(m_projectionType)
{
case MultiRegions::eDiscontinuousGalerkin:
case MultiRegions::eDiscontinuous:
{
// Just copy over array
int npoints = GetNpoints();
......
......@@ -165,7 +165,7 @@ namespace Nektar
}
else if(ProjectStr == "DisContinuous")
{
m_projectionType = MultiRegions::eDiscontinuousGalerkin;
m_projectionType = MultiRegions::eDiscontinuous;
}
else
{
......
......@@ -85,9 +85,13 @@ namespace Nektar
{
m_projectionType = MultiRegions::eGalerkin;
}
else if((ProjectStr == "MixedCGDG")||(ProjectStr == "Mixed_CG_Discontinuous"))
{
m_projectionType = MultiRegions::eMixed_CG_Discontinuous;
}
else if(ProjectStr == "DisContinuous")
{
m_projectionType = MultiRegions::eDiscontinuousGalerkin;
m_projectionType = MultiRegions::eDiscontinuous;
}
else
{
......@@ -122,7 +126,7 @@ namespace Nektar
const Array<OneD, int> &vel_loc,
const Array<OneD, const Array<OneD, NekDouble> > &pInarray,
Array<OneD, Array<OneD, NekDouble> > &pOutarray,
NekDouble m_time,
NekDouble time,
Array<OneD, NekDouble> &pWk)
{
int i,j;
......@@ -131,7 +135,7 @@ namespace Nektar
Array<OneD, Array<OneD, NekDouble> > velocity(VelDim);
Array<OneD, NekDouble > Deriv;
m_nConvectiveFields = nConvectiveFields;
ASSERTL1(nConvectiveFields == pInarray.num_elements(),"Number of convective fields and Inarray are not compatible");
for(i = 0; i < VelDim; ++i)
{
......@@ -147,6 +151,25 @@ namespace Nektar
}
}
DoAdvection(pFields,velocity,pInarray,pOutarray,time,pWk);
}
void AdvectionTerm::DoAdvection(Array<OneD, MultiRegions::ExpListSharedPtr> &pFields,
const Array<OneD, const Array<OneD, NekDouble> > &velocity,
const Array<OneD, const Array<OneD, NekDouble> > &pInarray,
Array<OneD, Array<OneD, NekDouble> > &pOutarray,
NekDouble time,
Array<OneD, NekDouble> &pWk)
{
int i;
int nqtot = pFields[0]->GetTotPoints();
int VelDim = velocity.num_elements();
Array<OneD, NekDouble > Deriv;
m_nConvectiveFields = pInarray.num_elements();
// Set up Derivative work space;
if(pWk.num_elements())
{
......@@ -161,7 +184,7 @@ namespace Nektar
for(i=0; i< m_nConvectiveFields; ++i)