Commit 0805461d authored by Spencer Sherwin's avatar Spencer Sherwin Committed by Chris Cantwell
Browse files

Update for first attempt of filter history point and removed another transposition

parent f47c8616
......@@ -1869,7 +1869,7 @@ namespace Nektar
return sqrt(err);
}
Array<OneD, NekDouble> ExpList::v_HomogeneousEnergy (void)
Array<OneD, const NekDouble> ExpList::v_HomogeneousEnergy (void)
{
ASSERTL0(false,
"This method is not defined or valid for this class type");
......@@ -1887,7 +1887,7 @@ namespace Nektar
}
Array<OneD, unsigned int> ExpList::v_GetZIDs(void)
Array<OneD, const unsigned int> ExpList::v_GetZIDs(void)
{
ASSERTL0(false,
"This method is not defined or valid for this class type");
......@@ -1896,7 +1896,7 @@ namespace Nektar
return NoModes;
}
Array<OneD, unsigned int> ExpList::v_GetYIDs(void)
Array<OneD, const unsigned int> ExpList::v_GetYIDs(void)
{
ASSERTL0(false,
"This method is not defined or valid for this class type");
......
......@@ -470,7 +470,7 @@ namespace Nektar
/// This function calculates the energy associated with each one of the modes
/// of a 3D homogeneous nD expansion
Array<OneD, NekDouble> HomogeneousEnergy (void)
Array<OneD, const NekDouble> HomogeneousEnergy (void)
{
return v_HomogeneousEnergy();
}
......@@ -479,7 +479,7 @@ namespace Nektar
/// numbers in z-direction associated
/// with the 3D homogenous expansion. Required if a
/// parellelisation is applied in the Fourier direction
Array<OneD, unsigned int> GetZIDs(void)
Array<OneD, const unsigned int> GetZIDs(void)
{
return v_GetZIDs();
}
......@@ -495,7 +495,7 @@ namespace Nektar
/// numbers in y-direction associated
/// with the 3D homogenous expansion. Required if a
/// parellelisation is applied in the Fourier direction
Array<OneD, unsigned int> GetYIDs(void)
Array<OneD, const unsigned int> GetYIDs(void)
{
return v_GetYIDs();
}
......@@ -1165,10 +1165,10 @@ namespace Nektar
virtual NekDouble v_L2(void);
virtual NekDouble v_L2(const Array<OneD, const NekDouble> &soln);
virtual Array<OneD, NekDouble> v_HomogeneousEnergy(void);
virtual Array<OneD, const NekDouble> v_HomogeneousEnergy(void);
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void);
virtual Array<OneD, unsigned int> v_GetZIDs(void);
virtual Array<OneD, unsigned int> v_GetYIDs(void);
virtual Array<OneD, const unsigned int> v_GetZIDs(void);
virtual Array<OneD, const unsigned int> v_GetYIDs(void);
// 1D Scaling and projection
virtual void v_PhysInterp1DScaled(const NekDouble scale, const Array<OneD, NekDouble> &inarray, Array<OneD, NekDouble> &outarray);
......
......@@ -427,7 +427,7 @@ namespace Nektar
return sqrt(err);
}
Array<OneD, NekDouble> ExpList3DHomogeneous1D::v_HomogeneousEnergy(void)
Array<OneD, const NekDouble> ExpList3DHomogeneous1D::v_HomogeneousEnergy(void)
{
int cnt = 0, cnt1 = 0;
int ncoeffs_per_plane = m_planes[0]->GetNcoeffs();
......
......@@ -130,7 +130,7 @@ namespace Nektar
virtual NekDouble v_L2(void);
virtual NekDouble v_L2(const Array<OneD, const NekDouble> &soln);
virtual Array<OneD, NekDouble> v_HomogeneousEnergy(void);
virtual Array<OneD, const NekDouble> v_HomogeneousEnergy(void);
virtual void v_GetPeriodicEdges(
vector<map<int,int> > &periodicVertices,
......
......@@ -860,80 +860,83 @@ namespace Nektar
m_planes[i]->PhysDeriv(inarray + i*nP_pts ,tmp2 = out_d0 + i*nP_pts , tmp3 = out_d1 + i*nP_pts );
}
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
if(out_d2 != NullNekDouble1DArray)
{
if(m_WaveSpace)
{
temparray = inarray;
}
else
{
HomogeneousFwdTrans(inarray,temparray);
}
NekDouble sign = -1.0;
NekDouble beta;
//Half Mode
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe)
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierSingleMode ||
m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe || m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
beta = sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
if(m_WaveSpace)
{
temparray = inarray;
}
else
{
HomogeneousFwdTrans(inarray,temparray);
}
Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
}
else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
beta = -sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
NekDouble sign = -1.0;
NekDouble beta;
Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
}
//Fully complex
else
{
for(int i = 0; i < m_planes.num_elements(); i++)
//Half Mode
if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeRe)
{
beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
sign = -1.0*sign;
beta = sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
}
else if(m_homogeneousBasis->GetBasisType() == LibUtilities::eFourierHalfModeIm)
{
beta = -sign*2*M_PI*(m_transposition->GetK(0))/m_lhom;
Vmath::Smul(nP_pts,beta,temparray,1,outarray,1);
}
}
if(m_WaveSpace)
{
out_d2 = outarray;
}
else
{
HomogeneousBwdTrans(outarray,out_d2);
}
}
else
{
ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,"Parallelisation in the homogeneous direction implemented just for Fourier basis");
if(m_WaveSpace)
{
ASSERTL0(false,"Semi-phyisical time-stepping not implemented yet for non-Fourier basis");
//Fully complex
else
{
for(int i = 0; i < m_planes.num_elements(); i++)
{
beta = -sign*2*M_PI*(m_transposition->GetK(i))/m_lhom;
Vmath::Smul(nP_pts,beta,tmp1 = temparray + i*nP_pts,1,tmp2 = outarray + (i-int(sign))*nP_pts,1);
sign = -1.0*sign;
}
}
if(m_WaveSpace)
{
out_d2 = outarray;
}
else
{
HomogeneousBwdTrans(outarray,out_d2);
}
}
else
{
StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXYtoZ);
ASSERTL0(m_comm->GetColumnComm()->GetSize() == 1,"Parallelisation in the homogeneous direction implemented just for Fourier basis");
for(int i = 0; i < nP_pts; i++)
if(m_WaveSpace)
{
StdSeg.PhysDeriv(temparray + i*m_planes.num_elements(), tmp2 = outarray + i*m_planes.num_elements());
ASSERTL0(false,"Semi-phyisical time-stepping not implemented yet for non-Fourier basis");
}
else
{
StdRegions::StdSegExp StdSeg(m_homogeneousBasis->GetBasisKey());
m_transposition->Transpose(inarray,temparray,false,LibUtilities::eXYtoZ);
for(int i = 0; i < nP_pts; i++)
{
StdSeg.PhysDeriv(temparray + i*m_planes.num_elements(), tmp2 = outarray + i*m_planes.num_elements());
}
m_transposition->Transpose(outarray,out_d2,false,LibUtilities::eZtoXY);
Vmath::Smul(nT_pts,2.0/m_lhom,out_d2,1,out_d2,1);
}
m_transposition->Transpose(outarray,out_d2,false,LibUtilities::eZtoXY);
Vmath::Smul(nT_pts,2.0/m_lhom,out_d2,1,out_d2,1);
}
}
}
......@@ -1058,7 +1061,7 @@ namespace Nektar
return m_transposition;
}
Array<OneD, unsigned int> ExpListHomogeneous1D::v_GetZIDs(void)
Array<OneD, const unsigned int> ExpListHomogeneous1D::v_GetZIDs(void)
{
return m_transposition->GetPlanesIDs();
}
......
......@@ -244,7 +244,7 @@ namespace Nektar
virtual LibUtilities::TranspositionSharedPtr v_GetTransposition(void);
virtual Array<OneD, unsigned int> v_GetZIDs(void);
virtual Array<OneD, const unsigned int> v_GetZIDs(void);
virtual ExpListSharedPtr &v_GetPlane(int n)
{
......
......@@ -76,6 +76,21 @@ namespace Nektar
m_outputFrequency = atoi(pParams.find("OutputFrequency")->second.c_str());
}
m_session->MatchSolverInfo("Homogeneous","1D",m_isHomogeneous1D,false);
if(m_isHomogeneous1D)
{
if (pParams.find("OutputPlane") == pParams.end())
{
m_outputPlane = 0;
}
else
{
m_outputPlane = atoi(pParams.find("OutputPlane")->second.c_str());
}
}
ASSERTL0(pParams.find("Points") != pParams.end(),
"Missing parameter 'Points'.");
m_historyPointStream.str(pParams.find("Points")->second);
......@@ -110,6 +125,17 @@ namespace Nektar
while (!m_historyPointStream.fail())
{
m_historyPointStream >> gloCoord[0] >> gloCoord[1] >> gloCoord[2];
if(m_isHomogeneous1D) // overwrite with plane z
{
NekDouble Z = (pFields[0]->GetHomogeneousBasis()->GetZ())[m_outputPlane];
if(fabs(gloCoord[2] - Z) > NekConstants::kVertexTheSameDouble)
{
cout <<"Reseting History point from " << gloCoord[2] <<
" to " << Z << endl;
}
gloCoord[2] = Z;
}
if (!m_historyPointStream.fail())
{
SpatialDomains::VertexComponentSharedPtr vert
......@@ -136,8 +162,30 @@ namespace Nektar
gloCoord[2]);
idList[i] = pFields[0]->GetExpIndex(gloCoord);
if (idList[i] != -1) {
procList[i] = vRank;
if (idList[i] != -1)
{
if(m_isHomogeneous1D)
{
int j;
Array<OneD, const unsigned int> IDs = pFields[0]->GetZIDs();
for(j = 0; j < IDs.num_elements(); ++j)
{
if(IDs[j] == m_outputPlane)
{
break;
}
}
if(j != IDs.num_elements())
{
m_outputPlane == j;
procList[i] = vRank;
}
}
else
{
procList[i] = vRank;
}
}
}
vComm->AllReduce(procList, LibUtilities::ReduceMax);
......@@ -187,7 +235,14 @@ namespace Nektar
m_outputStream << m_session->GetVariable(i) <<",";
}
m_outputStream << ") at points:" << endl;
if(m_isHomogeneous1D)
{
m_outputStream << ") at points:";
}
else
{
m_outputStream << ") at points:" << endl;
}
for (i = 0; i < m_historyPoints.size(); ++i)
{
......@@ -204,6 +259,11 @@ namespace Nektar
m_outputStream << gloCoord[2];
m_outputStream << endl;
}
if(m_isHomogeneous1D)
{
m_outputStream << "(in Wavespace)" << endl;
}
}
v_Update(pFields, time);
}
......@@ -232,15 +292,31 @@ namespace Nektar
// Pull out data values field by field
for (j = 0; j < numFields; ++j)
{
if(pFields[j]->GetPhysState() == false)
if(m_isHomogeneous1D)
{
pFields[j]->BwdTrans(pFields[j]->GetCoeffs(),pFields[j]->UpdatePhys());
if(pFields[j]->GetPhysState() == false)
{
pFields[j]->GetPlane(m_outputPlane)->BwdTrans(pFields[j]->GetPlane(m_outputPlane)->GetCoeffs(),pFields[j]->GetPlane(m_outputPlane)->UpdatePhys());
}
pFields[j]->GetPlane(m_outputPlane)->PutPhysInToElmtExp();
for (k = 0, x = m_historyList.begin(); x != m_historyList.end(); ++x, ++k)
{
(*x).first->GetCoords(gloCoord[0], gloCoord[1], gloCoord[2]);
data[m_historyLocalPointMap[k]*numFields+j] = pFields[j]->GetPlane(m_outputPlane)->GetExp((*x).second)->PhysEvaluate(gloCoord);
}
}
pFields[j]->PutPhysInToElmtExp();
for (k = 0, x = m_historyList.begin(); x != m_historyList.end(); ++x, ++k)
else
{
(*x).first->GetCoords(gloCoord[0], gloCoord[1], gloCoord[2]);
data[m_historyLocalPointMap[k]*numFields+j] = pFields[j]->GetExp((*x).second)->PhysEvaluate(gloCoord);
if(pFields[j]->GetPhysState() == false)
{
pFields[j]->BwdTrans(pFields[j]->GetCoeffs(),pFields[j]->UpdatePhys());
}
pFields[j]->PutPhysInToElmtExp();
for (k = 0, x = m_historyList.begin(); x != m_historyList.end(); ++x, ++k)
{
(*x).first->GetCoords(gloCoord[0], gloCoord[1], gloCoord[2]);
data[m_historyLocalPointMap[k]*numFields+j] = pFields[j]->GetExp((*x).second)->PhysEvaluate(gloCoord);
}
}
}
......
......@@ -74,6 +74,8 @@ namespace Nektar
SpatialDomains::VertexComponentVector m_historyPoints;
unsigned int m_index;
unsigned int m_outputFrequency;
unsigned int m_outputPlane; // plane to take history point from if using a homogeneous1D expansion
bool m_isHomogeneous1D;
std::string m_outputFile;
std::ofstream m_outputStream;
std::stringstream m_historyPointStream;
......
......@@ -1291,7 +1291,7 @@ namespace Nektar
std::vector<std::vector<NekDouble> > &fielddata)
{
ASSERTL1(fielddefs.size() == fielddata.size(),
"Length of fielddefs and fielddata incompatible");
"Length of fielddefs and fielddata incompatible");
TiXmlDocument doc;
TiXmlDeclaration * decl = new TiXmlDeclaration("1.0", "utf-8", "");
......@@ -1307,10 +1307,10 @@ namespace Nektar
ASSERTL1(fielddata[f].size() > 0,
"Fielddata vector must contain at least one value.");
int datasize = CheckFieldDefinition(fielddefs[f]);
ASSERTL1(fielddata[f].size() == fielddefs[f]->m_fields.size()
* datasize, "Invalid size of fielddata vector.");
* datasize, "Invalid size of fielddata vector.");
//---------------------------------------------
// Write ELEMENTS
......
......@@ -55,7 +55,6 @@ int main(int argc, char *argv[])
int datalen2 = fielddata2[i].size()/fielddef2[i]->m_fields.size();
ASSERTL0(datalen1*2 == datalen2,"Data per fields is note compatible");
// Determine the number of coefficients per element
int ncoeffs;
......@@ -143,7 +142,7 @@ int main(int argc, char *argv[])
fielddef2[i]->m_numModes[2] += 2;
fielddef2[i]->m_homogeneousZIDs.push_back(2);
fielddef2[i]->m_homogeneousZIDs.push_back(3);
// check to see if any field in fielddef1[i]->m_fields is
// not defined in fielddef2[i]->m_fields
for(int k = 0; k < fielddef1[i]->m_fields.size(); ++k)
......
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