Commit 645a7313 authored by David Moxey's avatar David Moxey
Browse files

Merge branch 'master' into feature/petsc

parents 45d2c0e8 c73d0600
......@@ -35,6 +35,5 @@ IF( NEKTAR_USE_VTK )
INCLUDE (${VTK_USE_FILE})
MARK_AS_ADVANCED(VTK_DIR)
ADD_DEFINITIONS(-DNEKTAR_USING_VTK)
INCLUDE_DIRECTORIES(${VTK_INCLUDE_DIRS})
ENDIF( NEKTAR_USE_VTK )
......@@ -73,6 +73,7 @@ FOREACH(dir ${dir_list})
ENDFOREACH(dir ${dir_list})
ADD_NEKTAR_TEST(Helmholtz1D_CG_P8)
ADD_NEKTAR_TEST(Helmholtz1D_CG_P8_periodic)
ADD_NEKTAR_TEST(Helmholtz1D_CG_P8_RBC)
ADD_NEKTAR_TEST(Helmholtz1D_HDG_P8)
ADD_NEKTAR_TEST(Helmholtz1D_HDG_P8_RBC)
......@@ -112,6 +113,7 @@ ADD_NEKTAR_TEST(Deriv3D_Homo1D)
ADD_NEKTAR_TEST(Deriv3D_Homo2D)
IF (NEKTAR_USE_MPI)
ADD_NEKTAR_TEST(Helmholtz1D_CG_P8_periodic_par3)
ADD_NEKTAR_TEST(Helmholtz2D_CG_P7_Modes_AllBCs_xxt_full)
ADD_NEKTAR_TEST(Helmholtz2D_CG_P7_Modes_AllBCs_xxt_sc)
ADD_NEKTAR_TEST(Helmholtz2D_CG_P7_Modes_AllBCs_iter_sc_par3)
......
......@@ -47,11 +47,15 @@ int main(int argc, char *argv[])
factors[StdRegions::eFactorLambda] = vSession->GetParameter("Lambda");
const SpatialDomains::ExpansionMap &expansions = graph1D->GetExpansions();
LibUtilities::BasisKey bkey0 = expansions.begin()->second->m_basisKeyVector[0];
cout << "Solving 1D Helmholtz: " << endl;
cout << " Communication: " << vComm->GetType() << endl;
cout << " Solver type : " << vSession->GetSolverInfo("GlobalSysSoln") << endl;
cout << " Lambda : " << factors[StdRegions::eFactorLambda] << endl;
cout << " No. modes : " << bkey0.GetNumModes() << endl;
if (vComm->GetRank() ==0)
{
cout << "Solving 1D Helmholtz: " << endl;
cout << " Communication: " << vComm->GetType() << endl;
cout << " Solver type : " << vSession->GetSolverInfo("GlobalSysSoln") << endl;
cout << " Lambda : " << factors[StdRegions::eFactorLambda] << endl;
cout << " No. modes : " << bkey0.GetNumModes() << endl;
}
//----------------------------------------------
//----------------------------------------------
......
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>Helmholtz 1D CG with P=8, periodic BCs</description>
<executable>Helmholtz1D</executable>
<parameters>Helmholtz1D_P8_periodic.xml</parameters>
<files>
<file description="Session File">Helmholtz1D_P8_periodic.xml</file>
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-9">4.00482e-10</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-9">1.64741e-10</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<test>
<description>Helmholtz 1D CG with P=8, periodic BCs</description>
<executable>Helmholtz1D</executable>
<parameters>Helmholtz1D_P8_periodic.xml</parameters>
<processes>3</processes>
<files>
<file description="Session File">Helmholtz1D_P8_periodic.xml</file>
</files>
<metrics>
<metric type="L2" id="1">
<value tolerance="1e-9">4.00482e-10</value>
</metric>
<metric type="Linf" id="2">
<value tolerance="1e-9">1.64741e-10</value>
</metric>
</metrics>
</test>
<?xml version="1.0" encoding="utf-8"?>
<NEKTAR xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:noNamespaceSchemaLocation="http://www.nektar.info/schema/nektar.xsd">
<GEOMETRY DIM="1" SPACE="2">
<VERTEX>
<V ID="0"> 0 0 0.0 </V>
<V ID="1"> 2 0 0.0 </V>
<V ID="2"> 4 0 0.0 </V>
<V ID="3"> 6 0 0.0 </V>
<V ID="4"> 8 0 0.0 </V>
<V ID="5"> 10 0 0.0 </V>
</VERTEX>
<ELEMENT>
<S ID="0"> 0 1 </S>
<S ID="1"> 1 2 </S>
<S ID="2"> 2 3 </S>
<S ID="3"> 3 4 </S>
<S ID="4"> 4 5 </S>
</ELEMENT>
<COMPOSITE>
<C ID="0"> S[0-4] </C>
<C ID="1"> V[0] </C>
<C ID="2"> V[5] </C>
</COMPOSITE>
<DOMAIN> C[0] </DOMAIN>
</GEOMETRY>
<EXPANSIONS>
<E COMPOSITE="C[0]" NUMMODES="16" FIELDS="u" TYPE="MODIFIED" />
</EXPANSIONS>
<CONDITIONS>
<PARAMETERS>
<P> Lambda = 2 </P>
</PARAMETERS>
<VARIABLES>
<V ID="0"> u </V>
</VARIABLES>
<BOUNDARYREGIONS>
<B ID="0"> C[1] </B>
<B ID="1"> C[2] </B>
</BOUNDARYREGIONS>
<BOUNDARYCONDITIONS>
<REGION REF="0">
<P VAR="u" VALUE="[1]" />
</REGION>
<REGION REF="1">
<P VAR="u" VALUE="[0]" />
</REGION>
</BOUNDARYCONDITIONS>
<FUNCTION NAME="Forcing">
<E VAR="u" VALUE="-(Lambda + PI*PI)*cos(PI*x)" />
</FUNCTION>
<FUNCTION NAME="ExactSolution">
<E VAR="u" VALUE="cos(PI*x)" />
</FUNCTION>
</CONDITIONS>
</NEKTAR>
......@@ -547,6 +547,10 @@ namespace Nektar
weight = StdSegData::getNumberOfCoefficients(na);
bndWeight = StdSegData::getNumberOfBndCoefficients(na);
break;
case 'V':
weight = 1;
bndWeight = 1;
break;
default:
break;
}
......@@ -711,6 +715,10 @@ namespace Nektar
weight = StdSegData::getNumberOfCoefficients(na);
bndWeight = StdSegData::getNumberOfBndCoefficients(na);
break;
case 'V':
weight = 1;
bndWeight = 1;
break;
default:
break;
}
......@@ -1120,6 +1128,12 @@ namespace Nektar
// Based on entity type, check if in this partition
switch (vIt->second.type)
{
case 'V':
if (vVertices.find(vIt->second.list[j]) == vVertices.end())
{
continue;
}
break;
case 'E':
if (vEdges.find(vIt->second.list[j]) == vEdges.end())
{
......
......@@ -53,7 +53,7 @@
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#ifndef MAX_PARAM
#define MAX_PARAM 4 // default maximum number of parameters to support
#define MAX_PARAM 5 // default maximum number of parameters to support
#endif
namespace Nektar
......
......@@ -1373,6 +1373,12 @@ namespace Nektar
"Edge normal not computed.");
return x->second;
}
const StdRegions::NormalVector &Expansion2D::v_GetSurfaceNormal(
const int id) const
{
return v_GetEdgeNormal(id);
}
void Expansion2D::v_NegateEdgeNormal(const int edge)
{
......
......@@ -169,6 +169,7 @@ namespace Nektar
virtual bool v_EdgeNormalNegated(const int edge);
virtual void v_SetUpPhysNormals(const int edge);
const StdRegions::NormalVector &v_GetEdgeNormal(const int edge) const;
const StdRegions::NormalVector &v_GetSurfaceNormal(const int id) const;
};
inline ExpansionSharedPtr Expansion2D::GetEdgeExp(int edge, bool SetUpNormal)
......
......@@ -214,24 +214,39 @@ namespace Nektar
}
}
// STEP 2: Order the periodic vertices next
// This allows to give corresponding DOF's the same
// global ID
// STEP 2: Order the periodic vertices next. This allows to give
// corresponding DOF's the same global ID
for(pIt = periodicVerts.begin(); pIt != periodicVerts.end(); pIt++)
{
meshVertId = pIt->first;
meshVertId2 = pIt->second[0].id;
ASSERTL0(vertReorderedGraphVertId.count(meshVertId) == 0,
"This periodic boundary vertex has been specified before");
ASSERTL0(vertReorderedGraphVertId.count(meshVertId) == 0,
"This periodic boundary vertex has been specified before");
vertReorderedGraphVertId[meshVertId] = graphVertId;
vertReorderedGraphVertId[meshVertId2] = graphVertId++;
if (!pIt->second[0].isLocal)
{
if (vertReorderedGraphVertId.count(meshVertId) == 0)
{
vertReorderedGraphVertId[meshVertId] = graphVertId++;
}
}
else
{
if (vertReorderedGraphVertId.count(meshVertId) == 0 &&
vertReorderedGraphVertId.count(meshVertId2) == 0)
{
vertReorderedGraphVertId[meshVertId] = graphVertId;
vertReorderedGraphVertId[meshVertId2] = graphVertId++;
}
else if ((vertReorderedGraphVertId.count(meshVertId) == 0 &&
vertReorderedGraphVertId.count(meshVertId2) != 0) ||
(vertReorderedGraphVertId.count(meshVertId) != 0 &&
vertReorderedGraphVertId.count(meshVertId2) == 0))
{
ASSERTL0(false,
"Numbering error for periodic boundary conditions!");
}
}
}
// STEP 3: List the other vertices
// STEP 4: Set up simple map based on vertex and edge id's
for(i = 0; i < locExpVector.size(); ++i)
......
......@@ -546,18 +546,17 @@ namespace Nektar
SpatialDomains::MeshGraph1D>(m_graph);
SpatialDomains::BoundaryRegionCollection::const_iterator it;
LibUtilities::CommSharedPtr vComm =
LibUtilities::CommSharedPtr vComm =
m_session->GetComm()->GetRowComm();
int region1ID;
int region2ID;
int i, region1ID, region2ID;
SpatialDomains::BoundaryConditionShPtr locBCond;
map<int,int> BregionToVertMap;
map<int,int> BregionToVertMap;
// Construct list of all periodic Region and their global vertex on
// this process.
// Construct list of all periodic Region and their global vertex on
// this process.
for (it = bregions.begin(); it != bregions.end(); ++it)
{
locBCond = GetBoundaryCondition(bconditions, it->first, variable);
......@@ -575,44 +574,40 @@ namespace Nektar
map<int,int>::iterator iit;
set<int> islocal;
if(vComm->GetSize() != 1) // share map between all processors
{
// number of bregions
int nregions = 0;
for(iit = BregionToVertMap.begin(); iit != BregionToVertMap.end(); ++iit)
{
nregions = (nregions > iit->first)? nregions:iit->first;
}
vComm->AllReduce(nregions, LibUtilities::ReduceMax);
Array<OneD, int> bregmap(nregions,-1);
for(iit = BregionToVertMap.begin(); iit != BregionToVertMap.end(); ++iit)
{
bregmap[iit->first] = iit->second;
islocal.insert(iit->first);
}
int n = vComm->GetSize();
int p = vComm->GetRank();
vComm->AllReduce(bregmap, LibUtilities::ReduceSum);
Array<OneD, int> nregions(n, 0);
nregions[p] = BregionToVertMap.size();
vComm->AllReduce(nregions, LibUtilities::ReduceSum);
for(int i = 0; i < nregions; ++i)
{
if(bregmap[i] >=0 )
{
BregionToVertMap[i] = bregmap[i];
}
}
int totRegions = Vmath::Vsum(n, nregions, 1);
Array<OneD, int> regOffset(n, 0);
for (i = 1; i < n; ++i)
{
regOffset[i] = regOffset[i-1] + nregions[i-1];
}
else
Array<OneD, int> bregmap(totRegions, 0);
Array<OneD, int> bregid (totRegions, 0);
for(i = regOffset[p], iit = BregionToVertMap.begin();
iit != BregionToVertMap.end(); ++iit, ++i)
{
// in serial case set all verties on boundary as being local
for(iit = BregionToVertMap.begin(); iit != BregionToVertMap.end(); ++iit)
{
islocal.insert(iit->first);
}
bregid [i] = iit->first;
bregmap[i] = iit->second;
islocal.insert(iit->first);
}
vComm->AllReduce(bregmap, LibUtilities::ReduceSum);
vComm->AllReduce(bregid, LibUtilities::ReduceSum);
for (int i = 0; i < totRegions; ++i)
{
BregionToVertMap[bregid[i]] = bregmap[i];
}
// Construct list of all periodic pairs local to this process.
for (it = bregions.begin(); it != bregions.end(); ++it)
......@@ -640,7 +635,7 @@ namespace Nektar
PeriodicEntity ent(BregionToVertMap[region2ID],
StdRegions::eNoOrientation,
islocal.count(region2ID) != 0);
m_periodicVerts[BregionToVertMap[region1ID]].push_back(ent);
}
}
......
......@@ -1255,8 +1255,8 @@ namespace Nektar
// Find nearest element
if (!elmtIdDist.empty())
{
NekDouble min_d = elmtIdDist[0].second;
int min_id = elmtIdDist[0].first;
NekDouble min_d = elmtIdDist[0].first;
int min_id = elmtIdDist[0].second;
for (int i = 1; i < elmtIdDist.size(); ++i)
{
......@@ -1266,10 +1266,6 @@ namespace Nektar
}
}
// retrieve local coordinates of chosen point
(*m_exp)[min_id]->GetGeom()->ContainsPoint(gloCoords,
locCoords,
tol, resid);
return min_id;
}
else {
......@@ -1313,24 +1309,6 @@ namespace Nektar
}
/**
* The operation is evaluated locally by the elemental
* function StdRegions#StdExpansion#GetSurfaceNormal.
*/
void ExpList::GetSurfaceNormal(Array<OneD, NekDouble> &SurfaceNormal,
const int k)
{
int i;
Array<OneD, NekDouble> normals;
for(i = 0; i < (*m_exp).size(); ++i)
{
//e_SN = SurfaceNormal + m_phys_offset[i];
normals = (*m_exp)[i]->GetSurfaceNormal()[k];
Vmath::Vcopy(normals.num_elements(), &normals[0], 1, &SurfaceNormal[0] + m_phys_offset[i], 1);
}
}
/**
* Configures geometric info, such as tangent direction, on each
* expansion.
......
......@@ -347,12 +347,6 @@ namespace Nektar
Array<OneD, NekDouble> &outarray,
int BndID);
/// This function calculates Surface Normal vector of a smooth
/// manifold.
MULTI_REGIONS_EXPORT void GetSurfaceNormal(
Array<OneD,NekDouble> &SurfaceNormal,
const int k);
/// Apply geometry information to each expansion.
MULTI_REGIONS_EXPORT void ApplyGeomInfo();
......
......@@ -70,9 +70,10 @@ namespace Nektar
void Forcing::Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr>& fields,
const Array<OneD, Array<OneD, NekDouble> >& inarray,
Array<OneD, Array<OneD, NekDouble> >& outarray)
Array<OneD, Array<OneD, NekDouble> >& outarray,
const NekDouble& time)
{
v_Apply(fields, inarray, outarray);
v_Apply(fields, inarray, outarray, time);
}
......@@ -114,6 +115,26 @@ namespace Nektar
return vForceList;
}
void Forcing::EvaluateTimeFunction(
LibUtilities::SessionReaderSharedPtr pSession,
std::string pFieldName,
Array<OneD, NekDouble>& pArray,
const std::string& pFunctionName,
NekDouble pTime)
{
ASSERTL0(pSession->DefinesFunction(pFunctionName),
"Function '" + pFunctionName + "' does not exist.");
LibUtilities::EquationSharedPtr ffunc =
pSession->GetFunction(pFunctionName, pFieldName);
Array<OneD, NekDouble> x0(1,0.0);
Array<OneD, NekDouble> x1(1,0.0);
Array<OneD, NekDouble> x2(1,0.0);
ffunc->Evaluate(x0, x1, x2, pTime, pArray);
}
void Forcing::EvaluateFunction(
Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
......@@ -139,11 +160,11 @@ namespace Nektar
Array<OneD, NekDouble> x0(nq);
Array<OneD, NekDouble> x1(nq);
Array<OneD, NekDouble> x2(nq);
pFields[0]->GetCoords(x0, x1, x2);
LibUtilities::EquationSharedPtr ffunc =
pSession->GetFunction(pFunctionName, pFieldName);
pSession->GetFunction(pFunctionName, pFieldName);
ffunc->Evaluate(x0, x1, x2, pTime, pArray);
}
else if (vType == LibUtilities::eFunctionTypeFile)
......
......@@ -80,7 +80,8 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT void Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
SOLVER_UTILS_EXPORT static std::vector<ForcingSharedPtr> Load(
const LibUtilities::SessionReaderSharedPtr& pSession,
......@@ -107,15 +108,24 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)=0;
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)=0;
void EvaluateFunction(
Array<OneD, MultiRegions::ExpListSharedPtr> pFields,
LibUtilities::SessionReaderSharedPtr pSession,
std::string pFieldName, Array<OneD, NekDouble>& pArray,
LibUtilities::SessionReaderSharedPtr pSession,
std::string pFieldName,
Array<OneD, NekDouble>& pArray,
const std::string& pFunctionName,
NekDouble pTime = NekDouble(0));
void EvaluateTimeFunction(
LibUtilities::SessionReaderSharedPtr pSession,
std::string pFieldName,
Array<OneD, NekDouble>& pArray,
const std::string& pFunctionName,
NekDouble pTime = NekDouble(0));
};
}
}
......
......@@ -110,7 +110,8 @@ namespace SolverUtils
void ForcingBody::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
for (int i = 0; i < m_NumVariable; i++)
{
......
......@@ -79,7 +79,8 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
private:
ForcingBody(const LibUtilities::SessionReaderSharedPtr& pSession);
......
......@@ -80,7 +80,8 @@ namespace SolverUtils
void ForcingNoise::v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray)
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time)
{
for (int i = 0; i < m_NumVariable; i++)
{
......
......@@ -79,7 +79,8 @@ namespace SolverUtils
SOLVER_UTILS_EXPORT virtual void v_Apply(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
Array<OneD, Array<OneD, NekDouble> > &outarray);
Array<OneD, Array<OneD, NekDouble> > &outarray,
const NekDouble &time);
private:
ForcingNoise(const LibUtilities::SessionReaderSharedPtr& pSession);
......