Commit 56e4de5e authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Modifications for 3D Substepping in Inc NavierStokes Solver


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@4030 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent ca88815f
......@@ -450,8 +450,7 @@ namespace Nektar
*
* According to their boundary region, the separate segmental boundary
* expansions are bundled together in an object of the class
* MultiRegions#ExpList1D. The list of expansions of the Dirichlet
* boundary regions are listed first in the array #m_bndCondExpansions.
* MultiRegions#ExpList1D.
*
* \param graph2D A mesh, containing information about the domain and
* the spectral/hp element expansion.
......
This diff is collapsed.
......@@ -128,16 +128,9 @@ namespace Nektar
const SpatialDomains::BoundaryConditions &bcs,
const std::string &variable);
bool SameTypeOfBoundaryConditions(const DisContField3D &In);
void SetUpDG();
/// Populates the list of boundary condition expansions.
void SetBoundaryConditionExpansion(
const SpatialDomains::MeshGraphSharedPtr &graph3D,
const SpatialDomains::BoundaryConditions &bcs,
const std::string variable,
Array<OneD, ExpListSharedPtr> &bndCondExpansions,
Array<OneD, SpatialDomains::BoundaryConditionShPtr>
&bndConditions);
bool SameTypeOfBoundaryConditions(const DisContField3D &In);
/// Generates a map of periodic faces in the mesh.
void GetPeriodicFaces(
......@@ -167,31 +160,8 @@ namespace Nektar
const NekDouble x2_in = NekConstants::kNekUnsetDouble,
const NekDouble x3_in = NekConstants::kNekUnsetDouble);
/**
* \brief This method extracts the "forward" and "backward" trace
* data from the array \a field and puts the data into output
* vectors \a Fwd and \a Bwd.
*
* We first define the convention which defines "forwards" and
* "backwards". First an association is made between the face of
* each element and its corresponding face in the trace space
* using the mapping #m_traceMap. The element can either be
* left-adjacent or right-adjacent to this trace face (see
* Expansion2D::GetLeftAdjacentElementExp). Boundary faces are
* always left-adjacent since left-adjacency is populated first.
*
* If the element is left-adjacent we extract the face trace data
* from \a field into the forward trace space \a Fwd; otherwise,
* we place it in the backwards trace space \a Bwd. In this way,
* we form a unique set of trace normals since these are always
* extracted from left-adjacent elements.
*
* \param field is a NekDouble array which contains the 3D data
* from which we wish to extract the backward and forward
* orientated trace/face arrays.
*
* \return Updates a NekDouble array \a Fwd and \a Bwd
*/
bool IsLeftAdjacentFace(const int n, const int e);
virtual void v_GetFwdBwdTracePhys(
Array<OneD,NekDouble> &Fwd,
Array<OneD,NekDouble> &Bwd);
......@@ -207,6 +177,10 @@ namespace Nektar
virtual void v_AddTraceIntegral(
const Array<OneD, const NekDouble> &Fn,
Array<OneD, NekDouble> &outarray);
virtual void v_AddFwdBwdTraceIntegral(
const Array<OneD, const NekDouble> &Fwd,
const Array<OneD, const NekDouble> &Bwd,
Array<OneD, NekDouble> &outarray);
virtual const Array<OneD, const MultiRegions::ExpListSharedPtr>
&v_GetBndCondExpansions();
......@@ -220,6 +194,10 @@ namespace Nektar
virtual ExpListSharedPtr &v_GetTrace(void)
{
if(m_trace == NullExpListSharedPtr)
{
SetUpDG();
}
return m_trace;
}
......
......@@ -300,9 +300,9 @@ namespace Nektar
//m_fields[i]->SetWaveSpace(false);
}
}
}
// Half mode stability analysis
else if(m_HalfMode)
else if(m_HalfMode)
{
const LibUtilities::PointsKey PkeyZ(m_npointsZ, LibUtilities::eFourierSingleModeSpaced);
......@@ -398,7 +398,18 @@ namespace Nektar
::AllocateSharedPtr(*firstfield, m_graph,
m_session->GetVariable(i));
}
//}
if(m_projectionType == MultiRegions::eMixed_CG_Discontinuous)
{
/// Setting up the normals
m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(m_spacedim);
for(i = 0; i < m_spacedim; ++i)
{
m_traceNormals[i] = Array<OneD, NekDouble> (GetTraceNpoints());
}
m_fields[0]->GetTrace()->GetNormals(m_traceNormals);
}
break;
}
default:
......
......@@ -869,63 +869,111 @@ namespace Nektar
NekDouble UnsteadySystem::v_GetTimeStep(int ExpOrder, NekDouble CFL, NekDouble TimeStability)
{
ASSERTL0(false, "v_GetTimeStep is not implemented in the base class (UnsteadySystem). Check if your equation class has its own implementation");
return 0.0;
}
Array<OneD,NekDouble> UnsteadySystem::GetStdVelocity(const Array<OneD, Array<OneD,NekDouble> > inarray)
{
// Checking if the problem is 2D
ASSERTL0(m_expdim==2,"Method not implemented for 1D and 3D");
ASSERTL0(m_expdim>=2,"Method not implemented for 1D");
int nTotQuadPoints = GetTotPoints();
int n_element = m_fields[0]->GetExpSize(); // number of element in the mesh
int nvel = inarray.num_elements();
int npts = 0;
NekDouble pntVelocity;
// Getting the standard velocity vector on the 2D normal space
Array<OneD, Array<OneD, NekDouble> > stdVelocity(2);
Array<OneD, Array<OneD, NekDouble> > stdVelocity(nvel);
Array<OneD, NekDouble> stdV(n_element,0.0);
for (int i = 0; i < 2; ++i)
for (int i = 0; i < nvel; ++i)
{
stdVelocity[i] = Array<OneD, NekDouble>(nTotQuadPoints);
}
for(int el = 0; el < n_element; ++el)
{
Array<OneD, const NekDouble> jac = m_fields[0]->GetExp(el)->GetGeom2D()->GetJac();
Array<TwoD, const NekDouble> gmat = m_fields[0]->GetExp(el)->GetGeom2D()->GetGmat();
int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
if(m_fields[0]->GetExp(el)->GetGeom2D()->GetGtype() == SpatialDomains::eDeformed)
{
// d xi/ dx = gmat = 1/J * d x/d xi
for(int i=0; i<n_points; i++)
if(nvel == 2)
{
for(int el = 0; el < n_element; ++el)
{
int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
Array<OneD, const NekDouble> jac = m_fields[0]->GetExp(el)->GetGeom2D()->GetJac();
Array<TwoD, const NekDouble> gmat = m_fields[0]->GetExp(el)->GetGeom2D()->GetGmat();
if(m_fields[0]->GetExp(el)->GetGeom2D()->GetGtype() == SpatialDomains::eDeformed)
{
stdVelocity[0][i] = gmat[0][i]*inarray[0][i] + gmat[2][i]*inarray[1][i];
stdVelocity[1][i] = gmat[1][i]*inarray[0][i] + gmat[3][i]*inarray[1][i];
for(int i=0; i<n_points; i++)
{
stdVelocity[0][i] = gmat[0][i]*inarray[0][i] + gmat[2][i]*inarray[1][i];
stdVelocity[1][i] = gmat[1][i]*inarray[0][i] + gmat[3][i]*inarray[1][i];
}
}
}
else
{
else
{
for(int i=0; i<n_points; i++)
{
stdVelocity[0][i] = gmat[0][0]*inarray[0][i] + gmat[2][0]*inarray[1][i];
stdVelocity[1][i] = gmat[1][0]*inarray[0][i] + gmat[3][0]*inarray[1][i];
}
}
for(int i=0; i<n_points; i++)
{
stdVelocity[0][i] = gmat[0][0]*inarray[0][i] + gmat[2][0]*inarray[1][i];
stdVelocity[1][i] = gmat[1][0]*inarray[0][i] + gmat[3][0]*inarray[1][i];
pntVelocity = sqrt(stdVelocity[0][i]*stdVelocity[0][i] + stdVelocity[1][i]*stdVelocity[1][i]);
if(pntVelocity>stdV[el])
{
stdV[el] = pntVelocity;
}
}
}
NekDouble pntVelocity;
for(int i=0; i<n_points; i++)
{
pntVelocity = sqrt(stdVelocity[0][i]*stdVelocity[0][i] + stdVelocity[1][i]*stdVelocity[1][i]);
if(pntVelocity>stdV[el])
}
else
{
for(int el = 0; el < n_element; ++el)
{
int n_points = m_fields[0]->GetExp(el)->GetTotPoints();
Array<OneD, const NekDouble> jac = m_fields[0]->GetExp(el)->GetGeom3D()->GetJac();
Array<TwoD, const NekDouble> gmat = m_fields[0]->GetExp(el)->GetGeom3D()->GetGmat();
if(m_fields[0]->GetExp(el)->GetGeom3D()->GetGtype() == SpatialDomains::eDeformed)
{
stdV[el] = pntVelocity;
for(int i=0; i<n_points; i++)
{
stdVelocity[0][i] = gmat[0][i]*inarray[0][i] + gmat[3][i]*inarray[1][i] + gmat[6][i]*inarray[2][i];
stdVelocity[1][i] = gmat[1][i]*inarray[0][i] + gmat[4][i]*inarray[1][i] + gmat[7][i]*inarray[2][i];
stdVelocity[2][i] = gmat[2][i]*inarray[0][i] + gmat[5][i]*inarray[1][i] + gmat[8][i]*inarray[2][i];
}
}
else
{
Array<OneD, const NekDouble> jac = m_fields[0]->GetExp(el)->GetGeom3D()->GetJac();
Array<TwoD, const NekDouble> gmat = m_fields[0]->GetExp(el)->GetGeom3D()->GetGmat();
for(int i=0; i<n_points; i++)
{
stdVelocity[0][i] = gmat[0][0]*inarray[0][i] + gmat[3][0]*inarray[1][i] + gmat[6][0]*inarray[2][i];
stdVelocity[1][i] = gmat[1][0]*inarray[0][i] + gmat[4][0]*inarray[1][i] + gmat[7][0]*inarray[2][i];
stdVelocity[2][i] = gmat[2][0]*inarray[0][i] + gmat[5][0]*inarray[1][i] + gmat[8][0]*inarray[2][i];
}
}
for(int i=0; i<n_points; i++)
{
pntVelocity = sqrt(stdVelocity[0][i]*stdVelocity[0][i] + stdVelocity[1][i]*stdVelocity[1][i] + stdVelocity[2][i]*stdVelocity[2][i]);
if(pntVelocity>stdV[el])
{
stdV[el] = pntVelocity;
}
}
}
}
......
......@@ -1088,6 +1088,23 @@ namespace Nektar
void VelocityCorrectionScheme::AddDuDt(const Array<OneD, const Array<OneD, NekDouble> > &N, NekDouble Aii_Dt)
{
switch(m_velocity.num_elements())
{
case 1:
ASSERTL0(false,"Velocity correction scheme not designed to have just one velocity component");
break;
case 2:
AddDuDt2D(N,Aii_Dt);
break;
case 3:
AddDuDt3D(N,Aii_Dt);
break;
}
}
void VelocityCorrectionScheme::AddDuDt2D(const Array<OneD, const Array<OneD, NekDouble> > &N, NekDouble Aii_Dt)
{
int i,n;
ASSERTL0(m_velocity.num_elements() == 2," Routine currently only set up for 2D");
......@@ -1170,7 +1187,7 @@ namespace Nektar
Blas::Dscal(nq,1.0/Aii_Dt,&ubc[0],1);
Blas::Dscal(nq,1.0/Aii_Dt,&vbc[0],1);
// subtrace off du/dt derivative (can ignore aii_dt for moment)
// subtrace off du/dt derivative
Pbc->NormVectorIProductWRTBase(ubc,vbc,Pvals);
Vmath::Vsub(ncoeffs,Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),1, Pvals,1, Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),1);
......@@ -1185,9 +1202,114 @@ namespace Nektar
{
ASSERTL0(false,"Unknown USERDEFINEDTYPE in pressure boundary condition");
}
}
}
}
void VelocityCorrectionScheme::AddDuDt3D(const Array<OneD, const Array<OneD, NekDouble> > &N, NekDouble Aii_Dt)
{
int i,n;
ASSERTL0(m_velocity.num_elements() == 3," Routine currently only set up for 3D");
Array<OneD, const SpatialDomains::BoundaryConditionShPtr > PBndConds;
Array<OneD, MultiRegions::ExpListSharedPtr> PBndExp,UBndExp,VBndExp,WBndExp;
PBndConds = m_pressure->GetBndConditions();
PBndExp = m_pressure->GetBndCondExpansions();
UBndExp = m_fields[m_velocity[0]]->GetBndCondExpansions();
VBndExp = m_fields[m_velocity[1]]->GetBndCondExpansions();
WBndExp = m_fields[m_velocity[2]]->GetBndCondExpansions();
StdRegions::StdExpansionSharedPtr elmt;
StdRegions::StdExpansion2DSharedPtr Pbc;
int maxpts = 0;
int cnt,elmtid,nq,offset, boundary,ncoeffs;
// find the maximum values of points
for(n = 0; n < PBndConds.num_elements(); ++n)
{
for(i = 0; i < PBndExp[n]->GetExpSize(); ++i)
{
maxpts = max(maxpts, PBndExp[n]->GetExp(i)->GetTotPoints());
}
}
Array<OneD, NekDouble> N1(maxpts), N2(maxpts), N3(maxpts);
Array<OneD, NekDouble> ubc(maxpts),vbc(maxpts),wbc(maxpts);
Array<OneD, NekDouble> Pvals(maxpts);
Array<OneD, NekDouble> Nu,Nv,Nw,Ptmp;
for(cnt = n = 0; n < PBndConds.num_elements(); ++n)
{
SpatialDomains::BndUserDefinedType type = PBndConds[n]->GetUserDefined();
if(type == SpatialDomains::eHigh)
{
for(i = 0; i < PBndExp[n]->GetExpSize(); ++i,cnt++)
{
// find element and face of this expansion.
elmtid = m_pressureBCtoElmtID[cnt];
elmt = m_fields[0]->GetExp(elmtid);
offset = m_fields[0]->GetPhys_Offset(elmtid);
Nu = N[0] + offset;
Nv = N[1] + offset;
Nw = N[2] + offset;
Pbc = boost::dynamic_pointer_cast<StdRegions::StdExpansion2D> (PBndExp[n]->GetExp(i));
nq = Pbc->GetTotPoints();
ncoeffs = Pbc->GetNcoeffs();
boundary = m_pressureBCtoTraceID[cnt];
// Get velocity bc
UBndExp[n]->GetExp(i)->BwdTrans(UBndExp[n]->GetCoeffs() + UBndExp[n]->GetCoeff_Offset(i),ubc);
VBndExp[n]->GetExp(i)->BwdTrans(VBndExp[n]->GetCoeffs() + VBndExp[n]->GetCoeff_Offset(i),vbc);
WBndExp[n]->GetExp(i)->BwdTrans(WBndExp[n]->GetCoeffs() + WBndExp[n]->GetCoeff_Offset(i),wbc);
// Get edge values and put into Nu,Nv
elmt->GetFacePhysVals(boundary,Pbc,Nu,N1);
elmt->GetFacePhysVals(boundary,Pbc,Nv,N2);
elmt->GetFacePhysVals(boundary,Pbc,Nw,N3);
// Take different as Forward Euler but N1,N2
// actually contain the integration of the
// previous steps from the time integration
// scheme.
Vmath::Vsub(nq,ubc,1,N1,1,ubc,1);
Vmath::Vsub(nq,vbc,1,N2,1,vbc,1);
Vmath::Vsub(nq,wbc,1,N3,1,wbc,1);
// Divide by aii_Dt to get correct Du/Dt. This is
// because all coefficients in the integration
// scheme are normalised so u^{n+1} has unit
// coefficient and N is already multiplied by
// local coefficient when taken from integration
// scheme
Blas::Dscal(nq,1.0/Aii_Dt,&ubc[0],1);
Blas::Dscal(nq,1.0/Aii_Dt,&vbc[0],1);
Blas::Dscal(nq,1.0/Aii_Dt,&wbc[0],1);
// subtrace off du/dt derivative
Pbc->NormVectorIProductWRTBase(ubc,vbc,wbc,Pvals);
Vmath::Vsub(ncoeffs,Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),1, Pvals,1, Ptmp = PBndExp[n]->UpdateCoeffs()+PBndExp[n]->GetCoeff_Offset(i),1);
}
}
// setting if just standard BC no High order
else if(type == SpatialDomains::eNoUserDefined || type == SpatialDomains::eTimeDependent)
{
cnt += PBndExp[n]->GetExpSize();
}
else
{
ASSERTL0(false,"Unknown USERDEFINEDTYPE in pressure boundary condition");
}
}
}
// Evaluate divergence of velocity field.
void VelocityCorrectionScheme::SetUpPressureForcing(const Array<OneD, const Array<OneD, NekDouble> > &fields, Array<OneD, Array<OneD, NekDouble> > &Forcing, const NekDouble aii_Dt)
......@@ -1238,7 +1360,7 @@ namespace Nektar
}
}
void VelocityCorrectionScheme::FillHOPBCMap(const int HOPBCnumber)
void VelocityCorrectionScheme::FillHOPBCMap(const int HOPBCnumber)
{
// Count number of HBC conditions
......
......@@ -96,7 +96,9 @@ namespace Nektar
void SelectiveFrequencyDamping(void);
void AddDuDt(const Array<OneD, const Array<OneD, NekDouble> > &N, NekDouble Aii_Dt);
void AddDuDt (const Array<OneD, const Array<OneD, NekDouble> > &N, NekDouble Aii_Dt);
void AddDuDt2D(const Array<OneD, const Array<OneD, NekDouble> > &N, NekDouble Aii_Dt);
void AddDuDt3D(const Array<OneD, const Array<OneD, NekDouble> > &N, NekDouble Aii_Dt);
protected:
......
Markdown is supported
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