Skip to content

GitLab

Commit 9cc907aa authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Updates related to substepping and getting VortexWaveInteractionSolver working again 


git-svn-id: https://gforge.sci.utah.edu/svn/nektar/trunk@3989 305cdda6-5ce1-45b3-a98d-dfc68c8b3305
parent dfb00652
Hide whitespace changes
Inline Side-by-side
library/LibUtilities/Interpreter/AnalyticExpressionEvaluator.cpp
View file @ 9cc907aa
......@@ -343,7 +343,13 @@ namespace Nektar
}
else
{
ASSERTL1(m_constant[it->second] == value, "Attempt to add numerically different constants under the same name: " + name);
if(m_constant[it->second] != value)
{
std::string errormsg("Attempt to add numerically different constants under the same name: ");
errormsg += name;
std::cout << errormsg << std::endl;
}
//ASSERTL1(m_constant[it->second] == value, "Attempt to add numerically different constants under the same name: " + name);
}
return it->second;
}
......
library/LibUtilities/TimeIntegration/TimeIntegrationScheme.cpp
View file @ 9cc907aa
......@@ -400,12 +400,12 @@ namespace Nektar
m_timeLevelOffset[1] = 0;
}
break;
case eBDFImplicitOrder2:
case eBDFImplicitOrder2:
{
NekDouble third = 1.0/3.0;
m_numsteps = 2;
m_numsteps = 2;
m_numstages = 1;
m_A = Array<OneD, Array<TwoD,NekDouble> >(1);
m_B = Array<OneD, Array<TwoD,NekDouble> >(1);
......@@ -678,7 +678,7 @@ namespace Nektar
m_timeLevelOffset[0] = 0;
}
break;
case eCNAB:
case eCNAB:
{
NekDouble secondth = 1.0/2.0;
m_numsteps = 4;
......@@ -1097,7 +1097,7 @@ namespace Nektar
// initial value. Initialise all other multi-step values
// and derivatives to zero
TimeIntegrationSolutionSharedPtr y_out =
MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(m_schemeKey,y_0,time,timestep);
MemoryManager<TimeIntegrationSolution>::AllocateSharedPtr(m_schemeKey,y_0,time,timestep);
// calculate the initial derivative, if is part of the
// solution vector of the current scheme
......@@ -1126,7 +1126,7 @@ namespace Nektar
}
return y_out;
}
}
TimeIntegrationScheme::ConstDoubleArray&
TimeIntegrationScheme::TimeIntegrate(const NekDouble timestep,
......@@ -1217,8 +1217,7 @@ namespace Nektar
// STEP 2: time-integrate for one step using the
// current scheme
TimeIntegrationSolutionSharedPtr solvector_out = MemoryManager<TimeIntegrationSolution>::
AllocateSharedPtr(GetIntegrationSchemeKey(),nvar,npoints); // output solution vector of the current scheme
TimeIntegrationSolutionSharedPtr solvector_out = MemoryManager<TimeIntegrationSolution>:: AllocateSharedPtr(GetIntegrationSchemeKey(),nvar,npoints); // output solution vector of the current scheme
// integrate
TimeIntegrate(timestep, solvector_in->GetSolutionVector(),
......@@ -1320,7 +1319,7 @@ namespace Nektar
// solution vector of the current scheme
//DoubleArray& y_n = solvector_out->GetValue ( curTimeLevels[n] );
//NekDouble t_n = solvector_out->GetValueTime( curTimeLevels[n] );
y_n = solvector_out->GetValue ( curTimeLevels[n] );
y_n = solvector_out->GetValue ( curTimeLevels[n] );
t_n = solvector_out->GetValueTime( curTimeLevels[n] );
// Set the calculated value in the master solution vector
......@@ -1331,8 +1330,9 @@ namespace Nektar
{
// Get the calculated derivative out of the output
// solution vector of the current scheme
//DoubleArray& dtFy_n = solvector_out->GetDerivative ( curTimeLevels[n] );
dtFy_n = solvector_out->GetDerivative ( curTimeLevels[n] );
// DoubleArray& dtFy_n =
// solvector_out->GetDerivative (curTimeLevels[n]);
dtFy_n = solvector_out->GetDerivative ( curTimeLevels[n] );
// Set the calculated derivative in the master
// solution vector
......@@ -1462,18 +1462,6 @@ namespace Nektar
// The stage values Y are a linear combination of:
// 1: the stage derivatives
/* if( i == 0 )
{
#if 1
// SJS Why is tmp memory being redefined here
// it seems it has been defined above.
for (k=0; k<nvar; k++)
{
tmp[k] = Array<OneD, NekDouble>(npoints,0.0);
}
#endif
}*/
if( i != 0 )
{
for(k = 0; k < nvar; k++)
......@@ -1510,7 +1498,7 @@ namespace Nektar
// 2: the imported multi-step solution of the
// previous time level
for( j = 0; j < m_numsteps; j++)
for(j = 0; j < m_numsteps; j++)
{
for(k = 0; k < nvar; k++)
{
......@@ -1525,20 +1513,20 @@ namespace Nektar
if(type == eDiagonallyImplicit)
{
if(m_numstages==1)
{
T= t_old[0]+timestep;
}
else
{
T= t_old[0];
for(int j=0; j<=i; ++j)
{
T += A(i,j)*timestep;
}
}
{
T= t_old[0]+timestep;
}
else
{
T= t_old[0];
for(int j=0; j<=i; ++j)
{
T += A(i,j)*timestep;
}
}
op.DoImplicitSolve(tmp, Y, T, A(i,i)*timestep);
for(k = 0; k < nvar; k++)
{
Vmath::Vsub(npoints,Y[k],1,tmp[k],1,F[i][k],1);
......@@ -1548,17 +1536,17 @@ namespace Nektar
else if(type == eIMEX)
{
if(m_numstages==1)
{
T= t_old[0]+timestep;
}
else
{
T= t_old[0];
for(int j=0; j<=i; ++j)
{
T += A(i,j)*timestep;
}
}
{
T= t_old[0]+timestep;
}
else
{
T= t_old[0];
for(int j=0; j<=i; ++j)
{
T += A(i,j)*timestep;
}
}
if(fabs(A(i,i)) > NekConstants::kNekZeroTol)
{
......@@ -1602,22 +1590,22 @@ namespace Nektar
}
if (m_numstages==1 && type == eIMEX)
{
t_new[0] = t_old[0]+timestep;
}
else
{
t_new[0] = B(0,0)*timestep;
for(j = 1; j < m_numstages; j++)
{
t_new[0] += B(0,j)*timestep;
}
for(j = 0; j < m_numsteps; j++)
{
t_new[0] += V(0,j)*t_old[j];
}
i_start = 1;
}
{
t_new[0] = t_old[0]+timestep;
}
else
{
t_new[0] = B(0,0)*timestep;
for(j = 1; j < m_numstages; j++)
{
t_new[0] += B(0,j)*timestep;
}
for(j = 0; j < m_numsteps; j++)
{
t_new[0] += V(0,j)*t_old[j];
}
i_start = 1;
}
}
for(i = i_start; i < m_numsteps; i++)
......@@ -1637,12 +1625,12 @@ namespace Nektar
y_new[i][k],1);
}
}
if(m_numstages != 1 || type != eIMEX)
{
t_new[i] = B(i,0)*timestep;
}
if(m_numstages != 1 || type != eIMEX)
{
t_new[i] = B(i,0)*timestep;
}
for(j = 1; j < m_numstages; j++)
{
for(k = 0; k < nvar; k++)
......@@ -1657,10 +1645,10 @@ namespace Nektar
y_new[i][k],1);
}
}
if(m_numstages != 1 || type != eIMEX)
{
t_new[i] += B(i,j)*timestep;
}
if(m_numstages != 1 || type != eIMEX)
{
t_new[i] += B(i,j)*timestep;
}
}
// 2: the imported multi-step solution of the previous
......@@ -1672,11 +1660,10 @@ namespace Nektar
Vmath::Svtvp(npoints,V(i,j),y_old[j][k],1,
y_new[i][k],1,y_new[i][k],1);
}
if(m_numstages != 1 || type != eIMEX)
{
t_new[i] += V(i,j)*t_old[j];
}
if(m_numstages != 1 || type != eIMEX)
{
t_new[i] += V(i,j)*t_old[j];
}
}
}
......
library/LibUtilities/TimeIntegration/TimeIntegrationScheme.h
View file @ 9cc907aa
......@@ -812,6 +812,12 @@ namespace Nektar
}
}
// sets the soln Vector
inline void SetSolVector(const int Offset, const DoubleArray& y)
{
m_solVector[Offset] = y;
}
// Rotate the solution vector
// (i.e. updating without calculating/inserting new values)
inline void RotateSolutionVector()
......
library/LocalRegions/QuadExp.cpp
View file @ 9cc907aa
......@@ -1241,7 +1241,6 @@ namespace Nektar
DNekMatSharedPtr QuadExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
{
DNekMatSharedPtr returnval;
switch(mkey.GetMatrixType())
{
case StdRegions::eHybridDGHelmholtz:
......@@ -1255,7 +1254,6 @@ namespace Nektar
default:
returnval = StdQuadExp::v_GenMatrix(mkey);
}
return returnval;
}
......
library/LocalRegions/TriExp.cpp
View file @ 9cc907aa
......@@ -1185,7 +1185,6 @@ namespace Nektar
DNekMatSharedPtr TriExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
{
DNekMatSharedPtr returnval;
switch(mkey.GetMatrixType())
{
case StdRegions::eHybridDGHelmholtz:
......@@ -1200,6 +1199,7 @@ namespace Nektar
returnval = StdTriExp::v_GenMatrix(mkey);
break;
}
return returnval;
}
......
library/MultiRegions/DisContField2D.cpp
View file @ 9cc907aa
......@@ -44,7 +44,8 @@ namespace Nektar
DisContField2D::DisContField2D(void):
ExpList2D(),
m_bndCondExpansions(),
m_bndConditions()
m_bndConditions(),
m_trace(NullExpListSharedPtr)
{
}
......@@ -62,7 +63,6 @@ namespace Nektar
m_boundaryEdges (In.m_boundaryEdges),
m_perEdgeToExpMap (In.m_perEdgeToExpMap)
{
}
DisContField2D::DisContField2D(
......@@ -73,7 +73,8 @@ namespace Nektar
const bool DeclareCoeffPhysArrays):
ExpList2D(pSession,graph2D,DeclareCoeffPhysArrays,variable),
m_bndCondExpansions(),
m_bndConditions()
m_bndConditions(),
m_trace(NullExpListSharedPtr)
{
SpatialDomains::BoundaryConditions bcs(m_session, graph2D);
......@@ -138,7 +139,8 @@ namespace Nektar
const std::string &variable,
const bool SetUpJustDG,
const bool DeclareCoeffPhysArrays) :
ExpList2D(In,DeclareCoeffPhysArrays)
ExpList2D(In,DeclareCoeffPhysArrays),
m_trace(NullExpListSharedPtr)
{
// Set up boundary conditions for this variable.
SpatialDomains::BoundaryConditions bcs(m_session, graph2D);
......@@ -269,101 +271,102 @@ namespace Nektar
*/
void DisContField2D::SetUpDG()
{
ExpList1DSharedPtr trace;
SpatialDomains::MeshGraph2DSharedPtr graph2D =
boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(
m_graph);
// Set up matrix map
m_globalBndMat = MemoryManager<GlobalLinSysMap>::
AllocateSharedPtr();
// Set up Trace space
trace = MemoryManager<ExpList1D>::AllocateSharedPtr(
m_bndCondExpansions, m_bndConditions, *m_exp,
graph2D, m_periodicEdges);
m_trace = boost::dynamic_pointer_cast<ExpList>(trace);
m_traceMap = MemoryManager<AssemblyMapDG>::
AllocateSharedPtr(m_session, graph2D, trace, *this,
m_bndCondExpansions, m_bndConditions,
m_periodicEdges);
Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
// Scatter trace segments to 2D elements. For each element, we find
// the trace segment associated to each edge. The element then
// retains a pointer to the trace space segments, to ensure
// uniqueness of normals when retrieving from two adjoining elements
// which do not lie in a plane.
SpatialDomains::Geometry1DSharedPtr ElmtSegGeom;
SpatialDomains::Geometry1DSharedPtr TraceSegGeom;
for (int i = 0; i < m_exp->size(); ++i)
if(m_trace == NullExpListSharedPtr) // check for multiple calls
{
for (int j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
ExpList1DSharedPtr trace;
SpatialDomains::MeshGraph2DSharedPtr graph2D =
boost::dynamic_pointer_cast<SpatialDomains::MeshGraph2D>(
m_graph);
// Set up matrix map
m_globalBndMat = MemoryManager<GlobalLinSysMap>::
AllocateSharedPtr();
// Set up Trace space
trace = MemoryManager<ExpList1D>::AllocateSharedPtr(m_bndCondExpansions, m_bndConditions, *m_exp, graph2D, m_periodicEdges);
m_trace = boost::dynamic_pointer_cast<ExpList>(trace);
m_traceMap = MemoryManager<AssemblyMapDG>::
AllocateSharedPtr(m_session, graph2D, trace, *this,
m_bndCondExpansions, m_bndConditions,
m_periodicEdges);
Array<OneD, Array<OneD, StdRegions::StdExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
// Scatter trace segments to 2D elements. For each element, we find
// the trace segment associated to each edge. The element then
// retains a pointer to the trace space segments, to ensure
// uniqueness of normals when retrieving from two adjoining elements
// which do not lie in a plane.
SpatialDomains::Geometry1DSharedPtr ElmtSegGeom;
SpatialDomains::Geometry1DSharedPtr TraceSegGeom;
for (int i = 0; i < m_exp->size(); ++i)
{
LocalRegions::Expansion2DSharedPtr exp2d =
boost::dynamic_pointer_cast<
LocalRegions::Expansion2D>((*m_exp)[i]);
LocalRegions::Expansion1DSharedPtr exp1d =
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(elmtToTrace[i][j]);
LocalRegions::ExpansionSharedPtr exp =
boost::dynamic_pointer_cast<
LocalRegions::Expansion> (elmtToTrace[i][j]);
exp2d->SetEdgeExp (j, exp );
exp1d->SetAdjacentElementExp(j, exp2d);
for (int j = 0; j < (*m_exp)[i]->GetNedges(); ++j)
{
LocalRegions::Expansion2DSharedPtr exp2d =
boost::dynamic_pointer_cast<
LocalRegions::Expansion2D>((*m_exp)[i]);
LocalRegions::Expansion1DSharedPtr exp1d =
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(elmtToTrace[i][j]);
LocalRegions::ExpansionSharedPtr exp =
boost::dynamic_pointer_cast<
LocalRegions::Expansion> (elmtToTrace[i][j]);
exp2d->SetEdgeExp (j, exp );
exp1d->SetAdjacentElementExp(j, exp2d);
}
}
}
// Set up physical normals
SetUpPhysNormals();
// Set up information for parallel jobs.
for (int i = 0; i < m_trace->GetExpSize(); ++i)
{
LocalRegions::Expansion1DSharedPtr traceEl =
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(m_trace->GetExp(i));
int offset = m_trace->GetPhys_Offset(i);
// Set up physical normals
SetUpPhysNormals();
if (m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
// Set up information for parallel jobs.
for (int i = 0; i < m_trace->GetExpSize(); ++i)
{
traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
traceEl->GetLeftAdjacentElementEdge());
LocalRegions::Expansion1DSharedPtr traceEl =
boost::dynamic_pointer_cast<
LocalRegions::Expansion1D>(m_trace->GetExp(i));
int offset = m_trace->GetPhys_Offset(i);
if (m_traceMap->GetTraceToUniversalMapUnique(offset) < 0)
{
traceEl->GetLeftAdjacentElementExp()->NegateEdgeNormal(
traceEl->GetLeftAdjacentElementEdge());
}
}
}
int cnt, n, e;
// Identify boundary edges
for(cnt = 0, n = 0; n < m_bndCondExpansions.num_elements(); ++n)
{
if (m_bndConditions[n]->GetBoundaryConditionType() !=
SpatialDomains::ePeriodic)
int cnt, n, e;
// Identify boundary edges
for(cnt = 0, n = 0; n < m_bndCondExpansions.num_elements(); ++n)
{
for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
if (m_bndConditions[n]->GetBoundaryConditionType() !=
SpatialDomains::ePeriodic)
{
m_boundaryEdges.insert(m_trace->GetOffset_Elmt_Id(
m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e)));
for(e = 0; e < m_bndCondExpansions[n]->GetExpSize(); ++e)
{
m_boundaryEdges.insert(m_trace->GetOffset_Elmt_Id(
m_traceMap->GetBndCondTraceToGlobalTraceMap(cnt+e)));
}
}
cnt += m_bndCondExpansions[n]->GetExpSize();
}
cnt += m_bndCondExpansions[n]->GetExpSize();
}
// Set up information for periodic boundary conditions.
for (n = 0; n < m_exp->size(); ++n)
{
for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
// Set up information for periodic boundary conditions.
for (n = 0; n < m_exp->size(); ++n)
{
map<int,int>::iterator it = m_periodicEdges.find(
(*m_exp)[n]->GetGeom2D()->GetEid(e));
if (it != m_periodicEdges.end())
for (e = 0; e < (*m_exp)[n]->GetNedges(); ++e)
{
m_perEdgeToExpMap[it->first] = make_pair(n, e);
map<int,int>::iterator it = m_periodicEdges.find(
(*m_exp)[n]->GetGeom2D()->GetEid(e));
if (it != m_periodicEdges.end())
{
m_perEdgeToExpMap[it->first] = make_pair(n, e);
}
}
}
}
......
library/MultiRegions/DisContField2D.h
View file @ 9cc907aa
......@@ -201,6 +201,11 @@ namespace Nektar
virtual ExpListSharedPtr &v_GetTrace()
{
if(m_trace == NullExpListSharedPtr)
{
SetUpDG();
}
return m_trace;
}
......
library/SolverUtils/EquationSystem.cpp
View file @ 9cc907aa
......@@ -114,8 +114,8 @@ namespace Nektar
m_useFFT = false;
m_dealiasing = false;
m_SingleMode = false;
m_HalfMode = false;
m_MultipleModes = false;
m_HalfMode = false;
m_MultipleModes = false;
m_HomogeneousType = eNotHomogeneous;
......@@ -648,7 +648,7 @@ namespace Nektar
}
EvaluateFunction(fieldStr, m_forces, "BodyForce");
if(m_SingleMode || m_HalfMode)
if(m_SingleMode || m_HalfMode)
{
for(int i=0; i< v_GetForceDimension(); ++i)
{
......@@ -799,16 +799,15 @@ namespace Nektar
cout << "- Field " << pFieldName << ": "
<< "from file " << filename << endl;
}
ImportFld(filename,m_fields);
/*
#if 0
ImportFld(filename,m_fields);
#else
std::vector<SpatialDomains::FieldDefinitionsSharedPtr> FieldDef;
std::vector<std::vector<NekDouble> > FieldData;
Array<OneD, NekDouble> vCoeffs(m_fields[0]->GetNcoeffs());
m_graph->Import(filename,FieldDef,FieldData);
int idx = -1;
// Loop over all the expansions
for(int i = 0; i < FieldDef.size(); ++i)
......@@ -822,13 +821,25 @@ namespace Nektar
}
}
ASSERTL1(idx >= 0, "Field " + pFieldName + " not found.");
m_fields[0]->ExtractDataToCoeffs(FieldDef[i], FieldData[i],