Commit a5b6b4d2 authored by Dave Moxey's avatar Dave Moxey
Browse files

Merge branch 'fix/FC3DH1Defficiency' into 'master'

Switch to Fourier points in FieldConvert with homogeneous expansions

This MR switches the point distribution of homogeneous expansions in FieldConvert to FourierEvenlySpaced.

The main reason for that is a performance issue recently reported when the number of Fourier modes is large. This was caused by the expansion creating a PolyEvenlySpaced distribution with many points, which scales very poorly due to the derivative matrix (for 1000 points this takes around 2 hours). This cost is also independent of the number of mpi processes...

This is also related to issue #10.

See merge request !627
parents fda40e6f 3a3fabcd
......@@ -15,6 +15,7 @@ v4.3.2
**FieldConvert**:
- Fix appearence of duplicate messages when running in parallel (!626)
- Fix issue with efficiency when using large number of 3DH1D planes (!627)
**Packaging**:
- Fixes for DEB package dependencies (!630)
......
......@@ -302,12 +302,12 @@ namespace Nektar
int retval = MPI_Sendrecv(pSendData.get(),
(int) pSendData.num_elements(),
MPI_DOUBLE,
pRecvProc,
pSendProc,
0,
pRecvData.get(),
(int) pRecvData.num_elements(),
MPI_DOUBLE,
pSendProc,
pRecvProc,
0,
m_comm,
&status);
......
......@@ -1217,7 +1217,7 @@ namespace Nektar
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
outfile << " <Points>" << endl;
outfile << " <DataArray type=\"Float32\" "
outfile << " <DataArray type=\"Float64\" "
<< "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
outfile << " ";
for (i = 0; i < ntot; ++i)
......
......@@ -367,9 +367,24 @@ namespace Nektar
int expansion,
int istrip)
{
// If there is only one plane (e.g. HalfMode), we write a 2D plane.
if (m_planes.num_elements() == 1)
{
m_planes[0]->WriteVtkPieceHeader(outfile, expansion);
return;
}
// If we are using Fourier points, output extra plane to fill domain
int outputExtraPlane = 0;
if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
&& m_homogeneousBasis->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
outputExtraPlane = 1;
}
int i, j;
int nquad0 = (*m_exp)[expansion]->GetNumPoints(0);
int nquad1 = m_planes.num_elements();
int nquad1 = m_planes.num_elements() + outputExtraPlane;
int ntot = nquad0*nquad1;
int ntotminus = (nquad0-1)*(nquad1-1);
......@@ -379,6 +394,20 @@ namespace Nektar
coords[2] = Array<OneD,NekDouble>(ntot);
GetCoords(expansion,coords[0],coords[1],coords[2]);
if (outputExtraPlane)
{
// Copy coords[0] and coords[1] to extra plane
Array<OneD,NekDouble> tmp;
Vmath::Vcopy (nquad0, coords[0], 1,
tmp = coords[0] + (nquad1-1)*nquad0, 1);
Vmath::Vcopy (nquad0, coords[1], 1,
tmp = coords[1] + (nquad1-1)*nquad0, 1);
// Fill coords[2] for extra plane
NekDouble z = coords[2][nquad0*m_planes.num_elements()-1] +
(coords[2][nquad0] - coords[2][0]);
Vmath::Fill(nquad0, z, tmp = coords[2] + (nquad1-1)*nquad0, 1);
}
outfile << " <Piece NumberOfPoints=\""
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
......
......@@ -109,7 +109,7 @@ namespace Nektar
(*m_exp).push_back(m_planes[0]->GetExp(j));
}
for(n = 1; n < m_homogeneousBasis->GetNumPoints(); ++n)
for(n = 1; n < m_planes.num_elements(); ++n)
{
m_planes[n] = MemoryManager<ExpList2D>::AllocateSharedPtr(*plane_zero,False);
for(j = 0; j < nel; ++j)
......@@ -371,10 +371,19 @@ namespace Nektar
return;
}
// If we are using Fourier points, output extra plane to fill domain
int outputExtraPlane = 0;
if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
&& m_homogeneousBasis->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
outputExtraPlane = 1;
}
int i,j,k;
int nq0 = (*m_exp)[expansion]->GetNumPoints(0);
int nq1 = (*m_exp)[expansion]->GetNumPoints(1);
int nq2 = m_planes.num_elements();
int nq2 = m_planes.num_elements() + outputExtraPlane;
int ntot = nq0*nq1*nq2;
int ntotminus = (nq0-1)*(nq1-1)*(nq2-1);
......@@ -384,6 +393,20 @@ namespace Nektar
coords[2] = Array<OneD,NekDouble>(ntot);
GetCoords(expansion,coords[0],coords[1],coords[2]);
if (outputExtraPlane)
{
// Copy coords[0] and coords[1] to extra plane
Array<OneD,NekDouble> tmp;
Vmath::Vcopy (nq0*nq1, coords[0], 1,
tmp = coords[0] + (nq2-1)*nq0*nq1, 1);
Vmath::Vcopy (nq0*nq1, coords[1], 1,
tmp = coords[1] + (nq2-1)*nq0*nq1, 1);
// Fill coords[2] for extra plane
NekDouble z = coords[2][nq0*nq1*m_planes.num_elements()-1] +
(coords[2][nq0*nq1] - coords[2][0]);
Vmath::Fill(nq0*nq1, z, tmp = coords[2] + (nq2-1)*nq0*nq1, 1);
}
NekDouble DistStrip;
m_session->LoadParameter("DistStrip", DistStrip, 0);
// Reset the z-coords for homostrips
......@@ -396,7 +419,7 @@ namespace Nektar
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
outfile << " <Points>" << endl;
outfile << " <DataArray type=\"Float32\" "
outfile << " <DataArray type=\"Float64\" "
<< "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
outfile << " ";
for (i = 0; i < ntot; ++i)
......
......@@ -367,7 +367,7 @@ namespace Nektar
<< ntot << "\" NumberOfCells=\""
<< ntotminus << "\">" << endl;
outfile << " <Points>" << endl;
outfile << " <DataArray type=\"Float32\" "
outfile << " <DataArray type=\"Float64\" "
<< "NumberOfComponents=\"3\" format=\"ascii\">" << endl;
outfile << " ";
for (i = 0; i < ntot; ++i)
......
......@@ -872,8 +872,37 @@ namespace Nektar
int nq = (*m_exp)[expansion]->GetTotPoints();
int npoints_per_plane = m_planes[0]->GetTotPoints();
// If we are using Fourier points, output extra plane to fill domain
int outputExtraPlane = 0;
Array<OneD, NekDouble> extraPlane;
if ( m_homogeneousBasis->GetBasisType() == LibUtilities::eFourier
&& m_homogeneousBasis->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
outputExtraPlane = 1;
// Get extra plane data
if (m_StripZcomm->GetSize() == 1)
{
extraPlane = m_phys + m_phys_offset[expansion];
}
else
{
// Determine to and from rank for communication
int size = m_StripZcomm->GetSize();
int rank = m_StripZcomm->GetRank();
int fromRank = (rank+1) % size;
int toRank = (rank == 0) ? size-1 : rank-1;
// Communicate using SendRecv
extraPlane = Array<OneD, NekDouble>(nq);
Array<OneD, NekDouble> send (nq,
m_phys + m_phys_offset[expansion]);
m_StripZcomm->SendRecv(toRank, send,
fromRank, extraPlane);
}
}
// printing the fields of that zone
outfile << " <DataArray type=\"Float32\" Name=\""
outfile << " <DataArray type=\"Float64\" Name=\""
<< var << "\">" << endl;
outfile << " ";
for (int n = 0; n < m_planes.num_elements(); ++n)
......@@ -884,6 +913,14 @@ namespace Nektar
outfile << (fabs(phys[i]) < NekConstants::kNekZeroTol ? 0 : phys[i]) << " ";
}
}
if (outputExtraPlane)
{
for(i = 0; i < nq; ++i)
{
outfile << (fabs(extraPlane[i]) < NekConstants::kNekZeroTol ?
0 : extraPlane[i]) << " ";
}
}
outfile << endl;
outfile << " </DataArray>" << endl;
}
......
......@@ -734,7 +734,7 @@ namespace Nektar
int npoints_per_line = m_lines[0]->GetTotPoints();
// printing the fields of that zone
outfile << " <DataArray type=\"Float32\" Name=\""
outfile << " <DataArray type=\"Float64\" Name=\""
<< var << "\">" << endl;
outfile << " ";
for (int n = 0; n < m_lines.num_elements(); ++n)
......
......@@ -156,7 +156,7 @@ struct Field {
// Choose points to be at evenly spaced points at
// nplanes points
const LibUtilities::PointsKey
Pkey(nplanes, LibUtilities::ePolyEvenlySpaced);
Pkey(nplanes, LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey Bkey(btype, nplanes, Pkey);
......@@ -207,11 +207,11 @@ struct Field {
// Choose points to be at evenly spaced points at
// nplanes points
const LibUtilities::PointsKey
PkeyY(nylines, LibUtilities::ePolyEvenlySpaced);
PkeyY(nylines, LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey BkeyY(btype1, nylines, PkeyY);
const LibUtilities::PointsKey
PkeyZ(nzlines, LibUtilities::ePolyEvenlySpaced);
PkeyZ(nzlines, LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey BkeyZ(btype2, nzlines, PkeyZ);
if(m_declareExpansionAsContField)
......@@ -309,7 +309,7 @@ struct Field {
// Choose points to be at evenly spaced points at
// nplanes points
const LibUtilities::PointsKey
Pkey(nplanes, LibUtilities::ePolyEvenlySpaced);
Pkey(nplanes, LibUtilities::eFourierEvenlySpaced);
const LibUtilities::BasisKey Bkey(btype, nplanes, Pkey);
......
......@@ -68,6 +68,9 @@ int main(int argc, char* argv[])
("nprocs", po::value<int>(),
"Used to define nprocs if running serial problem to mimic "
"parallel run.")
("npz", po::value<int>(),
"Used to define number of partitions in z for Homogeneous1D "
"expansions for parallel runs.")
("onlyshape", po::value<string>(),
"Only use element with defined shape type i.e. -onlyshape "
" Tetrahedron")
......
......@@ -253,6 +253,13 @@ void InputXml::Process(po::variables_map &vm)
boost::lexical_cast<string>(vm["part-only-overlapping"].as<int>()));
}
if(vm.count("npz"))
{
cmdArgs.push_back("--npz");
cmdArgs.push_back(
boost::lexical_cast<string>(vm["npz"].as<int>()));
}
int argc = cmdArgs.size();
const char **argv = new const char*[argc];
for (int i = 0; i < argc; ++i)
......
......@@ -463,6 +463,76 @@ void ProcessEquiSpacedOutput::SetHomogeneousConnectivity(void)
int npts = m_f->m_fieldPts->GetNpoints();
int nptsPerPlane = npts/nPlanes;
if ( m_f->m_exp[0]->GetHomogeneousBasis()->GetBasisType() ==
LibUtilities::eFourier
&& m_f->m_exp[0]->GetHomogeneousBasis()->GetPointsType() ==
LibUtilities::eFourierEvenlySpaced)
{
// Write points with extra plane
Array<OneD, Array<OneD, NekDouble> > pts;
m_f->m_fieldPts->GetPts(pts);
Array<OneD, Array<OneD, NekDouble> > newPts(pts.num_elements());
for (int i = 0; i < pts.num_elements(); i++)
{
newPts[i] = Array<OneD, NekDouble> (npts + nptsPerPlane);
// Copy old points
Vmath::Vcopy(npts, pts[i], 1, newPts[i], 1);
// Get new plane
Array<OneD, NekDouble> extraPlane (nptsPerPlane, newPts[i]+npts);
if (m_f->m_comm->GetColumnComm()->GetSize() == 1)
{
if ( i == 2)
{
// z-coordinate
Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
pts[i], 1, extraPlane, 1);
}
else
{
Vmath::Vcopy(nptsPerPlane, pts[i], 1, extraPlane, 1);
}
}
else
{
// Determine to and from rank for communication
int size = m_f->m_comm->GetColumnComm()->GetSize();
int rank = m_f->m_comm->GetColumnComm()->GetRank();
int fromRank = (rank+1) % size;
int toRank = (rank == 0) ? size-1 : rank-1;
// Communicate using SendRecv
Array<OneD, NekDouble> send (nptsPerPlane, pts[i]);
m_f->m_comm->GetColumnComm()->SendRecv(toRank, send,
fromRank, extraPlane);
// Adjust z-coordinate of last process
if (i == 2 && (rank == size - 1) )
{
Vmath::Sadd(nptsPerPlane, m_f->m_exp[0]->GetHomoLen(),
extraPlane, 1, extraPlane, 1);
}
}
}
m_f->m_fieldPts->SetPts(newPts);
// Now extend plane connectivity
vector<Array<OneD, int> > oldConn;
Array<OneD, int> conn;
m_f->m_fieldPts->GetConnectivity(oldConn);
vector<Array<OneD, int> > newConn = oldConn;
int connPerPlane = oldConn.size()/nPlanes;
for (int i = 0; i < connPerPlane; i++)
{
conn = Array<OneD, int> (oldConn[i].num_elements());
for (int j = 0; j < conn.num_elements(); j++)
{
conn[j] = oldConn[i][j] + npts;
}
newConn.push_back(conn);
}
m_f->m_fieldPts->SetConnectivity(newConn);
nPlanes++;
npts += nptsPerPlane;
}
// Get maximum number of points per direction in the mesh
int maxN = 0;
for(int i = 0; i < nel; ++i)
......
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